Logistic回归实例2

2744557306

# logistic回归
实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。
" `1 |! x; u  T
0 I# I/ d5 g+ a) L* ?0 X/ k) bR软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。
* `! c" _- H) d5 @3 L: f3 A5 m. {" h

有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。
& c6 g2 K# _# d! J/ v
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)有完全相同的结果，但效率却低得多。
/ ~6 A9 I9 t5 K
2）二项分布


  _6 U3 G; Y- m8 C  L- C/ Qlogistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。
# Z! ~7 }1 ]5 n; i2 d
  C7 [% \; O! m& b$ ^) h& ]Logistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据, c& C+ i* Z\" R, Q1 _
   
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) ): Q& u, Y' K0 ]& f\" l
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  3 ~/ Q. W# G8 v: D2 q5 Q5 J: `
   
4. 
   % n: V1 T- l/ C
5. #2、建模
   3 s) n4 Z\" R/ r3 W1 N( ]
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)1 |\" e! a0 B8 Y6 M\" N
   
7. & |2 {  g- e$ X7 ^7 U* a& H: v
   
8. #3、模型评估
   1 k' U5 F& O0 J& Y! S
9. summary(glm.sol)

复制代码

1. ## ' k8 T; r; A# U. Y5 p
   
2. ## Call:7 P# V/ j' W/ R9 N% V
   
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)
   ) @' d' w\" ~9 O8 i4 R* y: @3 l
4. ## : d/ C# X6 I. [+ I\" O/ P\" G$ m
   
5. ## Deviance Residuals: $ \( G- }* n- T$ B! A) z! e
   
6. ##       1        2        3        4        5        6  ( a% [( q, v  o, S9 w
   
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
   2 {  Z- j: b. N* A: a
8. ## , h& w$ N( \7 H
   
9. ## Coefficients:
   , d: a& w: [# @, K3 g
10. ##             Estimate Std. Error z value Pr(>|z|)   
    - _) z7 A1 I9 V5 J/ t
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***0 X9 f4 E9 N* Z: |* C\" g
    
12. ## x             1.2459     0.1119   11.13   <2e-16 ***
    . C: M5 b2 Q# G\" C( u
13. ## ---
    . ]% g0 y5 i4 w/ g  n
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 11 G5 n, ~( q3 |, Q% a, v
    
15. ##
    ! O5 _9 \' Y\" S8 o4 G! R
16. ## (Dispersion parameter for binomial family taken to be 1)
    , {$ q\" T$ X/ O
17. ##
    ; A: r$ B  J+ n7 v7 a
18. ##     Null deviance: 250.4866  on 5  degrees of freedom
    \" g/ u$ P( ^0 |% `' o! }
19. ## Residual deviance:   9.3526  on 4  degrees of freedom
    . [, o$ }- b$ B  t: C( S; }* I
20. ## AIC: 34.093
    8 @  O( e: d* l
21. ## / o2 t( @7 h) f& p0 e; J1 W
    
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率, W* n6 X* A: j% a4 ]' |, H
   
2. 
   ; E8 h  _  }$ D1 D
3. #4、预测
   9 j7 L' ?) S2 {, |3 _
4. pre <- predict(glm.sol, data.frame(x=3.5))
   2 T: U7 J! j1 `
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1
   ! j0 O# N! f$ Y7 E3 u4 E% W3 p
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b12 b\" s# F( U( m& ?% b6 F
   
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x * f5 Y  C7 \: e- l% c
   
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：& b- |6 k0 s! Y4 v& ?- l+ g
   
2. d <- seq(0, 5, length=100)  W) d( Y1 J1 l0 `4 z, }. o& p
   
3. pre <- predict(glm.sol, data.frame(x=d))
   4 m* \( N$ K. F& N4 ?
4. p <- exp(pre)/(1+exp(pre)); |  S1 P  D! w  ?: w! z/ `
   
5. norell$y <- norell$success/norell$n2 k1 @\" O6 X! q6 R
   
6. plot(norell$x, norell$y)
   * \) H1 n: |5 n& w( Z/ o8 T: t
7. lines(d, p)

复制代码

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

复制代码