Logistic回归实例2

2744557306

# logistic回归
实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。
, i  ]7 l3 n+ H0 [' I. k
R软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。
8 _! |  k' j  ~' l+ W9 Y

有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。
5 r- n- v0 H% n7 J- }; v$ S
9 ]) p* G8 ?/ l! D( x! \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）二项分布


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

' X- f3 w: Z: JLogistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据8 U7 }2 V) h% J
   
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )) g7 u1 b; ~, ~) p/ d8 @; V1 x% P
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  ( S7 R0 D9 w9 A9 y7 T+ n' \6 D. |
   
4. 
   6 @/ f/ r$ V& C0 \9 }3 q
5. #2、建模7 y# U/ D& e- o, g! L
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)
   2 Y& n0 G7 L# }- \  N
7. 4 u/ j: q4 Q: b/ R
   
8. #3、模型评估5 }! i/ a  N: \0 {  S
   
9. summary(glm.sol)

复制代码

1. ## 4 t) c3 |5 m% ?
   
2. ## Call:
   + v4 a0 z& h( K! N2 B, t
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)
   7 I% A! x6 Z1 J6 s0 a6 L9 D  O
4. ##
   4 f) b! w/ ~% h3 G# k+ i
5. ## Deviance Residuals: 9 ^5 ^4 d2 T' z9 q: K
   
6. ##       1        2        3        4        5        6  2 ^% h$ n6 H( ^% {% q
   
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
   $ q0 n; i% \1 p# Z/ `
8. ## ) p1 s0 p0 P3 e% c
   
9. ## Coefficients:
   * K, @  H+ w% ^! a, L3 I
10. ##             Estimate Std. Error z value Pr(>|z|)    7 {+ F! s+ B1 s$ j1 X
    
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***
    * q4 R  x$ P+ O, ]0 k4 b) ~
12. ## x             1.2459     0.1119   11.13   <2e-16 ***
    & R* u! V9 ~! O9 c2 C+ t
13. ## ---
    . [% D; X7 O# |/ F: N
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- q5 P& r2 ~: b# S\" v6 w\" P
    
15. ## \" }1 Q( b9 y\" G
    
16. ## (Dispersion parameter for binomial family taken to be 1)+ B\" T6 G2 f+ M) i2 Y( m+ E1 K
    
17. ## 2 n. [8 r) p. g1 G( X' k
    
18. ##     Null deviance: 250.4866  on 5  degrees of freedom- |4 i3 e3 a. D0 M+ {) o% N# P
    
19. ## Residual deviance:   9.3526  on 4  degrees of freedom; F& x$ W9 `. j* I% [7 t6 D5 Q
    
20. ## AIC: 34.0930 t, l5 F\" V1 ~6 N# ]
    
21. ##
    1 W' y\" q! J% N6 W
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率' `/ l- T# }& s! K
   
2. % I7 N( i( U4 K4 s4 T3 j8 U
   
3. #4、预测7 L9 `7 D\" `4 [9 [9 l4 w3 M
   
4. pre <- predict(glm.sol, data.frame(x=3.5))( ?  q) B, u+ G3 S7 e# N
   
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1
   # W& C$ S# B# g; l; A
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b18 c/ g2 E+ H; d7 M7 F1 B0 o+ h6 m
   
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x , `: X: Y3 `8 }* G
   
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：
   - N' H7 g7 l8 R0 g% \% w
2. d <- seq(0, 5, length=100)
   2 ^& q) e  i/ l( p
3. pre <- predict(glm.sol, data.frame(x=d))- P& z5 u5 n8 U  G
   
4. p <- exp(pre)/(1+exp(pre))\" u* Q8 e/ G' H7 R( j: w
   
5. norell$y <- norell$success/norell$n, o% c$ b- ?% p
   
6. plot(norell$x, norell$y)
   ) d5 P: ]- z. A: N+ ~, U
7. lines(d, p)

复制代码

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

复制代码