Logistic回归实例2

2744557306

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

0 Q2 W2 z: i& Y* o  F+ xR软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。
8 a4 ^* x. H5 e4 [9 I) y

有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。

" j7 I) O" O* k4 P# m4 }( [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）二项分布


logistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。
* ~4 r$ k) M: M  }* D( s
2 H! V# t1 m3 n8 W' \Logistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据
   + e) Q\" h1 F8 W
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )
   \" M8 K- n+ C4 Z\" p% @& y
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  , p\" d2 s( c. D. P; |
   
4. 
   / I' [- G# E9 m% m1 K
5. #2、建模* U8 N; Z: K/ z  ]4 _6 L
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)- o! T# J0 Q3 d
   
7. 9 O, f\" X# y* w; x+ z
   
8. #3、模型评估
   * o' ^. i. V$ W! F
9. summary(glm.sol)

复制代码

1. ##
   3 V7 P% T7 |. w$ Y
2. ## Call:
   * i9 t5 r- B) W; |% z4 w
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)* y% p\" C9 ~7 ~$ f+ Q
   
4. ##
   # o2 a2 t3 x3 J: C9 @4 T
5. ## Deviance Residuals:
   1 W7 Q/ K+ Y$ ~9 h/ u
6. ##       1        2        3        4        5        6  \" X$ A6 C6 X9 M$ X
   
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
   7 V- A/ H& n) T4 v/ W
8. ##
   & f' G\" v+ t# d6 i+ ~% `
9. ## Coefficients:* d, X, b( V. Q) @3 O4 p, f* G
   
10. ##             Estimate Std. Error z value Pr(>|z|)   
    ) o& ?/ D* D* Q) W, c5 k\" Q6 R
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***
    + Z6 `. j% Q6 y\" U; ~* J
12. ## x             1.2459     0.1119   11.13   <2e-16 ***/ ^# z, M' c$ m, Z6 k. W
    
13. ## ---
    % ^$ B; F2 J( J9 S* z) _0 T
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  G: R2 E$ `3 k
    
15. ## # w8 z8 t% P4 ~7 f- L
    
16. ## (Dispersion parameter for binomial family taken to be 1)
    + D  l# E0 f9 B5 B% s# K) G1 p) K
17. ## # b! ^( H& j& w- {\" z+ p\" X; W
    
18. ##     Null deviance: 250.4866  on 5  degrees of freedom+ n0 d6 t5 R5 r7 {
    
19. ## Residual deviance:   9.3526  on 4  degrees of freedom! @1 [! x& c; @3 f* Q9 Q\" {. h\" U
    
20. ## AIC: 34.093
    - j' d. z9 u9 v  g* H7 I
21. ##
    . b' H$ o. v, r( X( m' ^\" r2 X8 [
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率; u3 A/ }0 p) I! P; g
   
2. 
   - n2 ?0 B. G8 t/ a
3. #4、预测  V) ]9 V5 Q& E: \5 j+ W6 r
   
4. pre <- predict(glm.sol, data.frame(x=3.5))' r. G5 v6 Q( F9 \& T% D: S7 B
   
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1 + S2 u. M; M3 E9 a! G
   
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1
   4 @0 R6 I' ]; H# _1 _
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x 0 C7 G' q% R% J4 f% Q
   
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：5 l2 c! S( f4 `4 \/ `
   
2. d <- seq(0, 5, length=100), P9 w; c7 s2 ?  z: {
   
3. pre <- predict(glm.sol, data.frame(x=d))4 \0 m! G& U1 ^: n  z! y. I
   
4. p <- exp(pre)/(1+exp(pre))$ {1 z) |: B! x3 p
   
5. norell$y <- norell$success/norell$n
   1 _0 _6 K/ n+ y9 g6 F/ V\" Y
6. plot(norell$x, norell$y)+ x: s! s/ v/ l& E1 m& B- u
   
7. lines(d, p)

复制代码

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

复制代码