Logistic回归实例2

2744557306

# logistic回归
实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。
/ F1 f5 r; e  d
$ ^& v, C; c: b0 N; AR软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。

! y) x4 t" ^! D! H1 N* O9 v
2 N, l* c3 n! q+ A1 C有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。
9 h2 x8 _! w3 l. C' p9 N" _, j
3 e& T4 W" r& b: Z+ s( D1）正态分布 正态分布族的使用方法： 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)有完全相同的结果，但效率却低得多。
% F( a1 I, n  C. ~6 S" E
8 A* ]6 y$ W/ X; L. W2）二项分布


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

1. #1、加载数据2 Q9 g5 T$ l5 J
   
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )! \8 d) n3 }4 ?\" {( A) }, i: s
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  + J3 E7 P/ ~/ A\" S7 \0 I
   
4. 
   : v5 Z( O5 m/ I3 ]6 O
5. #2、建模
   + s\" r1 O2 A- S1 w  h, \( u2 D
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)9 U, j: l+ S( z7 d! Z- z
   
7. 
   3 l2 D  H. k0 o$ [
8. #3、模型评估
   5 j/ ^\" N! c0 J, }
9. summary(glm.sol)

复制代码

1. ##
   - J! s+ A\" R0 H1 f& H
2. ## Call:7 G- g7 j4 H/ R: H' Z
   
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)4 [0 V$ ^3 A9 z7 X5 y
   
4. ##
   ; B# b% D  q3 m7 q
5. ## Deviance Residuals:
   $ G/ v7 Y* m' {/ n* G0 o& l4 j* E
6. ##       1        2        3        4        5        6    {& C: R5 k8 F$ [\" D\" @3 J
   
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  ; z! t3 |3 \( b( s. a2 u
   
8. ##
   3 u; ?) n) \1 w' J1 G( w1 J* n
9. ## Coefficients:
   4 M& b\" w# \* J+ x, H
10. ##             Estimate Std. Error z value Pr(>|z|)    5 ]# @/ }$ u\" u) \2 W
    
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***0 r7 ^; Y7 [5 \# K+ M
    
12. ## x             1.2459     0.1119   11.13   <2e-16 ***
    $ _* X* h( G9 W1 E
13. ## ---6 V( q, }3 Q. W3 G, Q! g1 B. |
    
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1: I8 V+ g7 ~+ T2 c, j
    
15. ## 4 d+ C  u; i2 @9 o+ [' X' o% O
    
16. ## (Dispersion parameter for binomial family taken to be 1). s1 w8 T: c0 h' l) e* N) S
    
17. ##
    ' Z: o) q5 q( v% r* L# }* `+ o  F
18. ##     Null deviance: 250.4866  on 5  degrees of freedom
    - m* n0 }$ g! a+ p7 c/ E
19. ## Residual deviance:   9.3526  on 4  degrees of freedom' b' G) {: Y# p5 V0 M\" M& ?) k
    
20. ## AIC: 34.093& Q) ]; N  D: M& a7 q  B
    
21. ##
    % P5 t/ s* S! R0 o5 n+ v
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率9 P, E  C- b* T
   
2. 
   * S' G& {2 e3 u! H* a& b
3. #4、预测  o9 O4 O1 q' z* w
   
4. pre <- predict(glm.sol, data.frame(x=3.5))
   1 g. F2 T, z4 X+ U
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1
   * D: ^; s  q% Y5 w( w% e1 D
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1
   ' c9 W( K5 I. G6 b: v1 L, Y
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x 4 C; Q$ l. H\" k. P
   
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：
   \" i5 z  Y8 K8 f' |2 I: @, J
2. d <- seq(0, 5, length=100)8 p& {$ q# p6 r$ ^
   
3. pre <- predict(glm.sol, data.frame(x=d)): T/ A6 R7 K% Y# X
   
4. p <- exp(pre)/(1+exp(pre))
   ! g\" q* o0 i2 s8 _7 s# k0 W
5. norell$y <- norell$success/norell$n
   4 r, {& f: d; G; X% G) z2 ^
6. plot(norell$x, norell$y)3 ~\" k. H0 R) v, |/ F# f\" c1 ~% G1 }
   
7. lines(d, p)

复制代码

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

复制代码

. l: m6 l9 G5 F2 ~