Logistic回归实例2

2744557306

# logistic回归
实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。

R软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。
: ?+ }' j8 r' D+ H; p4 S5 |

有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。

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)有完全相同的结果，但效率却低得多。

) M) {6 Q  r/ N8 v) |" @& i) ~2）二项分布
3 L3 v2 t  Q1 r! {  l  l4 v

logistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。

' E  u. c7 g2 d# }4 i; \5 xLogistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
# Y7 L* o1 k4 \) c, }& z! l9 N& N实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据
   ) c9 z3 w! h& {6 N+ U9 c
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )
   6 n) z1 y: L) O2 Q$ j
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  
   + k2 b0 d& F) P2 j
4. . s3 o% b6 g/ x% n/ S7 a! h) z% R
   
5. #2、建模4 p7 ~% H$ P8 T/ t
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)9 `; s1 {  ]; `  e9 [
   
7. 
   0 w4 y2 Y% ?; ^2 M2 K$ K
8. #3、模型评估7 z2 b$ v  C! S& q2 r9 X/ d
   
9. summary(glm.sol)

复制代码

1. ##
   5 _; V6 ], A! y2 t: R3 k
2. ## Call:
   \" R( l* m% Y. k1 y
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)2 _4 X! ?! W5 Z  A0 ]  V6 k# H& Z
   
4. ##
   * C5 @9 q6 ]+ T5 b
5. ## Deviance Residuals: , a' L# f: z9 S7 v0 ~5 @
   
6. ##       1        2        3        4        5        6  : A! q$ v7 A0 j3 U
   
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  8 D- I, M% p6 O$ R
   
8. ## 3 `: K' P$ G( M) [5 k
   
9. ## Coefficients:$ e$ r# H0 @# W( f- [5 m; w0 s
   
10. ##             Estimate Std. Error z value Pr(>|z|)    1 g8 A& V\" u0 _; i5 U# o  t# n
    
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***1 I% R6 G\" d$ ]7 z7 D8 ]
    
12. ## x             1.2459     0.1119   11.13   <2e-16 ***
    ! u* n' G% B- ]- U* O1 d3 P
13. ## ---
    6 J. g5 }8 ^# u. p4 O1 V! Q
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1. g; k5 J- l3 x7 J( `5 [5 O
    
15. ##   I% h' v& F; g6 ^9 p' c7 O6 c7 R
    
16. ## (Dispersion parameter for binomial family taken to be 1)
    3 x2 D# X1 u& p$ {  Y
17. ##
    ( ~. z\" r& o. Y
18. ##     Null deviance: 250.4866  on 5  degrees of freedom* F8 B$ y\" O! D7 t2 ~5 }3 m+ Z
    
19. ## Residual deviance:   9.3526  on 4  degrees of freedom6 t0 Q( U7 C# X\" W' c* E! `: `+ X
    
20. ## AIC: 34.093( K\" O; [- G( R! y! x9 I' f
    
21. ## : m2 }$ O% G* o8 d( G
    
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率; Z- E& C& w; e# N% F
   
2. 
   ; |) i. H( s1 {3 B
3. #4、预测
   ; m- R; S' b. S
4. pre <- predict(glm.sol, data.frame(x=3.5))6 l\" g! v& a% [/ a, t& F4 H  z
   
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1 $ V/ v3 b+ M$ a8 U7 I
   
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1
   % D0 E) L* j1 Q5 X0 l
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x # ?1 I: m- V: @, |% t
   
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：' @- L% @- Q* W- w! u7 E& S
   
2. d <- seq(0, 5, length=100)
   : u2 E6 x0 N# _$ {) k' N
3. pre <- predict(glm.sol, data.frame(x=d))
     Q4 i) q  z\" g6 H7 x4 v\" |/ M5 S7 O4 l
4. p <- exp(pre)/(1+exp(pre))
   6 D9 B5 e/ Q7 b8 g9 [% |
5. norell$y <- norell$success/norell$n! o4 E1 m+ t; \\" U( d' \# z
   
6. plot(norell$x, norell$y)8 P* W* j7 ~# t4 l7 R
   
7. lines(d, p)

复制代码

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

复制代码