Logistic回归实例2

2744557306

# logistic回归
6 z- }$ @! k. B% }- R6 S3 f实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。
7 g' h- E4 O; I, a/ u5 h5 [
R软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。
, R) i( s# i2 n! s" ~" k. l

5 j! E! K6 H( g有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。
' ~2 j3 X8 G& M: z
8 H6 V* S0 e8 j1 P( z0 u) ?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)有完全相同的结果，但效率却低得多。
% y! I5 l! ^, W$ [
2）二项分布
. h' g9 t2 _3 p7 T
- m0 h: j" ~: k) W: b
logistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。
; [& b! V' k- ]0 w: B
Logistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据
   6 `* R% e. l' p, x  G
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )0 o' D- R2 c% g2 }5 H( D
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  
   ( T4 M; g! E9 y/ f; ]6 O
4. 
   / [8 `- G. t/ u
5. #2、建模& |1 J  q\" X3 j# l  m3 X
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)* s2 ^5 O& i9 K9 O: |\" l8 m8 j
   
7. 
   $ O  Y3 D: J: G4 U' i' e
8. #3、模型评估
   * x: j$ k* g5 [* o- H
9. summary(glm.sol)

复制代码

1. ##
   + N5 Y4 I& e$ m# f( Y5 v
2. ## Call:  Y2 Q+ f& T; h' J+ F
   
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)3 g3 @, p\" M: c% B, x/ Q) g
   
4. ## , ?* n# `3 g\" Y! S
   
5. ## Deviance Residuals:
   5 g7 {6 m% q: q, D* \1 R8 O; L0 d
6. ##       1        2        3        4        5        6  + \( J) _0 G\" z# G+ i/ T1 E
   
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  ) F: Q$ m- @# r5 i
   
8. ## 3 h+ z1 H8 h- _3 p. O; [# |' @
   
9. ## Coefficients:6 b: F% D% S\" I6 ?' E: z* D4 A
   
10. ##             Estimate Std. Error z value Pr(>|z|)   
    2 ~: z3 N$ s# `5 p. N' L
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***
    1 K, D9 V0 M- O7 N
12. ## x             1.2459     0.1119   11.13   <2e-16 ***
    ' W( w+ f9 D5 f% b% P0 P) G
13. ## ---/ f$ p+ X2 e' T2 [3 V9 e
    
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 11 b5 _: r% u& x* z& h
    
15. ##
    ) G4 S0 {% M4 A
16. ## (Dispersion parameter for binomial family taken to be 1)9 O3 ^0 o* G# J( w  [. v
    
17. ##
    \" b7 T  G3 p8 j; ?# ~4 s0 R# I! i
18. ##     Null deviance: 250.4866  on 5  degrees of freedom& o$ d7 X6 p! O* x
    
19. ## Residual deviance:   9.3526  on 4  degrees of freedom\" A% h! Z5 R- _/ J' e, k
    
20. ## AIC: 34.093
    ! C; G1 ?. l8 a3 X* T$ b& R
21. ##
    * F$ Y! t6 Z5 L
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率
   & n8 E% E4 T) \  Y
2. $ f  C# ^3 P. E: v5 ~2 ~1 V5 n
   
3. #4、预测' Y2 h* b* z8 ^6 z. L
   
4. pre <- predict(glm.sol, data.frame(x=3.5)); ]) z( k! `9 A/ o+ `
   
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1
   3 K8 k/ t' W' m+ Z/ I
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1
   8 N$ E9 S& `) G8 N
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x : b  d3 n' D5 F
   
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：/ r6 z  t$ O% \, f7 C  J
   
2. d <- seq(0, 5, length=100)
   ' T! }2 S& y/ }9 @- ?% W
3. pre <- predict(glm.sol, data.frame(x=d))# L3 z  k7 x1 }9 ?
   
4. p <- exp(pre)/(1+exp(pre))
     `8 m2 T: O0 h- |0 @; Q
5. norell$y <- norell$success/norell$n
   2 |1 d! Y: v* {$ y4 ^
6. plot(norell$x, norell$y)
   * m5 t! W7 w+ @* u' s  S
7. lines(d, p)

复制代码

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

复制代码