Logistic回归实例2

2744557306

# logistic回归
( a5 E/ G% j7 K4 h- g  _4 L( E实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。
9 E' i  j* k4 o$ v1 y4 Z. n  e
* ~" S/ A: K& r5 l+ S" ~0 Z  ]. r9 U. gR软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。
: k! C* L( u: N+ _
6 Z! t& S# f2 N1 c; t% o" b, F
有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。

4 S7 C, V% \, {# C. f% M1）正态分布 正态分布族的使用方法： 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)有完全相同的结果，但效率却低得多。
  P( o/ H3 D1 M0 C
& G! e5 \5 ^0 ~2）二项分布


* m  j4 R/ K4 ^+ M; Ulogistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。

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

1. #1、加载数据
   \" t* ~. H: X& n# C  m# k1 {7 b- l( u
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )
   \" \( t# s! x0 I5 e$ P/ F/ P) x' {
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  
   ) l2 W- V$ }) Q( A( }; Q9 A5 A4 ]
4. 
   # d, e4 j5 L\" t5 @! |6 [
5. #2、建模
   0 z( u) I3 {: [  }
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)) Y5 ?, L3 _. H% K
   
7. 
   $ m2 R5 W- H/ o, o- E0 p7 |
8. #3、模型评估3 r% L* Z\" O4 X8 k: w0 j( R
   
9. summary(glm.sol)

复制代码

1. ## . o  E8 \4 x2 n! m9 V
   
2. ## Call:
   9 Q/ m. s3 X# t& w7 w5 T3 |8 z
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)( V( O+ |& M- b
   
4. ## ; X, J; y1 c  ?3 _. V# G; f
   
5. ## Deviance Residuals: 5 J: }+ @& B- E4 _% u
   
6. ##       1        2        3        4        5        6  
     T8 ]* J6 d6 ~8 ]) J' O\" O: ?, `0 {
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  - d\" m1 T* R6 a' u+ i
   
8. ## & }* e5 |2 A. l8 h! ^7 c
   
9. ## Coefficients:# A* a6 J8 A: Y6 Y
   
10. ##             Estimate Std. Error z value Pr(>|z|)    2 W\" f9 m8 `9 ~
    
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***
    ! M4 B7 b2 C4 ?3 @. _- l8 g
12. ## x             1.2459     0.1119   11.13   <2e-16 ***
    4 {4 z- |/ P+ L, a  S3 x3 P. ^
13. ## ---
    \" G/ N+ z! F+ p$ S6 P1 R$ Y
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    9 T  O4 l! i& p\" ?+ Q\" {
15. ## - }3 h6 {# ?9 v: ~& X' j0 j0 _
    
16. ## (Dispersion parameter for binomial family taken to be 1)
    8 z/ k+ O! b* }2 ], u
17. ##
    ) e5 I- [8 V/ H; ~/ E1 J. ]
18. ##     Null deviance: 250.4866  on 5  degrees of freedom, ~1 a5 _. {$ V
    
19. ## Residual deviance:   9.3526  on 4  degrees of freedom* _1 a# g. g7 N1 u3 ~& d/ k( ]% P
    
20. ## AIC: 34.093
    7 j  B, c1 _\" V. m) r: @
21. ## - f3 G; f3 m/ ^( F
    
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率( r' d  \8 f) g4 F8 X4 d$ {, S* e% W' E
   
2. . F& x  J9 _8 n- A3 c7 S* l
   
3. #4、预测
   - j8 `0 ~- T8 _
4. pre <- predict(glm.sol, data.frame(x=3.5))2 }* J' p- i: R3 i* V, L
   
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1 5 a; e4 S7 w3 X! M: t8 E( v
   
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1
   # W, `6 M0 ?. k& U1 X
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x $ F9 d+ A$ I& _: z1 o0 a& T+ H/ `
   
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：1 w8 i( ?7 @1 x
   
2. d <- seq(0, 5, length=100)
   \" K; f- ~. t7 q2 k/ Y
3. pre <- predict(glm.sol, data.frame(x=d))9 n. T& l  j: i1 d! T
   
4. p <- exp(pre)/(1+exp(pre))
     D7 t# D\" a2 r2 D\" v
5. norell$y <- norell$success/norell$n
   - i% ^* f7 g2 a# q! g! D' a# m
6. plot(norell$x, norell$y)
   1 x6 q. |3 O- W' g2 @2 e
7. lines(d, p)

复制代码

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

复制代码

* s/ b, ]$ n4 F' J, C  I