Logistic回归实例2

2744557306

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

R软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。
. _. e9 P8 v- Z" l

有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。
' }) [" a$ t4 ^' _; T8 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）二项分布


logistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。
" V! @* d; I- p" m9 S. W- L* O) L
) O& A' r6 S& B' y* M" ?- k# DLogistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
. r0 w1 Z* {1 T* w  F实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据. G, q- o/ [% G
   
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )% T5 {% {: M- x9 L\" H
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  2 |+ J( @, G8 d( c: A
   
4. 
   4 y0 ^$ t2 Y: F0 P0 D
5. #2、建模
   & T4 H! J4 L4 R/ O( \, T% U
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)
   % y, e\" c\" T& B( r; B; W4 n
7. ) ?, N! v5 d( l- ]& c
   
8. #3、模型评估
   4 ]- @; P1 s$ D3 i\" Z& \! i' {
9. summary(glm.sol)

复制代码

1. ## 4 D\" j! }. X+ _4 M9 B
   
2. ## Call:
   \" A/ U: P! L' V- [& t3 @$ t4 C
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)
   \" g# t: k' {% a+ c\" }
4. ## 1 Z\" g  {8 M7 s3 J) B, H
   
5. ## Deviance Residuals: \" E0 V# u5 x& t2 ~3 ^
   
6. ##       1        2        3        4        5        6  2 w9 p\" s6 V  u0 n: D1 H0 }$ U4 `
   
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
   % O7 I7 e7 C7 i& c
8. ## $ k- @9 }3 [- @
   
9. ## Coefficients:
   8 |. A3 y6 L# Q\" x
10. ##             Estimate Std. Error z value Pr(>|z|)    5 ~4 |9 ~6 }9 t  j9 w6 _\" d
    
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***9 u) Z9 g9 I5 J6 p! ~
    
12. ## x             1.2459     0.1119   11.13   <2e-16 ***8 J\" L+ v/ K; w' e- q
    
13. ## ---6 B; E+ j1 J3 L\" u( d+ s( `& n1 N! a
    
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    3 x# _' h- n- I! Y, p! c- G' K
15. ## 2 i. B1 j% L  X4 L9 E$ s
    
16. ## (Dispersion parameter for binomial family taken to be 1)1 ]0 G# `0 O9 f6 Y4 m- x
    
17. ##
    ; H( {% ~/ i3 A, c1 m* r
18. ##     Null deviance: 250.4866  on 5  degrees of freedom9 w9 h0 K5 ?% {3 A  B' W
    
19. ## Residual deviance:   9.3526  on 4  degrees of freedom
    5 q0 k3 n8 `; u\" [; E- J. K  H1 w
20. ## AIC: 34.093# }; C0 l& e0 N2 y; y
    
21. ## % Q) P# y: t2 m8 w: }, o
    
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率+ Y9 @3 }! U# M3 u8 M
   
2. 8 u/ s! C8 }2 g+ M  W/ o
   
3. #4、预测
   ) u- q% @0 \5 J; i- }  T
4. pre <- predict(glm.sol, data.frame(x=3.5))
   / `7 R+ Y* g  U- l3 R  k
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1 & ^# L) v! P. o8 w5 R) M
   
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1
   1 Z1 l; @7 o\" Q# u/ Z
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x
   8 S! F4 n, m) `\" C* A
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：# N5 K\" G6 j\" ?* ?+ H# s. @3 t
   
2. d <- seq(0, 5, length=100)
     X# c* R# t+ L$ `8 {
3. pre <- predict(glm.sol, data.frame(x=d))
     X: I: B  C: \3 e( u
4. p <- exp(pre)/(1+exp(pre))
   5 r  k3 S2 x3 L# P3 t7 A! h
5. norell$y <- norell$success/norell$n0 a4 Z% ~4 h) i3 c4 u( K' A
   
6. plot(norell$x, norell$y)
   ) R; A7 |9 P% q5 B, @
7. lines(d, p)

复制代码

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

复制代码

& }& Y: G1 [. `) h) {