Logistic回归实例2

2744557306

# logistic回归
0 ?& x: L6 d- f& B& v实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。
) u0 s% q% ^5 @6 [+ i
+ a' D9 ?* r1 K$ w: JR软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。
; o: b1 Q  |& g/ G9 L0 u; R
/ s( G9 u8 J. m8 `! v1 _/ v
有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。

1）正态分布 正态分布族的使用方法： fm <- glm(formula, family=gaussian(link=identity), data=data.frame) 其中，link=identity可以不写，因为正态分布的连接函数缺省值是恒等（identity）。事实上，整个参数family=gaussian也可以不写，因为分布族的缺省值就是正态分布。 注意：正态分布的广义线性模型实际上与线性模型是相同的，也就是 fm <- glm(formula, family=gaussian, data=data.frame) 与线性模型 fm <- lm(formula, data=data.frame)有完全相同的结果，但效率却低得多。

2）二项分布

3 S9 r4 Z! Y6 v+ c
logistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。
1 ~( a0 y8 k6 q
Logistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据
   - ?+ d$ t+ f( J% A9 q1 w0 ]
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )8 @( G; s\" p0 B
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  ; i( \) S- |; A9 \2 z
   
4. , e- c$ a  t$ z# X; f+ z
   
5. #2、建模
   ' u6 K4 D( p; N- x! Z
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)' f# B7 D4 \# [# q/ f
   
7. 
   * k& O( s8 J( ]! H
8. #3、模型评估
   ' K! ^- |: }  F5 \' m; J\" G/ C3 o
9. summary(glm.sol)

复制代码

1. ##
   1 I2 A1 K2 b0 @8 \: W5 h
2. ## Call:8 ]) g\" M6 w* f; V' Y$ n
   
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)
   5 z: \% h$ f7 P8 y( I
4. ##
   2 `  }: z; p4 k
5. ## Deviance Residuals:
   3 K3 ], U. Z% \5 q/ X- m4 d; ~: b
6. ##       1        2        3        4        5        6  
   ( F3 J2 C! O9 j: e  e4 o
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  + t% X/ E; t+ B\" C& H3 H
   
8. ##
   9 Y\" r3 D, h$ e) T
9. ## Coefficients:
   . ^) z, R8 T5 s, q7 K
10. ##             Estimate Std. Error z value Pr(>|z|)   
    % a* E2 H( o. m\" g\" S; b
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***8 ]. s: B0 M7 t% i0 X( q
    
12. ## x             1.2459     0.1119   11.13   <2e-16 ***
    / U0 y+ P7 T7 e0 v+ u; c* Y4 l
13. ## ---
    5 }/ U0 Q' y8 K2 |4 t. u) d
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12 r' M3 E- f5 D$ y
    
15. ##
    \" B& o; E0 H. [! ]
16. ## (Dispersion parameter for binomial family taken to be 1)# F% N  @1 u5 y: G% J& d; K( k
    
17. ## \" C9 L; _1 I0 f, C* t) b! j
    
18. ##     Null deviance: 250.4866  on 5  degrees of freedom
    ) m* C1 B% I  @7 Z( `- p
19. ## Residual deviance:   9.3526  on 4  degrees of freedom6 ~9 l7 _# ?, `) ?  T6 y
    
20. ## AIC: 34.093
    : V# z* Z/ D! \4 H# ]% m
21. ##
    6 g  z: `$ U* W! t! P8 M) w1 i* Y$ b$ ^
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率
   9 p\" n2 `/ U/ u\" M8 {
2. % ]6 \. k6 n3 i7 w. X& x% |5 n, |
   
3. #4、预测2 x! L( m( q1 _/ V4 |+ g, E. O
   
4. pre <- predict(glm.sol, data.frame(x=3.5))
   1 m- m- e5 d$ O
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1 ; s) ~1 O1 q' X% l7 f% c. V
   
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1
   0 K. [2 `3 z, W9 k
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x 2 S4 Z1 }3 d9 N
   
2. ##   -3.301035    1.245937

复制代码

1. (X <- -glm.sol$coefficients[[1]]/glm.sol$coefficients[[2]])

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：1 {% O0 ], d/ {- ]1 l
   
2. d <- seq(0, 5, length=100)
   & C* |( q\" Q3 k, |. N6 y
3. pre <- predict(glm.sol, data.frame(x=d))
   + }\" {! v/ b8 Q* ~
4. p <- exp(pre)/(1+exp(pre))
   & M( j+ \4 g! h- g) b3 n5 Z5 A, q7 ?
5. norell$y <- norell$success/norell$n/ X$ n9 ~! q- s  p$ R* h
   
6. plot(norell$x, norell$y)! t3 X  D. i\" L
   
7. lines(d, p)

复制代码

1. #其中，d是给出曲线横坐标的点，pre是计算预测值，p是相应的预测概率。用plot函数和lines给出散点图和对应的预测曲线。

复制代码