Logistic回归实例2

2744557306

# logistic回归
实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。
; Z$ F* Q& v* k) ]
8 `6 J( C8 n5 [5 I: k& A, Q3 vR软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。
) H. ?) E' `' N9 u
) @1 \3 K8 N6 W% W3 Q
有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。

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)有完全相同的结果，但效率却低得多。
& [4 |1 @% n& P  o3 ~
- ~5 x* A8 V1 U/ y2）二项分布

. H6 [/ @8 ^1 V% C: `
  y4 D( Y: M2 U$ ^+ Z1 @logistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。

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

1. #1、加载数据
   ; d2 ]- J0 g$ G
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )% q) K% H% a! |1 p/ ^: C' R
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  * ~0 v% ?6 f& A6 D. e2 s
   
4. 
   2 A$ @1 U' U2 x% W4 Y$ u, T9 f$ H
5. #2、建模\" U& }% J* r2 u. U* u
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)* p6 P' D7 r/ g  M5 q
   
7. ( D8 R' _2 P\" m# B% @
   
8. #3、模型评估' {% c/ J; k* p; o1 ~4 T* K
   
9. summary(glm.sol)

复制代码

1. ##
   + ^( a+ d: e0 N* R- q
2. ## Call:# U9 ?% u0 e, s0 i1 X
   
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)\" u; m- Z- Z9 f
   
4. ## 1 E0 Z5 ^, ?8 K) G
   
5. ## Deviance Residuals: 0 x% B0 Y; ?& l# Z5 t) C
   
6. ##       1        2        3        4        5        6  
   ! C# V: R/ v% ^5 Z
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
   8 k! l* A( F0 t
8. ## ) N9 `! N/ s1 R; P
   
9. ## Coefficients:+ [2 ^+ U0 z* o' Q' n
   
10. ##             Estimate Std. Error z value Pr(>|z|)   
    - q9 d$ x% c  f
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***, X( \9 g# C1 k- Z
    
12. ## x             1.2459     0.1119   11.13   <2e-16 ***
    ! |7 @7 g- K\" K/ t6 l
13. ## ---
    2 v7 V; ]* W: T
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14 Q4 ?8 _) l7 P8 W. G% [$ r
    
15. ##
    , Y6 G* u; [7 r1 B9 {# F
16. ## (Dispersion parameter for binomial family taken to be 1)
    8 H3 K) h& g1 l7 T' e
17. ## ! n; G) @  T, q( ^; z4 c  `( \
    
18. ##     Null deviance: 250.4866  on 5  degrees of freedom7 ?\" k7 D  `* f
    
19. ## Residual deviance:   9.3526  on 4  degrees of freedom
    , |/ f6 T6 p% A: f2 [( l
20. ## AIC: 34.093
    . E  p  X$ M\" B4 m; F# j6 k
21. ## ) y6 \0 H: z2 Y# I
    
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率0 v% \4 d: O# g& ~1 A
   
2. ! z, r4 r4 ]: W& {0 u
   
3. #4、预测: ]' |  Q% x3 w* t+ ]5 n' f
   
4. pre <- predict(glm.sol, data.frame(x=3.5))
   # U! X  m8 T$ |; G2 E
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1 + A$ L7 @, d  K, A* q3 E
   
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1
   . V  g( j! y0 B0 q7 }  x
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x
   4 H0 _! I2 M3 ^! V
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：: ?* a* t) e  {( i( D7 T; s2 H
   
2. d <- seq(0, 5, length=100)
   $ C* e2 l$ e, G& }4 _6 l9 |
3. pre <- predict(glm.sol, data.frame(x=d))
   / n% G; A* R0 ?, Z# t- P\" |
4. p <- exp(pre)/(1+exp(pre))! `4 T% G& I& X/ C
   
5. norell$y <- norell$success/norell$n
   . E% s0 X! f- u) E
6. plot(norell$x, norell$y)
   0 v% o* b# t9 r1 d' p9 W0 v9 f* L
7. lines(d, p)

复制代码

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

复制代码