Logistic回归实例2

2744557306

# logistic回归
实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。
  `/ q- o3 ~. |) a2 O! O0 q2 K3 y) g
R软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。

, P: N2 n3 m$ j
有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。
5 L( U3 X: X. S& I
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)有完全相同的结果，但效率却低得多。

/ ?! t. z& ?! x9 F; |2）二项分布
: z' k  z7 W) l
4 Y6 s3 K) ]: I) `& f9 _4 M- f* u9 r
logistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。
/ _4 A* j/ J+ j
6 K9 B* _0 J- W0 X+ y6 y! S2 _Logistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
9 l& n" u# }- P" j实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据
   7 V3 V9 J\" `1 u, b/ E0 @! N3 \
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )# R4 s  G- p6 x6 ?, X\" G: W( D2 S
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  
   - }* M2 }% f! r\" b  S' o
4. / D4 j8 H; g/ H- q7 _, R4 o  ^
   
5. #2、建模2 d1 E6 E8 S. ^/ @  ]4 ~$ S
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)
   % @, m) Q+ c3 e
7. 5 }0 z1 p$ p9 \$ V+ R2 K
   
8. #3、模型评估7 f* I& g$ X' I  H
   
9. summary(glm.sol)

复制代码

1. ## $ Q2 q6 [4 e& p
   
2. ## Call:
   ' X! I: q: v( y
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)1 B' l' @\" w  H5 _- d) w& f
   
4. ##
   & @) y  f- L6 m) Y# G8 q. d8 y
5. ## Deviance Residuals:   B- ?. W) O4 J0 r: c
   
6. ##       1        2        3        4        5        6  
   4 G\" P! U6 Q. W+ y
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
   6 J8 @* [  U9 a2 ^
8. ## . f4 X) w- m2 g% P
   
9. ## Coefficients:
   2 ]; M* M! N, M7 o, Z& k
10. ##             Estimate Std. Error z value Pr(>|z|)    / t3 R6 T( H9 n6 T4 d
    
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***
    % J* |( A& K, S+ u
12. ## x             1.2459     0.1119   11.13   <2e-16 ***
    / f' [5 N2 r2 s  u+ i1 A
13. ## ---3 n' v1 h& t/ p9 z9 h1 Y) L
    
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ! U# T3 z) K4 O2 w5 Y# a9 ]
15. ##
    ' R, t2 n  J0 ^
16. ## (Dispersion parameter for binomial family taken to be 1)& p- m/ l) A0 M4 ?1 r7 o
    
17. ## $ w- X3 d1 X2 ^6 E8 [; [$ i
    
18. ##     Null deviance: 250.4866  on 5  degrees of freedom
    3 H8 ^! h; }2 }5 @
19. ## Residual deviance:   9.3526  on 4  degrees of freedom- L, b* ^1 K3 s8 w# G
    
20. ## AIC: 34.093$ Y8 _. l  U2 ~  L9 g\" n
    
21. ##
    # I6 R; V5 v) h% C, R
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率
   ; S; X0 [5 _: Z- h4 T
2. , V\" V, ~4 ]\" C1 @* a% B7 m
   
3. #4、预测) i3 r9 M, O3 |* F8 \6 o' }
   
4. pre <- predict(glm.sol, data.frame(x=3.5))$ K% s* e7 c4 t\" h4 @
   
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1 8 ?: d* L; k& ^
   
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1; Y# {# o8 r. [' ^
   
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x 3 C# W# Z, C, k! |( x
   
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：
   + H: z4 N1 s5 q/ H
2. d <- seq(0, 5, length=100)' z3 T( O9 ^2 @6 ^\" f' e! N
   
3. pre <- predict(glm.sol, data.frame(x=d))4 U/ d; Q0 p  _5 r4 S
   
4. p <- exp(pre)/(1+exp(pre))
   8 Q1 e) V* \7 x2 F$ L5 ]1 E
5. norell$y <- norell$success/norell$n
   + B1 ]0 J, r\" d7 P! P
6. plot(norell$x, norell$y)  M* B3 J. M2 V+ O: T
   
7. lines(d, p)

复制代码

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

复制代码

5 ?8 D* x4 n+ y; S5 ~