Logistic回归实例2

2744557306

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

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


8 M7 X" V# P5 [" t& Y2 z9 ^- a, g8 k有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。

4 e  M% {: b6 g" L3 M1）正态分布 正态分布族的使用方法： 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）二项分布
; U9 k  ~1 Z( L- V$ L
1 g) h- l! }6 F' l. {; M  _% M  k
logistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。

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

1. #1、加载数据$ B/ ]$ t& w1 }, V7 V$ |
   
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )
   ! Z# x6 |8 {8 t
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  # C. O* w& _& @. |# Y
   
4. ) ?8 p7 S4 b+ ~. i: ^/ u
   
5. #2、建模' A2 Y* G* y) w8 d$ }
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)& ~& y7 r/ R& q
   
7. 7 B  C+ E: J, s( ~
   
8. #3、模型评估
   * a' O5 n$ j' L) ?! d
9. summary(glm.sol)

复制代码

1. ## 7 O; H7 B/ A; m* Y: B4 D
   
2. ## Call:4 r% c6 X, |; p$ g: Z, e' O  D
   
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)
   * T; u\" g/ P+ B, b7 o1 Z
4. ##
   & n% ?6 u  s2 x( E5 M( N
5. ## Deviance Residuals: # x# a8 O$ g7 Y# c) k$ q
   
6. ##       1        2        3        4        5        6  
   # o' M- B5 x7 }) P; z. o1 T
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  ; l$ C4 ^& M! M5 v4 K
   
8. ##
   / g+ w3 b- d# J
9. ## Coefficients:
   ; F5 J+ U% N! m2 J* S+ d! ?
10. ##             Estimate Std. Error z value Pr(>|z|)    . O9 {& N' I( L# q) ~' Z$ y
    
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***9 D3 J6 G: H; e3 f
    
12. ## x             1.2459     0.1119   11.13   <2e-16 ***: u4 f( w. |3 X0 B$ |- |
    
13. ## ---# o1 X7 b* }' T8 O3 p3 ?& h* _6 C
    
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ; I: k2 i4 k! p' w/ V
15. ## - X1 T& z# l0 X
    
16. ## (Dispersion parameter for binomial family taken to be 1)& T, I/ b4 {( Q! W$ n0 T
    
17. ##
    4 \) n9 I  t+ s3 a! ?; [\" D
18. ##     Null deviance: 250.4866  on 5  degrees of freedom
    ) D6 A6 r( a' m* h. \% ^& v+ [; U0 U9 f
19. ## Residual deviance:   9.3526  on 4  degrees of freedom
    6 \' N& I\" _. J+ C# x5 x5 V- c8 g
20. ## AIC: 34.093+ r4 m4 A1 i0 v  g2 g* Z1 f' r
    
21. ##
    # _9 ~! w7 h, D9 _* N$ Z
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率
   ) a% }' }/ J( t+ N$ ]* ~: i6 S\" a' Z
2. 2 H% v+ `, v8 M( X% e% d) p- |
   
3. #4、预测( B\" P# H, _8 n' R
   
4. pre <- predict(glm.sol, data.frame(x=3.5))+ M7 K( }9 J$ M7 p2 A
   
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1
   3 G! d% s$ a3 l$ H' b0 a# ?
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1$ r! k1 l( e  a1 \+ ^
   
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x ; u  b6 ~# L9 V; \. z# ~- [
   
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：1 J\" Y3 p2 E4 v# @! n+ c
   
2. d <- seq(0, 5, length=100)
   ; G* [: G+ R' U
3. pre <- predict(glm.sol, data.frame(x=d))# S8 E* V( n( _8 n5 ^4 }2 @8 b
   
4. p <- exp(pre)/(1+exp(pre))
   ; F& S! \6 c2 l% N
5. norell$y <- norell$success/norell$n
   $ T# N) V$ z\" x' L0 C& U
6. plot(norell$x, norell$y)
   # H6 g8 n5 x  I
7. lines(d, p)

复制代码

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

复制代码

1 r7 {7 Y9 X% j. k$ n+ m3 y