Logistic回归实例2

2744557306

# logistic回归
6 f: v5 ~2 @6 i4 k实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。

: d2 w7 V( |2 YR软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。
2 m" ~0 ~5 E! T; M+ A1 X
0 m8 P- I, |/ C
有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。

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)有完全相同的结果，但效率却低得多。
1 G) S2 `; p2 h1 m
2 m( j! i& o" R2 v2）二项分布


: c5 W& @& K! r6 [) e# _- Y$ nlogistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。
# f) _6 w' g: j& L! o
Logistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
! h& @, i$ M0 C2 s实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据
   7 q% J+ O  s! C# H/ p$ ]4 Y4 }
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )
   ( a, j  y% O- o, K& _6 F7 ~+ H
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  ! Z9 ~* e5 Q  n* M3 ^: G/ D
   
4. 
     T7 ?, _. {* e5 R* g% B
5. #2、建模) q% `$ i5 g* R; Z& I\" S7 ?0 _
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)1 N1 P  y) F0 D) Z* n! S
   
7. 
   3 I4 }1 m8 f6 K* u* }8 e0 c
8. #3、模型评估0 D4 b, N& y( I8 t
   
9. summary(glm.sol)

复制代码

1. ## - p7 u+ J\" W1 m
   
2. ## Call:
   0 ^# D( O6 H; c- l
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)
   & q; L/ }6 k# v7 G0 a
4. ## 3 _  f& H7 _7 Q6 F/ g
   
5. ## Deviance Residuals: ' a# O. s: H& ^4 k
   
6. ##       1        2        3        4        5        6  1 ^5 Y\" r$ k! t, |9 J( @
   
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
   ) ~9 S5 @) @$ w  {& e9 A: j4 g
8. ##
   . Q# E) @! c3 F( l# E* G; V3 x4 V\" Y
9. ## Coefficients:
   \" F\" H+ R$ j- B; R; e5 x
10. ##             Estimate Std. Error z value Pr(>|z|)   
    ) C7 T/ L+ L5 q
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***
    ! @2 G8 {! ~  y3 w( g( U
12. ## x             1.2459     0.1119   11.13   <2e-16 ***
    ; W+ J( r, u$ ], M+ p; p5 g8 o
13. ## ---
    , `& M. h+ F0 t/ @) I& i
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 19 W8 V: c4 Y1 c. F% `7 O' @+ a8 S
    
15. ## 7 D3 t3 o& ]9 J) k: ?: }' U6 c
    
16. ## (Dispersion parameter for binomial family taken to be 1)* z4 \/ U( T8 U! P3 Z- _; z
    
17. ## ! K8 U8 s\" h7 C1 E
    
18. ##     Null deviance: 250.4866  on 5  degrees of freedom
    0 \8 S( _. v& V$ v) c3 i( C
19. ## Residual deviance:   9.3526  on 4  degrees of freedom
    - v2 |- l1 f! u; |6 c
20. ## AIC: 34.093
    + `( _; U, J) K+ O. C
21. ##
    + @3 h  ~( ?! y\" f, r7 L
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率
   ' m) ]& p% f* q1 f8 Y
2. % i8 F1 g& A5 l9 q4 c8 Q/ `2 [
   
3. #4、预测
   ( ^5 |: b  @0 o1 ?; g% z
4. pre <- predict(glm.sol, data.frame(x=3.5))- m1 O; \7 H2 U  \- o# S
   
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1
   & Y% ~, S* B+ e; r
2. ## 0.742642

复制代码

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

复制代码

1. ## (Intercept)           x
   \" q3 k( r. B* ~
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：
   9 S1 z, Q% d4 |! z# r0 J
2. d <- seq(0, 5, length=100)
   \" z2 Y, Q  K. a* N- Z, k
3. pre <- predict(glm.sol, data.frame(x=d))
   8 M) m8 t% a/ W) ~) E  l\" B: P# B
4. p <- exp(pre)/(1+exp(pre))+ f( o! i) k3 t
   
5. norell$y <- norell$success/norell$n4 w# M. n9 h6 R& i; O
   
6. plot(norell$x, norell$y)7 Z7 Y0 Z$ R( l- a
   
7. lines(d, p)

复制代码

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

复制代码