Logistic回归实例2

2744557306

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

R软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。
. H7 z& f# }* L' p% }0 }/ Z  a& y
' E7 {4 j; y" b
有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。

8 ~0 L" A! G  p# V9 O3 T# ^6 p1）正态分布 正态分布族的使用方法： 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)有完全相同的结果，但效率却低得多。
/ [. }" J/ i' y% w8 c) w. L
2）二项分布
& o9 l& R% F( O7 H

logistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。

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

1. #1、加载数据  h7 S! o* s! ?2 @
   
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )3 F: ?9 X# K8 R+ ~
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  
   + D* d7 w' i- ]# O
4. 6 ~' B& d9 O% S
   
5. #2、建模+ T+ d0 L' C5 I9 j\" i3 q
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)
   : y+ z3 e/ ]5 O, Z/ j. [) g
7. 
   9 ^6 g9 ]9 U) j4 q# Z/ G3 n
8. #3、模型评估
   ; [% H- b* _0 q, L* F
9. summary(glm.sol)

复制代码

1. ## 4 l1 U0 k. W; ?
   
2. ## Call:
   ! F! g3 }/ d; j  b4 J/ _$ Y7 Y0 w# ~
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)\" R3 ?* d. \  s* i# F
   
4. ## % W6 a$ o3 m% L: u7 v- V( K
   
5. ## Deviance Residuals:
   4 U3 y5 S8 m. s. |+ a6 J! r
6. ##       1        2        3        4        5        6  \" w1 j9 V, ]9 D: r8 T. v8 _
   
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
   # M! R' A9 C\" q$ j! f& j  s
8. ## 6 I* q# x1 {$ r3 _8 u8 c! Q
   
9. ## Coefficients:
   5 I% d7 h* e1 ~
10. ##             Estimate Std. Error z value Pr(>|z|)   
    . y8 P. Z4 ^% X, }
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***) ?7 a0 X/ d+ }$ I! M1 t6 b1 ~- q
    
12. ## x             1.2459     0.1119   11.13   <2e-16 ***
    ( v1 @0 n, }4 n  v5 {/ z/ h5 a
13. ## ---* T! G, I9 j4 M, e
    
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    & r( y3 l- U\" N) }% u7 l
15. ## 0 ?$ x) C& U* N1 O/ @! L% F
    
16. ## (Dispersion parameter for binomial family taken to be 1). ^; ^  S# _, n
    
17. ## ( t. f( \9 A# g! r( Z
    
18. ##     Null deviance: 250.4866  on 5  degrees of freedom8 J. j1 z\" U$ O2 [9 b
    
19. ## Residual deviance:   9.3526  on 4  degrees of freedom: ~# A& L) `& o' @% V! d4 r) f1 b2 g
    
20. ## AIC: 34.093
    5 m! @9 t( Z$ K% d% q
21. ##
    ) U5 M\" v\" Q; P+ F- e9 N
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率/ C' K/ l8 O! S% a( o/ X- T
   
2. 
   & y0 j/ l' a1 T' E8 {' S
3. #4、预测8 \$ A, u: A* b8 J
   
4. pre <- predict(glm.sol, data.frame(x=3.5))9 Y% S. N\" t5 t$ g6 b7 e
   
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1
   / Q  _7 ^. J6 H5 D7 O5 v: f7 g
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1: F7 }# j+ b* n- G$ J* }' B
   
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x . I/ @7 o* V3 a8 P, v- g1 u
   
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：
   \" v\" ]2 H0 v8 }, B* P/ ?% v
2. d <- seq(0, 5, length=100)
   8 _2 f: B: j4 Z/ z+ I2 t
3. pre <- predict(glm.sol, data.frame(x=d))! E  K( L/ I: ^7 E\" G
   
4. p <- exp(pre)/(1+exp(pre))$ x- D( o1 y& j# i+ K  Y
   
5. norell$y <- norell$success/norell$n
     X) H: B$ n\" c5 A7 |
6. plot(norell$x, norell$y)% H! y\" {1 A. t; b3 N+ }
   
7. lines(d, p)

复制代码

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

复制代码