Logistic回归实例2

2744557306

# logistic回归
( J0 i" g; u2 o" L/ q0 [0 ^9 w实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。
4 H) T: h& H; H
R软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。
8 E4 G! I) L' ~2 {) a# V, @( T% i

9 @5 _, e+ P# N- V, L7 O( x, p# g! o. h有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。

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

, z" \* w5 }0 x7 }, n# e
! l! D1 I7 p! N# ylogistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。
6 l* t7 n7 k6 J6 q8 C$ f0 Y
Logistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
0 L& S' B! ]: p. T5 e& j实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据
   ' U% f, K7 g- R3 w
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )- }* E; `7 E) a/ n) g4 S2 {
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  9 c. d0 q: \- T1 Y, w6 C
   
4. + }' N: o7 m. |/ t\" C2 \
   
5. #2、建模
   # g$ Z+ o* r/ I  v/ Z
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)
   ! S* D7 C3 a* B5 c
7. 
   1 Q$ q! M( H. d. \8 E
8. #3、模型评估' t2 z! \+ G: Z  @, D
   
9. summary(glm.sol)

复制代码

1. ##
   / _. y$ m8 h. g/ U0 u$ e. A
2. ## Call:  A: i# d; O. O( x
   
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)+ @4 ^# T  F- E& e
   
4. ##
   * P# p! O/ E6 V2 H, P1 M0 R
5. ## Deviance Residuals: ' |3 Q4 ~9 q1 v* w+ V
   
6. ##       1        2        3        4        5        6  & n( e% ]/ V6 J\" s, |! Z- d
   
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
   ; Z6 F* f5 x( K$ Q+ b2 Q9 L
8. ##
   3 ]' e: p. S, b5 y- l
9. ## Coefficients:  D% q7 T5 L- t3 i. n( h, i
   
10. ##             Estimate Std. Error z value Pr(>|z|)   
    3 n6 `! `\" Y0 f, w\" G% f& z* g
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***
    7 H  ~& Y1 w' _: a, X\" C4 @
12. ## x             1.2459     0.1119   11.13   <2e-16 ***
    5 y8 Z7 y; G* e* j: I. C( p% Q, v( ?
13. ## ---9 B; l& s. h' t! l' e% Y
    
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 11 z7 a, g( H& q. n
    
15. ## ) Y\" `4 p% ~( r0 Z
    
16. ## (Dispersion parameter for binomial family taken to be 1); L3 d- O2 g$ n7 T
    
17. ##
    \" M  `8 H: k7 {9 w, I2 W2 u3 }  p
18. ##     Null deviance: 250.4866  on 5  degrees of freedom' C4 a1 i$ X2 e) N
    
19. ## Residual deviance:   9.3526  on 4  degrees of freedom3 g. u2 q( d% `, }# I( u  i+ G3 J
    
20. ## AIC: 34.093
    ( |' b# G) G8 d# J
21. ## 8 E8 U0 H. A: _\" N
    
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率
   3 a; C( M+ w: N' A5 U
2. / O2 n/ K# j\" Z* B, G
   
3. #4、预测, \& @8 ^+ C2 n6 n4 B
   
4. pre <- predict(glm.sol, data.frame(x=3.5))
   5 M) U* k3 B! _
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1
   $ O6 `; U( z* {% a5 s3 L' `8 W6 U
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1% ]# `) J5 e+ N1 R: m$ Y
   
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x
   2 K* Y. v  a4 t5 X1 m5 ^
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：: u; S; K! W8 ^3 j  v6 U3 P
   
2. d <- seq(0, 5, length=100)
   & o8 u\" T, Q  Q# O\" @! |! Z, {
3. pre <- predict(glm.sol, data.frame(x=d))) C5 s* W# U* {/ t1 j\" M% U
   
4. p <- exp(pre)/(1+exp(pre))\" B# |\" I2 C) }/ X# [% Q\" }
   
5. norell$y <- norell$success/norell$n
   # K7 G7 P6 C! }, ~; I3 X* L
6. plot(norell$x, norell$y)% a- P# p7 O/ d9 ]  O! T  H
   
7. lines(d, p)

复制代码

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

复制代码

) y. U  `0 W, B; ]