Logistic回归实例2

2744557306

# logistic回归
实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。
1 c1 e; u) @* V$ g% L2 t
' C4 t- w. H, X; B! o* ]$ ^R软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。

- G6 o! Q4 F( Q
有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。
$ R0 Y- i0 z/ @/ H
3 d/ ~( ]9 c! o8 i3 ~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)有完全相同的结果，但效率却低得多。

! s6 ?; K+ A- Y) a# Y5 T* g4 w% {' Q7 Q2）二项分布
$ Z3 C  {( D9 h: U, y

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

1. #1、加载数据\" K! [' e  W2 |, y, V! Z0 P
   
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )1 A- V: g$ n\" w  X0 N9 S
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  \" ]6 U6 U# p% Q* ?+ m. a
   
4. 
   \" u! T1 A8 I5 Z
5. #2、建模
   # F& G; \) Y2 v) X. i# F; Q, d
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)/ J/ X; u9 h\" `. h8 U
   
7. % b  Z5 g  X( o& o( Y5 k2 W1 I. u
   
8. #3、模型评估+ r& |( X$ P3 q9 q+ T% x$ n& F. X
   
9. summary(glm.sol)

复制代码

1. ## 4 L9 b- _+ c* I. k
   
2. ## Call:) Z2 N$ T0 X: t* N
   
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)
   ( ?. O5 w9 P+ C/ r/ E3 e0 q, D
4. ##
   $ v5 M) ?5 p& Z
5. ## Deviance Residuals:
   6 i2 }( A' u9 `* u, U
6. ##       1        2        3        4        5        6  
   1 ^2 M( c2 V- U$ P; |& |' T
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
   8 D% R( w\" ~% {9 _! r
8. ## 7 k  e) g. w6 n, f/ R  w: Q
   
9. ## Coefficients:
   $ U  _* z  I+ e& R
10. ##             Estimate Std. Error z value Pr(>|z|)   
    8 O( B4 t% P# t2 c# m. V
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***+ o# ^# H8 j. \& `1 ?
    
12. ## x             1.2459     0.1119   11.13   <2e-16 ***/ D  A7 }+ ~6 e* _+ @
    
13. ## ---
    8 C9 e3 S' U0 }0 X$ Y) T9 r
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\" D- L4 c) E\" ~) P
    
15. ##
    6 @. ?, s9 j# [0 Q' i1 |- k+ x
16. ## (Dispersion parameter for binomial family taken to be 1)3 d2 k9 h5 ^5 {
    
17. ## 7 Y/ a/ |  g2 I1 x7 t: D/ C
    
18. ##     Null deviance: 250.4866  on 5  degrees of freedom\" r  i* r4 q2 L4 @; s7 s/ F7 x
    
19. ## Residual deviance:   9.3526  on 4  degrees of freedom
    / D6 n\" [) X1 s  P/ F8 v. t
20. ## AIC: 34.093+ r! e$ w  t* ]* m: A3 H
    
21. ##
    5 G* W: [; x\" d: t
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率
   2 X6 b! T* O4 `9 Z/ O/ A6 ?5 L
2. 
   8 Y2 N\" h6 A& L- p! {9 Q; l$ M2 n2 E* E
3. #4、预测/ J0 i% [  y* k6 E4 [' |6 o4 s+ y
   
4. pre <- predict(glm.sol, data.frame(x=3.5))
   3 I& U' _& ]# V% t5 w& y& [9 b
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1
   . @5 k  a( ]\" g7 \- K1 j
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1
   ) A! [! T; [5 T# ~
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x # f4 K2 i; L4 j! P: G1 t5 K
   
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：6 W  z6 ?; s& U% P
   
2. d <- seq(0, 5, length=100)2 R1 _. K8 r# [* E; U5 h3 K# F
   
3. pre <- predict(glm.sol, data.frame(x=d))
   9 m( N) g$ z& N/ J\" |
4. p <- exp(pre)/(1+exp(pre))# c7 h4 N$ d2 y( O$ z
   
5. norell$y <- norell$success/norell$n4 e8 }# y& R- u; {
   
6. plot(norell$x, norell$y)
   1 l$ z' _* o* v  H\" c- J
7. lines(d, p)

复制代码

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

复制代码