Logistic回归实例2

2744557306

# logistic回归
; U, ]$ L$ W" d/ S实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。

R软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。

/ b  K0 Z( u$ k
有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。

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）二项分布

$ b0 d) I# Z* @
, Q$ g- L. l# V8 q$ d! tlogistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。

# k% a) c1 U* W5 n3 {/ ~Logistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
! l/ X3 W5 @# p: g  R8 Q7 s7 e* q实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据& D5 a/ V* b, ^/ E
   
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )0 t* v: F# y  n6 z& V+ U& V% J
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  
   8 }2 e3 P0 d9 @  @: F
4. . A* z: t3 o: H5 I( l\" {
   
5. #2、建模+ Y& b  \- \3 p) x# ~
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)
   + M* p3 t; M% o# E  z3 D' k5 W4 N
7. 6 @( j! h! Y0 F4 v
   
8. #3、模型评估8 S% W  d, c8 R' [5 V3 n8 N
   
9. summary(glm.sol)

复制代码

1. ## 8 r3 t. x+ F4 V\" _% ~9 x
   
2. ## Call:' W& ^/ e* G/ u4 U! {
   
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)* l' K- G- N* L5 F
   
4. ##
   - [9 ^+ i2 \7 Z. ^\" g
5. ## Deviance Residuals:
   * L5 u* ]! d( A  {: g$ t; r
6. ##       1        2        3        4        5        6  
   7 \7 k$ x. y+ ]3 R9 d% E
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  6 X+ l& R4 [+ P8 U/ g; q
   
8. ##
   : b/ s  _$ `$ R4 ]& K
9. ## Coefficients:
   ' L* [- k  w% Q: K8 Z0 T5 t$ _, ~% V: Y
10. ##             Estimate Std. Error z value Pr(>|z|)    2 ^5 ~; a3 j+ Z
    
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***
    % o' D* X. W/ a, _
12. ## x             1.2459     0.1119   11.13   <2e-16 ***
    # @\" x1 j; `1 L4 c4 m
13. ## ---( a- |9 k\" F  F\" _5 ~' J, Y
    
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    7 C% U+ ]0 Z7 J\" Q0 x0 h/ P0 A
15. ## 3 q& z! N1 {5 Q# ^% U- s+ e! z
    
16. ## (Dispersion parameter for binomial family taken to be 1)5 I$ F9 y  f# Z  W8 g9 B* Q
    
17. ##
    0 T. U! {/ \, i4 F' [4 N. o
18. ##     Null deviance: 250.4866  on 5  degrees of freedom0 k, x+ r* _  X* f. ~, c
    
19. ## Residual deviance:   9.3526  on 4  degrees of freedom4 U1 n* Y\" m. o6 F6 b' S
    
20. ## AIC: 34.0936 ^& `; W0 h! v! J# w
    
21. ##
    4 I, i1 [( f( c' \- |* n
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率
   9 m7 x, t( a8 W- s3 o8 b
2. 
   ' X5 P, d- d; X
3. #4、预测
   \" E4 e+ x7 M, D- k6 h/ q7 O
4. pre <- predict(glm.sol, data.frame(x=3.5))
   , M  u+ ]) w7 L6 a% o
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1
   3 G: p7 p6 C8 x2 M. Q# y+ p0 W
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1' ~0 e- d( C) z+ f
   
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x
   ) M5 ]3 k' P7 i% z
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：
   . o* q( Y3 W; n8 Z+ z9 f
2. d <- seq(0, 5, length=100)' \1 L$ ^( W) f7 b
   
3. pre <- predict(glm.sol, data.frame(x=d))
   , c( b  b5 q' m7 U( X. K
4. p <- exp(pre)/(1+exp(pre))
   & I8 ~! L5 S  q\" e5 G$ ^0 {/ @/ c: `
5. norell$y <- norell$success/norell$n
   . i8 R: g7 }3 K0 O- ~
6. plot(norell$x, norell$y)5 x; ^\" z* u9 A+ ]1 d
   
7. lines(d, p)

复制代码

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

复制代码