Logistic回归实例2

2744557306

# logistic回归
7 U* _+ k5 u% Z0 J  U8 }5 p" s: Y实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。
% k6 ^  n+ c5 U+ q/ m# R9 w
R软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。
8 F4 i( U% p3 F$ x+ G

( D" m6 H9 Z+ `$ ^) l  a& D( `9 G有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。

! _% _  g+ G; M: w, [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)有完全相同的结果，但效率却低得多。
* ^& T; l8 C5 j6 `" Q9 X+ G: @
2 [3 F% G( @8 g2 y5 c- n1 c4 e& \4 U5 @2）二项分布
, ?' {4 J* p7 M/ O

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

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

1. #1、加载数据% G) g  k0 E+ O4 o% I# o
   
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )9 m6 z3 Y# M3 z+ ]/ u: ?9 O0 [
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  
   * E: y0 k* _0 X9 }; W
4. ; e) m9 s4 L0 L  ^0 l+ o) ?
   
5. #2、建模2 v2 u4 }1 k+ a\" G6 O7 ]2 T
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)7 s4 D% `' p) J2 W: X# W0 Y0 E
   
7. 
   ' c( k: v- k! i9 p6 p+ w
8. #3、模型评估
   : K  k0 \7 u. U$ D' ^( S! ^
9. summary(glm.sol)

复制代码

1. ## & ?$ p; x3 h7 ~/ T
   
2. ## Call:
   : t. W! d/ |3 \9 W8 \
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)- N3 H  p7 l% }% K& T
   
4. ## ! c/ r$ L3 y* B' w, a+ a' y; k\" k0 y2 }
   
5. ## Deviance Residuals:
   % h& F& n' `9 C, @, [) `8 f0 D
6. ##       1        2        3        4        5        6  
   ; w% e7 q5 T/ r+ d% K
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  - p$ |7 j0 c+ m6 f& U+ j
   
8. ##
   $ u/ P1 u' n3 p, N
9. ## Coefficients:
   6 b  g. e7 Z( w: V) l$ b8 r1 k
10. ##             Estimate Std. Error z value Pr(>|z|)   
    + |0 p% W! j3 Y$ |0 K- @% t9 O
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***
    1 l0 q6 n6 a: R- T1 M/ S
12. ## x             1.2459     0.1119   11.13   <2e-16 ***) m4 r% f0 U- N3 A$ A, k
    
13. ## ---) {0 B$ n4 h  W, i
    
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 10 X- v% y# @5 i2 ~0 X  H
    
15. ## 8 `; ?! N5 m2 B& c  L( |\" D* I; e
    
16. ## (Dispersion parameter for binomial family taken to be 1)  Z; g1 ~4 ?6 o' ~1 e* o
    
17. ##
    ; p) e6 W6 m# L
18. ##     Null deviance: 250.4866  on 5  degrees of freedom# p5 H0 W. {* A! e, f
    
19. ## Residual deviance:   9.3526  on 4  degrees of freedom2 k$ I. s, X- o
    
20. ## AIC: 34.093* Y) Z# [0 @0 ]* _. `6 V' ]: l
    
21. ##
    , n  R% }; S3 k8 I5 t# H
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率
   % e) J6 L; ~4 \) a' E# b8 Y8 ~
2. ) {, k- M; ]9 j! G
   
3. #4、预测0 {6 N* f- W2 ]; g& z2 ?
   
4. pre <- predict(glm.sol, data.frame(x=3.5))
   + x2 N/ Q6 s+ X0 z0 a8 e
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1 ' a' [/ u6 o7 C9 ]\" [
   
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b14 }8 y; q\" ]8 D. Y* G7 _
   
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x ( u) P7 r2 S6 `/ A6 {8 D  O) ~
   
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：; M9 e1 Y# a) h6 S6 H8 P
   
2. d <- seq(0, 5, length=100)  h0 C) e6 H9 S) j) W- _
   
3. pre <- predict(glm.sol, data.frame(x=d))
   ; q- A) O& r5 n' p\" I0 n
4. p <- exp(pre)/(1+exp(pre))6 l7 ^8 W6 w% a3 e' W
   
5. norell$y <- norell$success/norell$n5 y* k, }, w! Y9 S+ k
   
6. plot(norell$x, norell$y). }+ M' M& s( f' r( j% C& h' K
   
7. lines(d, p)

复制代码

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

复制代码