Logistic回归实例2

2744557306

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

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


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

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)有完全相同的结果，但效率却低得多。
0 z9 s/ g5 w' E1 f4 u  J
! v2 r4 q' Z# y# V. T4 u2）二项分布


( w0 G1 p' Q9 O2 C3 _: Q6 ?$ Hlogistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。
, `, W+ T! p2 k3 S+ g% `! t
Logistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
5 G' ~8 w& a6 s实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据
   ! r8 Z4 F+ p; }0 h% \0 e% |) F( f
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )& o  E) E) j- N
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  
   1 m2 D1 u& c& X& p\" @( K
4. 
   $ e9 M! o. m\" W
5. #2、建模, V  l; S1 S, \
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)4 ?: q! t$ k5 T0 ?2 G4 I9 S
   
7. 
   ( g( S\" c, U8 j  a, a- J; j! s
8. #3、模型评估
   0 L2 K% Q* B\" K& E3 G; ?1 H- h
9. summary(glm.sol)

复制代码

1. ##
   : O: _0 X8 p3 a0 c( x
2. ## Call:' U: @! a2 G& Y\" H
   
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)\" v3 |# Q2 p2 u7 E
   
4. ##
   4 X1 ^& o5 x1 T3 {
5. ## Deviance Residuals: ( `3 U: ~6 x, F( e: o4 b
   
6. ##       1        2        3        4        5        6  
   ' a) p& q0 Z5 |
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
   ! F+ o$ U0 t! \) w# G  {3 G$ B7 |
8. ##
   5 o  n' }6 V7 ^) {9 [; i8 \
9. ## Coefficients:: F9 S7 w6 @; F5 G\" F; s& q7 b
   
10. ##             Estimate Std. Error z value Pr(>|z|)    / j. T+ |0 ]+ E6 E/ a* E, G; S
    
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***# l7 ^8 D8 g; q9 H' v* p. Y7 b
    
12. ## x             1.2459     0.1119   11.13   <2e-16 ***# [4 c2 @. ]5 C8 T
    
13. ## ---0 v/ {  P$ p# P8 E, j
    
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1* c4 [' P# C, t
    
15. ## % V. D( q' R) W' M( f0 `, e- E
    
16. ## (Dispersion parameter for binomial family taken to be 1). v! R! \\" }) i8 `* G9 k
    
17. ## ) w! [8 J* r0 x* l, S# i: k2 G4 T
    
18. ##     Null deviance: 250.4866  on 5  degrees of freedom
    & @& X1 Q& B8 ~: z* ^% E
19. ## Residual deviance:   9.3526  on 4  degrees of freedom
    2 r0 i7 u9 h: N  `1 A9 T
20. ## AIC: 34.093- ~( M/ h3 U8 j& b7 J' j% m
    
21. ## , A* @/ Q- ]; G9 m
    
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率
   % g! [! M/ c4 X4 u
2. 
   5 T$ u# t, m; T' B
3. #4、预测7 A. s) u2 t8 z! r
   
4. pre <- predict(glm.sol, data.frame(x=3.5))
   ) A% G- e7 [: ?( O2 Y
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1 ( S# o: i+ @/ K; D6 @) D2 `
   
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1
   2 e* ]  F; n  V, c' t& T8 l
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x
   0 e! u3 p2 A4 P# z3 B
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：, \+ ~  b  F, n2 `% r1 H
   
2. d <- seq(0, 5, length=100)$ y' o$ v, u. e* j, _+ z9 B
   
3. pre <- predict(glm.sol, data.frame(x=d))$ p& q3 c; G+ z. U9 a
   
4. p <- exp(pre)/(1+exp(pre))( T' G9 c( ~& D$ C6 V/ L2 a
   
5. norell$y <- norell$success/norell$n+ m4 P4 H' c$ j6 p; y2 f
   
6. plot(norell$x, norell$y)
   1 W3 U( o& w9 R0 l: @
7. lines(d, p)

复制代码

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

复制代码

& e$ H, e" R' v0 N