Logistic回归实例2

2744557306

# logistic回归
实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。
' h. l0 H) `) Y0 ]6 @- n' a$ G% i( R
( K& b  L) X0 ~9 G; X* [R软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。
- [/ d% x9 h1 P# V* K

- o1 x8 t. Q) G6 z& I有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。
* k" Z* g( ]" L; R( f4 ^. [
' r$ Z3 d/ U4 D1）正态分布 正态分布族的使用方法： 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）二项分布


1 y4 l, z: f  K$ `8 C# Tlogistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。
5 i% W7 a( \/ I* o  a) a5 W8 ?
Logistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据
   7 {: m; t7 _* d( \
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )
   : n: L+ F( P- j1 J! h0 j# F0 |
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  2 O7 s. Z  Y3 f& r
   
4. 
   ' f  @& X% ?4 a
5. #2、建模
   0 [9 I  L/ I% H8 B
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)1 H1 k\" Q2 u$ e& C3 d
   
7. 
   ( z0 F; S7 S  M* n
8. #3、模型评估
   . ]0 }* z8 S# i$ y3 g% h& m\" h
9. summary(glm.sol)

复制代码

1. ##
   3 t; Y# C/ i. O- p; w2 H  b
2. ## Call:
   $ v, _5 h* {\" d% I6 ^# Y
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)* S$ w1 Y( J0 @4 d: Z/ j
   
4. ## ( H1 Z* Z& `+ h  |) V
   
5. ## Deviance Residuals: * e# J3 d7 P, b: _5 C( t( d\" ?
   
6. ##       1        2        3        4        5        6  % D* f0 d# X3 L# X  D; F  u- k
   
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
   3 N, g% @. n, `- V! j
8. ## + A' N$ I/ G! _5 _! P0 |\" V
   
9. ## Coefficients:, j5 X8 y\" M3 h# Q\" E
   
10. ##             Estimate Std. Error z value Pr(>|z|)   
    7 e$ c6 r) \3 f' }/ w
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***
    7 O/ X1 o! u5 o/ m
12. ## x             1.2459     0.1119   11.13   <2e-16 ***% j0 f2 ~5 W& T
    
13. ## ---1 O3 b9 @, B+ D6 D
    
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    + [! d  Y5 k/ D
15. ## : @\" ?! k$ g3 H/ }
    
16. ## (Dispersion parameter for binomial family taken to be 1)
    . I/ Q  V' s. r8 m# h+ l8 }& B2 J- ~
17. ## 1 m' W& ^! H2 s! ?' m& Q. O* N- }. u
    
18. ##     Null deviance: 250.4866  on 5  degrees of freedom
    3 d& q& X2 l0 b5 p3 Y. M
19. ## Residual deviance:   9.3526  on 4  degrees of freedom* b) c! e  Y# }7 I5 a# j; r1 l. C7 Z
    
20. ## AIC: 34.0933 w0 R8 W; h\" ]: ^: p4 z* g3 }
    
21. ##
    ( n3 ]& y# W$ I  u* _
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率
   + ~+ E3 K7 p$ g' \
2. / C0 p$ Z/ M/ |( c6 z  n
   
3. #4、预测) E: U1 V, ~! d
   
4. pre <- predict(glm.sol, data.frame(x=3.5))' s* O0 `8 t% a5 Y\" K4 w; }
   
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1
   1 P% F, Z0 U; v6 B
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1, z4 ]) K( d! Q+ O
   
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x
   # j  {9 J, ?; |  Q8 V7 @
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：
   $ B! b3 `% Y9 K/ R8 A8 S
2. d <- seq(0, 5, length=100)+ x  h. k# ~. M2 ]! `* d. P
   
3. pre <- predict(glm.sol, data.frame(x=d))
   ) q% K) S: [6 ~* R2 w  p6 @
4. p <- exp(pre)/(1+exp(pre))
   / a& e6 N3 \, N# ?+ b
5. norell$y <- norell$success/norell$n5 I* G; a- }8 x8 J
   
6. plot(norell$x, norell$y)1 [1 x8 C) K/ |7 y! G5 Z
   
7. lines(d, p)

复制代码

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

复制代码