Logistic回归实例2

2744557306

# logistic回归
, P* c% X$ b' \! \+ e1 D3 u实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。

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

有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。

8 Z1 R0 N2 `, h/ t1）正态分布 正态分布族的使用方法： 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回归模型属于广义线性模型。
# d6 H2 A3 i6 O# V4 ]) G8 _7 @$ ?
/ E+ U! z; N* y, i; LLogistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据
   ! U4 ?/ [% R0 s$ e% o
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )\" ^( I' u. _7 z: V
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  
   + x  |8 ?8 B) b) G& }
4. ' W1 m+ j1 F4 H- l* D\" X% l; B: \8 X
   
5. #2、建模% H0 \; \, V* I2 ^
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)- ?' U! r/ A; R. a& n
   
7. 
   # g! T7 ?2 `# w& j# S7 D# ]! G
8. #3、模型评估6 [% ]# G( w  ~$ e, P) p7 m
   
9. summary(glm.sol)

复制代码

1. ## 4 U. F, {& J- z+ L. I2 v9 d2 ]
   
2. ## Call:
   % y2 q9 D4 ~, [# P8 H% `- j
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)$ H\" T5 b9 u7 E* t$ p\" {- ~
   
4. ##
   ! e' O8 V+ \7 B  H$ @/ f\" l
5. ## Deviance Residuals: ( Q7 u: v3 ^: L: l
   
6. ##       1        2        3        4        5        6  
   , l# x6 B; x. A3 ]5 Z) j
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  ! ^# D5 C' y. b5 E! t
   
8. ##
   # z\" L. G# ^# Z
9. ## Coefficients:. h7 y! M\" K0 D
   
10. ##             Estimate Std. Error z value Pr(>|z|)   
    $ m- a- U9 Y: m
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***
    $ E. ]# ]6 u4 V7 J( B9 u\" J6 [
12. ## x             1.2459     0.1119   11.13   <2e-16 ***2 z, J& K/ L( M; h& _% _
    
13. ## ---
    . ?; u2 X6 s/ i8 _
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 19 F7 `) K0 y  A/ ?' A; g
    
15. ##
    1 H# P& g3 l, v( W  ]
16. ## (Dispersion parameter for binomial family taken to be 1)2 w% K5 b% Q9 ~( j; |: \. H; v
    
17. ## 1 ~8 r1 _; J( p. g6 J7 G7 z2 |% [
    
18. ##     Null deviance: 250.4866  on 5  degrees of freedom
    3 H/ b- @8 C( @6 @/ d6 s) Z7 G
19. ## Residual deviance:   9.3526  on 4  degrees of freedom( g- y\" Q  q8 C/ H\" y+ ?
    
20. ## AIC: 34.093
    ) x$ a9 ^  d9 C
21. ## . e0 p% l' k; |5 ?
    
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率
   5 s' d6 Z\" p# @\" h) F; |2 F
2. 
   & ~. n$ ^. |3 }1 B. Y% n6 k
3. #4、预测
   1 V( i$ h3 D9 V$ |% i% z
4. pre <- predict(glm.sol, data.frame(x=3.5))
   / D7 J7 e8 C+ E& r% G) ~
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1
   0 [/ ]* h9 ^9 e  N4 S
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b12 ^& a1 P0 G5 w9 L. T2 l) f
   
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x
   ) t* y# B\" n4 G5 i! @; ?
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：
   # ]9 g& w7 e6 `# ]( G
2. d <- seq(0, 5, length=100)
   + ?( B% ^! S% \: H0 B0 @# k
3. pre <- predict(glm.sol, data.frame(x=d))7 l$ y& E$ U0 h5 ^, R
   
4. p <- exp(pre)/(1+exp(pre))$ B) D8 Q' r\" |
   
5. norell$y <- norell$success/norell$n7 h7 W( y; k! Z7 g2 b9 k# ~0 p# h& N2 m
   
6. plot(norell$x, norell$y)0 I) ^  @' C4 Y: t: y! C0 |1 p
   
7. lines(d, p)

复制代码

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

复制代码