Logistic回归实例2

2744557306

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

$ K1 c( A, Z: L& sR软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。
6 F' F6 \5 o" o5 l! W* l

有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。
# E2 D0 q" ]: Q, k% c$ ?
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 L; x3 v7 X  F5 k# m  k6 \% |2）二项分布
" c; _( P8 v7 }) _/ n
$ C2 c. B# \1 }. f2 w2 _
logistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。
( X0 L/ @+ g6 ]8 N$ R- O& |
  {4 a. m, k) }% p' `% G: LLogistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
& m# R# l8 ~9 `+ I# u) E实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据' @) K7 `1 j, X, f
   
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )
   ) \( {  }- w9 [2 g+ n
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  
   / T9 p& \. }8 W: b1 X# Q: s
4. - g5 u6 T9 o2 ^- `4 h5 ~
   
5. #2、建模8 k4 A$ ?+ Q\" Z6 m& s$ I6 W
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)
   : K  R) D$ \8 J* d
7. 4 _, D& _, x4 p
   
8. #3、模型评估- l6 S+ ?( h+ `/ G8 K
   
9. summary(glm.sol)

复制代码

1. ## 3 Q+ [4 l6 w. \' F  D  z3 l; i
   
2. ## Call:& z7 U1 R! ?# e  N* W9 v$ L
   
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell); X* G( j4 R8 P: ~: b, ^
   
4. ## $ u: c  [' ?: k
   
5. ## Deviance Residuals:
   3 \' B5 B6 D. n( ]
6. ##       1        2        3        4        5        6  2 }% a; P+ p1 i
   
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
   ' n1 E% _4 `& h. S
8. ## 2 }7 J8 J2 _. O
   
9. ## Coefficients:9 D& `( }) ]4 E  v
   
10. ##             Estimate Std. Error z value Pr(>|z|)    % A# `% u% S  b9 A; m
    
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***
    : K% Z7 Q* A& `- C+ {; a; M3 r
12. ## x             1.2459     0.1119   11.13   <2e-16 ***, f; D\" J! V' D& I
    
13. ## ---- L) Z' Z5 {, w+ M
    
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1' i8 @% U% p\" p- |! B
    
15. ##
    , V) l: a6 G+ b$ M) N
16. ## (Dispersion parameter for binomial family taken to be 1)1 y8 z9 _0 i$ w7 G- v
    
17. ##
    5 Y( T$ h% e( A5 l
18. ##     Null deviance: 250.4866  on 5  degrees of freedom
    5 V6 m  s9 K- ^# X8 f5 M' G
19. ## Residual deviance:   9.3526  on 4  degrees of freedom0 g2 h% h1 {- w/ e# b. W8 F$ }8 B# y5 N/ U
    
20. ## AIC: 34.093& n2 M: c. N2 P
    
21. ## 0 P- w4 g, x7 @; m
    
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率; K8 B: c5 m) M8 H; i
   
2. 
   5 ~3 n+ V/ z4 u6 B! A
3. #4、预测  h% N( b4 T3 z8 E8 H
   
4. pre <- predict(glm.sol, data.frame(x=3.5))! Q; b5 b( q1 ^9 X% O; C
   
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1
   ! A& g+ F. u\" R8 i$ i8 Y8 z4 ?: Y7 G
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1
   ; v% _2 \) y- f: {, S+ [
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x
   0 C$ u; g4 {, }0 t+ B: ?\" Y
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：
   % F8 P! o$ j) R$ @- Y
2. d <- seq(0, 5, length=100)
   9 C) B0 N: E+ t, g( J) C+ U4 v$ M
3. pre <- predict(glm.sol, data.frame(x=d))& x9 u& p$ H! q+ \) Y
   
4. p <- exp(pre)/(1+exp(pre))
   9 H5 ?4 V) B+ E1 s& s
5. norell$y <- norell$success/norell$n
   # R8 @9 @0 A7 Y; k
6. plot(norell$x, norell$y)/ _# r8 @. F( I
   
7. lines(d, p)

复制代码

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

复制代码

- ~$ b; k9 s9 Z8 O: A) p