Logistic回归实例2

2744557306

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

R软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。
, l+ M4 Q* v5 A4 H' Q
7 e' Q% H# _' ~" W
0 f7 h/ V8 S4 c有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。
' j" R" a( @" G, q
: ^$ D" R: x- M' z1）正态分布 正态分布族的使用方法： 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）二项分布

& h5 E/ S2 F9 j: i: B: x. i
9 B# m" b9 _+ ^3 _; e3 {6 Ulogistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。

: N+ z0 M' T. u1 [7 oLogistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
9 B5 `, _8 s. A' A- Z$ ~5 Y4 [实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据
   6 h+ h\" }  Z* u' X9 [; y! c
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )# G$ Q3 S: u# k6 g7 B7 G
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  
   ! f$ `  o1 x3 n( n
4. 
   , m6 p; `1 t* ~9 s; ]. R\" p
5. #2、建模, ^1 A0 k6 g% H* p! I9 J
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)
   , U+ M, F: [' T
7. 1 K; z% [5 a6 ~' {+ `% n
   
8. #3、模型评估) ]/ h9 d# R\" Q  C* w+ V
   
9. summary(glm.sol)

复制代码

1. ## 7 P- D8 {( x2 C; I
   
2. ## Call:# m) i6 e8 Y  n* X' O& F
   
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)6 w; \& Y6 P7 P, b: f3 j6 Y6 ]
   
4. ## 7 E+ a. {2 b! z  m\" `( n$ h7 v
   
5. ## Deviance Residuals:
   , R* p0 y5 K$ e0 m7 ?2 d& H4 u' I& p$ l
6. ##       1        2        3        4        5        6  
   : d\" F! w$ Q5 {3 c! k
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  7 s$ k6 c; D0 b4 Z) l2 d
   
8. ## 1 J7 g- S& ~5 V2 a7 j
   
9. ## Coefficients:5 X8 R6 e\" x8 _$ R0 {0 l6 ^
   
10. ##             Estimate Std. Error z value Pr(>|z|)    * ?$ q$ \' e. C0 b
    
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***0 j$ B% i) w5 D2 \7 L
    
12. ## x             1.2459     0.1119   11.13   <2e-16 ***. y. k+ f: K; {# y% }6 T7 m' F
    
13. ## ---. b6 ?3 p$ d1 {
    
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    2 J, j! }' {3 u- `9 |4 C! |- f
15. ##
    ! U: g$ a2 q  k* L
16. ## (Dispersion parameter for binomial family taken to be 1)
    \" `. t2 F0 M. a8 {; t
17. ##
    / f* d% f2 Z( i  L7 p( |- P
18. ##     Null deviance: 250.4866  on 5  degrees of freedom  b9 @5 b1 g  _
    
19. ## Residual deviance:   9.3526  on 4  degrees of freedom$ o* F$ r9 l8 s2 \: z& I
    
20. ## AIC: 34.0936 I( N\" v; c+ j8 k( P9 F* y% U
    
21. ##
    ! n$ `. _. }2 d4 M2 c
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率( I' T( P+ u; Y, D/ p& o# W
   
2. 
   % J/ |- a2 G# z$ ]' |
3. #4、预测
   + ?  Y4 T* u; e) D6 Y
4. pre <- predict(glm.sol, data.frame(x=3.5))8 _) H; X4 Q+ t4 j
   
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1 ' H, [\" \% W! I: O& b
   
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1
   9 k% _& c  K' v: g5 F4 t
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x
   $ K5 g0 G$ E' ~4 c\" Q
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：
   - H' S9 `; ~* J& |8 G6 b7 f
2. d <- seq(0, 5, length=100)
   9 s! |4 g' [4 f* v8 u
3. pre <- predict(glm.sol, data.frame(x=d))
   $ X6 S2 R: y/ a0 ^) j
4. p <- exp(pre)/(1+exp(pre))' ?! \% d5 |9 I. x
   
5. norell$y <- norell$success/norell$n
   * T4 M2 }\" _1 ~. k! }
6. plot(norell$x, norell$y)
   4 A7 n  b3 p# ]# L1 w
7. lines(d, p)

复制代码

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

复制代码