Logistic回归实例2

2744557306

# logistic回归
; {8 k/ o$ M9 P: W3 L/ q实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。
5 T! C& ]5 ]5 o+ q
; C* e# k% |: w, \* I1 y3 w# M9 f* UR软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。

% [6 z0 Z3 `! A, Q$ x3 i8 M
/ V  R' S! @# _4 J有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。
* l; h- a( j% D5 \
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 u9 M% {% n; U+ n
; S# }6 f, a2 u$ V5 b! \& m2）二项分布

4 k0 T- A5 [! `( F$ a
! h/ p2 _& x8 v4 @$ Clogistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。

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

1. #1、加载数据
   & M+ O0 i' O0 \4 W' a
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )
   ' R/ s- b, }4 O. M/ Y6 M( T: T
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  
   + ?$ t9 A1 x5 X4 M- h  X. t7 g
4. 
   5 M  _7 K, u1 y: t9 q/ }$ P# L
5. #2、建模\" f, B\" @, A0 Z
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)7 m1 Y) ~/ W# M2 f+ Q: a
   
7. ; b/ ?9 J% W2 ~6 N$ e
   
8. #3、模型评估
   ; Y! _0 @2 `5 |. U6 x# h
9. summary(glm.sol)

复制代码

1. ##
   8 V) v  X2 K/ ?2 M% X+ |2 N\" b
2. ## Call:
   % a/ v0 M# I: P5 U9 b  u5 V/ C1 i
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)
   ; e' G5 e* H3 I. _, I5 r, D
4. ## 4 i: V3 a4 {) B4 W* z' `
   
5. ## Deviance Residuals: 1 G& K  y! @  R9 y# j% c
   
6. ##       1        2        3        4        5        6  + O- G* ]: g9 |
   
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
   8 S( M. T: \% T% o\" M0 g& X! {5 E
8. ## ( O) n6 o& a% _& q: s7 e
   
9. ## Coefficients:# U# s4 I. r( _( j\" z) P
   
10. ##             Estimate Std. Error z value Pr(>|z|)    $ k3 c2 ^1 ?, u5 H
    
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***
    4 S/ {2 N; w. F  c+ P' f
12. ## x             1.2459     0.1119   11.13   <2e-16 ***
    4 T2 T  K$ q& b3 @
13. ## ---4 e6 @& X4 `' f7 M$ d
    
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ) {$ h8 F8 D\" W& l0 Y# ~. [
15. ## 2 U) Z4 a: F- l; C4 ^& n4 N
    
16. ## (Dispersion parameter for binomial family taken to be 1)
    0 i2 ]) P4 z5 i8 _( G# ]+ `2 ?
17. ## # Q1 I5 E! U+ F) A: c
    
18. ##     Null deviance: 250.4866  on 5  degrees of freedom
    % M8 W! i+ h5 B, d9 T, E
19. ## Residual deviance:   9.3526  on 4  degrees of freedom
    ' a* a/ _& [1 @1 Q; {/ y( d; K\" A
20. ## AIC: 34.0939 g2 ~/ n+ J! N* @2 j  A& z
    
21. ## * j+ v9 u& B( p  v8 ^9 V
    
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率$ B) J* x\" P5 u0 m' ~
   
2. 1 G1 w. w+ Y8 e) D  ]  S% ~
   
3. #4、预测
   : |, H' T2 p* ?- Y2 [( M5 |6 {& w4 \
4. pre <- predict(glm.sol, data.frame(x=3.5)); G5 ]% x5 F9 i' `) o7 N
   
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1
   0 F3 H* e7 B' u2 X8 `$ \. i
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1
   - m$ ]- x! J( `8 w$ Z) P# _
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x
   7 L( O$ |6 _: l! F' M
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：0 F. C. w) k' x
   
2. d <- seq(0, 5, length=100)
   * I; G5 B6 t  ]. C' I; z
3. pre <- predict(glm.sol, data.frame(x=d))
   - N5 K2 k+ u8 h\" G/ B  {. T5 j
4. p <- exp(pre)/(1+exp(pre))
   5 L8 b3 y\" ~0 Z' M5 m' G4 D
5. norell$y <- norell$success/norell$n. ~- M: J  _% P/ n* W8 l
   
6. plot(norell$x, norell$y)
   4 t$ p5 Z- Z) u9 V* p* g
7. lines(d, p)

复制代码

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

复制代码