Logistic回归实例2

2744557306

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

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


有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。
) _- W* S- @; H! Z# R
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)有完全相同的结果，但效率却低得多。
: b* D8 V. A6 I0 S8 {
2）二项分布


logistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。

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

1. #1、加载数据# z3 o, s7 U; l, q
   
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )# [' ]9 Q4 K. Y+ y# W7 J
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  * Y8 P1 c8 f/ d+ R( k' j, P
   
4. 
   2 ?7 |: I6 r. b! X8 V3 j- f
5. #2、建模: v( n* s3 R) a) X# f0 S; U8 [0 K
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)
   + i5 r+ L% s( M; R
7. 
     K- `. m* r( w8 z; J  L0 u
8. #3、模型评估
   ) n( I\" \, d( A4 t. D* y
9. summary(glm.sol)

复制代码

1. ##
   % x8 ^) J& p$ i: Q( v  k) U# q
2. ## Call:
   ) a, d# X% |1 v4 a$ G
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)
   ' r% n5 M) ^. i$ S) P% C
4. ##
   , d0 P# |$ G* o
5. ## Deviance Residuals:   `! m2 Q& \9 h# ]; ~
   
6. ##       1        2        3        4        5        6  ) ?2 b% }7 F\" T, |0 f3 ~
   
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  . J  }# M1 U! N4 l: G. J4 L
   
8. ## % M\" r. q- v) o. t/ ~5 ]8 ?
   
9. ## Coefficients:
   , C1 q, r( t3 _6 b) m, D* y( I; j
10. ##             Estimate Std. Error z value Pr(>|z|)    \" t5 t; [\" z) A' x+ P6 x) T0 X/ P. N* m
    
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***1 _. \  N+ T5 V3 D! I\" |! ]% F
    
12. ## x             1.2459     0.1119   11.13   <2e-16 ***' A0 b7 l( F( t( _
    
13. ## ---
    . X  v! B- }, P0 F: }
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  T9 }: E: o\" q/ F; m
    
15. ## ' W3 {: O, X- D9 \
    
16. ## (Dispersion parameter for binomial family taken to be 1)) J! ~5 v0 |- A2 P8 Z# V) P; u- @
    
17. ##
    2 E4 Z9 r; k. j5 h\" N
18. ##     Null deviance: 250.4866  on 5  degrees of freedom) E4 m\" O* C* D5 W* v1 p2 E6 Z
    
19. ## Residual deviance:   9.3526  on 4  degrees of freedom. V$ S) n3 s2 U\" ~2 _& T
    
20. ## AIC: 34.093\" s/ w+ W& q8 f9 u* g  E- e
    
21. ##
    ) M7 G5 l# M- Y# |' K
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率
   ' m( i+ f. i5 \6 |
2. - V2 n& T& n4 O  o7 d/ j
   
3. #4、预测
   * x\" b9 e( D% o1 i( d( a
4. pre <- predict(glm.sol, data.frame(x=3.5))! H0 F$ m7 A: I+ }5 e
   
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1 ' H  ?8 L4 ^, `8 ?3 f8 k# ~% e4 X0 p
   
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1
   4 h4 {) g  o. f3 [\" P
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x
   . T! @5 ^+ X' H! j5 |
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：7 \8 I2 x9 R8 `% a5 K
   
2. d <- seq(0, 5, length=100): R1 y% k3 f\" Y- C
   
3. pre <- predict(glm.sol, data.frame(x=d))+ `* A9 `6 t\" {$ F7 Q8 H( _
   
4. p <- exp(pre)/(1+exp(pre))
   $ _+ g+ T! k% G& T2 Y* I( g0 E
5. norell$y <- norell$success/norell$n4 t  m+ j$ C( m. ~0 d1 f: X
   
6. plot(norell$x, norell$y)
   ( ~5 Z( I, K  [9 h, [
7. lines(d, p)

复制代码

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

复制代码