Logistic回归实例2

2744557306

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

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


7 y' }4 {; J! E1 v5 ]) J$ W% a5 y% D有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。

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)有完全相同的结果，但效率却低得多。

, k4 ~1 u' `, u, Z: h6 Q$ e2）二项分布
7 j7 n; U3 N7 H# K0 x, B
" L: C& h. s1 J% z$ h- F. R3 d
logistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。

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

1. #1、加载数据  u$ z3 F' o% @- S  T3 x
   
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )/ |( j& ?' u; {: E7 V
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  
   $ S3 T2 a; b( T) ]: \. P  s
4. ( a* e- D6 D, F1 `  ^( }* L7 m
   
5. #2、建模3 v1 y6 h# b& E: `* ~5 }; ?! i: C
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)
   - {7 l) q0 f, p2 H8 u5 ^# \/ O% [6 s
7. 6 r8 [8 Y6 j0 q9 L
   
8. #3、模型评估
   7 E$ C5 b\" ?* z
9. summary(glm.sol)

复制代码

1. ## 0 Y/ E6 d; N. w\" o& h/ N) L
   
2. ## Call:
   . d) m7 {1 t\" ~- p, W% w
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)' P. X7 P4 t4 m! G4 _
   
4. ##
   , r2 M! l8 Y' w1 O
5. ## Deviance Residuals: 8 n4 f8 t- P# c! T% ~( c7 K
   
6. ##       1        2        3        4        5        6  
   6 ~' v8 E4 K/ S& ?% q( q
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  + [# n7 N5 ]' l
   
8. ## * R, A' {5 b2 \/ c/ \
   
9. ## Coefficients:
   1 |6 F/ \  [5 f5 M4 c
10. ##             Estimate Std. Error z value Pr(>|z|)    # ]+ Y3 T5 w3 s2 E8 i
    
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 **** N8 {- W4 U) x, }9 J- j* _
    
12. ## x             1.2459     0.1119   11.13   <2e-16 ***+ \! S* h+ U/ o( N2 K
    
13. ## ---
    ' |+ C9 o5 n7 u9 g
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1/ k/ z8 H7 y  R! f: g
    
15. ## - l0 X9 I: o: f9 }9 z: Y
    
16. ## (Dispersion parameter for binomial family taken to be 1)% H5 y* A9 d4 O0 Z: i
    
17. ## ( D5 K# _- K+ H$ @- J6 c
    
18. ##     Null deviance: 250.4866  on 5  degrees of freedom
    $ m0 D7 i3 T8 ?8 H% \0 }6 E2 |
19. ## Residual deviance:   9.3526  on 4  degrees of freedom
    4 d* c/ m+ @3 m$ ]8 L( I( s
20. ## AIC: 34.093
    * A( _; a- }% W' s# U
21. ##
    / s  v9 I\" _! Y7 D! \
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率% E- i; F1 N( A
   
2. 
   8 V2 ]0 n* _2 z( ]+ I7 o
3. #4、预测
     `, ]\" p3 o' ]; @: ?1 Z
4. pre <- predict(glm.sol, data.frame(x=3.5))( W1 x# W( w# V- L/ Y
   
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1
   $ E+ ?& P9 u: T* k
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1
   : }) s6 `6 G9 n/ M9 e5 w
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x
   $ x* M9 i( m/ [( U9 n
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：
   2 O; _1 l5 B6 H
2. d <- seq(0, 5, length=100)7 a& _* I/ p# l2 D. }4 `' h8 w# X
   
3. pre <- predict(glm.sol, data.frame(x=d))
   3 N# E* k( }3 m1 v. s\" t- p( r
4. p <- exp(pre)/(1+exp(pre))# b3 ^8 U9 k8 ^- M\" R
   
5. norell$y <- norell$success/norell$n. {0 [- x! o* o- h4 f3 ?
   
6. plot(norell$x, norell$y)/ J) [/ I: |! j8 Q5 o
   
7. lines(d, p)

复制代码

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

复制代码