Logistic回归实例2

2744557306

# logistic回归
# U& x1 o" e; ~5 O; \实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。
) F; j8 E9 f, s/ V% h2 {
R软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。


有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。
* y' O& u& D: K$ G& s
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)有完全相同的结果，但效率却低得多。
+ T+ @/ G/ c4 ~( E5 J
8 q1 a( Y' K; \* J( f9 `$ {& b- ~  V2）二项分布
% e" l8 }7 b" [) [4 n; o( V  S
* \2 |: [  K$ [; ]: ^) I( w6 v2 O- T
logistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。

" u8 O% [! O9 E4 e6 b8 ]- fLogistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
* t- z2 v! _4 }# ~+ E% n实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据
   \" a6 M9 z. [7 ]/ J  J  D  D  H0 r, T: s
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )
   . r3 n/ ~* w$ E1 w
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  7 W0 S1 R8 j. A) q. n6 ~3 i
   
4. 0 N$ k/ Z' K3 M6 A\" q
   
5. #2、建模8 c8 Q& ?4 i3 @' l9 u
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)& T\" j/ S: {8 W6 A
   
7. 
   # ?' [3 l\" B, {0 Z( W
8. #3、模型评估
   ! O3 \- ~; z0 y% R$ k
9. summary(glm.sol)

复制代码

1. ##
   + v% |) ^0 Y! B) t1 T
2. ## Call:
   ' J5 ]* G' {! N- h2 Y& ?
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)
   - b4 c\" W& S! i1 p% u+ w9 V
4. ##
   ! i5 {: F: y) v( O9 c
5. ## Deviance Residuals: . |# ?4 F# s) N& N1 C
   
6. ##       1        2        3        4        5        6  
   ) ^0 j: _3 A6 c) a. N& ~4 {
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
   6 e: E, V% T0 c& S  p  T
8. ##
   8 j1 O& W/ _8 [9 A. `, ~  Q) x
9. ## Coefficients:) v  F. d5 g1 t2 [. C
   
10. ##             Estimate Std. Error z value Pr(>|z|)    : S# G) q\" I! x+ u- c' X- f
    
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***
    8 m1 Q- ?6 x8 ?/ D7 R\" i$ [
12. ## x             1.2459     0.1119   11.13   <2e-16 ***; T7 R% `; I& l: [) f1 @
    
13. ## ---* w9 ?2 q# f9 j- y- ^
    
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1+ c+ ]1 G3 h2 A8 p
    
15. ##
    : @0 W2 J5 j2 J
16. ## (Dispersion parameter for binomial family taken to be 1)
    7 D2 o\" D% Z9 X9 M
17. ## ( O  r5 G% E& i/ Z. b0 N: h
    
18. ##     Null deviance: 250.4866  on 5  degrees of freedom
    0 `; e( Y0 m8 p% a: b- g
19. ## Residual deviance:   9.3526  on 4  degrees of freedom
    % K. e) d/ q, p9 f
20. ## AIC: 34.0938 L0 X3 z: ^, R
    
21. ## * N& u2 _4 I+ t6 d& F4 U; A2 v
    
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率
   0 {- E# D# i3 z3 P/ ~
2. 
   ! `, G% y+ D& c$ P8 j- I$ O
3. #4、预测) l% }, J, m\" v/ l
   
4. pre <- predict(glm.sol, data.frame(x=3.5))
   , t* d) o& v% a* [
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1 5 I$ t4 ]7 u0 U. H; ]
   
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1
   ) ], m+ e4 R- A8 G+ ~/ K2 v
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x
   , I# O, ]& q* [\" N: U
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：. {; w/ p% x+ \
   
2. d <- seq(0, 5, length=100)
   * |% P0 D- s1 O# H
3. pre <- predict(glm.sol, data.frame(x=d))- e\" L. u* W% [  q& h2 l
   
4. p <- exp(pre)/(1+exp(pre))+ C( \; D& ^9 M  C. O) l! V
   
5. norell$y <- norell$success/norell$n
   2 E, }2 r8 I0 K& c9 Z$ L
6. plot(norell$x, norell$y)
   . A8 R# ?6 d6 _4 U0 x
7. lines(d, p)

复制代码

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

复制代码

: W8 [0 ?/ C0 d. v  T. @4 ]