Logistic回归实例2

2744557306

# logistic回归
实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。
  `( Y) S0 G% s1 [7 z
R软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。

3 c" R3 v8 @6 w- @# V: H
: r1 E& j( |  U9 _( J1 u有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。
: F5 ]$ l( s" s: s, H
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)有完全相同的结果，但效率却低得多。

* @7 M# R" Q7 N) n2）二项分布
2 {9 W* O9 L8 U* N0 m$ J
: S+ z) h# d3 \( @
logistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。

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

1. #1、加载数据
   3 t* q2 p2 h: }8 R
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )1 ]8 A6 s6 }0 c$ i+ q, ^5 \
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  
   + T( O( y1 p$ ?. y4 b* t, h- `
4. 
   \" x8 I\" a% W' T& G8 ]
5. #2、建模
   3 u& v6 `* t3 W* y
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)
   \" b* R) e$ T7 i- n
7. 
   ; F  T  o- a+ H; S* Y
8. #3、模型评估6 H; n( M9 c) U4 {
   
9. summary(glm.sol)

复制代码

1. ## 3 l' J, A: x. L! ~5 U4 P
   
2. ## Call:* j! Q1 Z% k% q0 |# X
   
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)- Z5 K7 L+ z/ r# m$ ^& R
   
4. ##
   \" |, C' V: m: c* p  B
5. ## Deviance Residuals: / c  @( d  m\" @/ L# v
   
6. ##       1        2        3        4        5        6  
   * i+ P\" w& W! H$ |- }
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  / m2 N& o- n4 }  B
   
8. ##
   6 X) W\" h) L* R& p1 Z
9. ## Coefficients:9 `0 r/ M. m0 W5 U' ]* [  |
   
10. ##             Estimate Std. Error z value Pr(>|z|)   
    ' [8 v+ p- e% }$ U
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***
    6 I$ c, v! o  c3 K
12. ## x             1.2459     0.1119   11.13   <2e-16 ***
    ; e, j! k% S' \8 U' o' J/ f
13. ## ---7 k0 U( \& N% I3 H$ \2 x# ]
    
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ! s: _, |3 q' c4 z# J7 G
15. ##
    * s) Z3 e4 o# |! K8 ?5 f+ m1 ]* A: ]( D- r
16. ## (Dispersion parameter for binomial family taken to be 1)+ a7 l; Q5 H\" b7 E! c4 s
    
17. ##
    + f% I. c( d+ y( `. N5 N8 Z- I
18. ##     Null deviance: 250.4866  on 5  degrees of freedom\" d. m; u, g# A( K
    
19. ## Residual deviance:   9.3526  on 4  degrees of freedom
    % x7 t, @( j% s+ B5 t
20. ## AIC: 34.093
    \" O' k+ \7 `$ B; o8 w, n/ {7 w
21. ## ( q/ F1 ^8 v\" X' L8 W. I\" K
    
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率9 J\" g2 ?3 |' P' V( D
   
2. 
   ' w4 b5 H! \$ i4 V$ `\" E6 Y
3. #4、预测
   7 o' p9 B9 }: `, G( y- N
4. pre <- predict(glm.sol, data.frame(x=3.5))% C: |\" ^4 S8 }& @% Y
   
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1
   \" s' Z# {  ^* n4 L8 g
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1
   2 i6 `0 A& C+ Q7 P5 ^
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x
   # ^& s' i5 M- o+ c
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：  [  `6 a# V5 j1 O, }. P
   
2. d <- seq(0, 5, length=100)3 ]& r\" N+ B9 [  I/ q5 I4 b
   
3. pre <- predict(glm.sol, data.frame(x=d))5 s+ X' i9 [0 |3 O5 z
   
4. p <- exp(pre)/(1+exp(pre)); q& D: b4 B9 B& c3 v/ L7 Z+ j
   
5. norell$y <- norell$success/norell$n
   / Z6 v0 k$ L8 S4 i
6. plot(norell$x, norell$y)# B2 h\" D- h! g
   
7. lines(d, p)

复制代码

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

复制代码