Logistic回归实例2

2744557306

# logistic回归
" n: b. G( z- m: F& y* A实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。

R软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。
4 ~6 a  u; _9 m
1 `1 z* N& A5 c
有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。

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

2 U. K7 A# h( Q& a2）二项分布


logistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。
4 e3 u$ K8 B# |) S+ l' Z. i8 _
+ Q! W. ?+ i1 d" d4 n4 ^" bLogistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
! V  S* Q. H5 W+ c8 x实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据9 u\" w( u\" k& V% h; h& K' t
   
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )6 R& j( p# j' q# W6 O  B4 z, l
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  
   8 T\" S+ E2 o8 |! H1 J
4. 
   2 o8 N3 ?) W4 N- w! O+ [- A+ T
5. #2、建模: `) ]  ^  A6 S; d+ H; ^
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)
   # n& p9 P& U) p. f
7. & ~. o# y6 Q+ j. q9 }
   
8. #3、模型评估5 b' Z6 a; _' o: F/ g% s, D( s
   
9. summary(glm.sol)

复制代码

1. ##
   : v/ M# k0 c3 r' r) o1 a. H
2. ## Call:( R* U! k/ K; l) x
   
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)
   2 D7 s6 J7 R$ Z5 p. |\" r
4. ## ( @: Z\" |8 F* ^( b' B/ n1 ~
   
5. ## Deviance Residuals:   H0 k' B6 w1 W  T8 `3 ~8 g; X
   
6. ##       1        2        3        4        5        6  
   0 |9 M* K2 r) \- @* ?5 \: M
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
   % `& ?8 j- R* X! V6 U2 x
8. ##
   ' @& s9 r6 w, C) k' x4 U! R
9. ## Coefficients:
   5 A  d+ J, d- [% Z7 H5 [9 W0 ?
10. ##             Estimate Std. Error z value Pr(>|z|)    2 a6 h  P% d3 q1 [/ T
    
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***! b* q5 E5 l( l6 a9 ]
    
12. ## x             1.2459     0.1119   11.13   <2e-16 ***; D, V' ^5 ], x9 W0 \' n
    
13. ## ---4 ~$ X/ ?$ I0 u. f. h: P( S0 [
    
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1; X; P7 h8 }8 J) ~+ j0 h! H
    
15. ##
    $ ]  |4 ~+ c7 k# c0 ^0 Y% ~3 N
16. ## (Dispersion parameter for binomial family taken to be 1)\" L/ r! g( W- i3 u4 `* ~( c
    
17. ## 9 s* [' [7 F+ J5 V, J
    
18. ##     Null deviance: 250.4866  on 5  degrees of freedom\" V: E  ~7 j+ t# l, w# C' h
    
19. ## Residual deviance:   9.3526  on 4  degrees of freedom. x2 L\" h4 _' n4 ?( d2 x  O
    
20. ## AIC: 34.093
    ! [) ~5 [+ o. s! K9 i# ~7 Q
21. ## / l4 ?: y% k& ~* L) a# b
    
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率5 o  P, U4 s/ R
   
2. ! U9 }+ m- G. v1 x5 ]# `3 H4 ^9 c# ]& {; Q
   
3. #4、预测
   % @1 v! D6 n: `; I. D' V! T+ Q8 s
4. pre <- predict(glm.sol, data.frame(x=3.5))
   / Z$ S; S. m$ ~* i/ ~
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1
   6 Y2 C3 D' j6 I2 A& s2 h
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1! Z\" J& Z7 y  `9 x3 h\" ~
   
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x
   % X. P7 \0 k- D9 U
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：, k1 s  E, n/ @
   
2. d <- seq(0, 5, length=100)
   6 p* v6 I. u: I0 m
3. pre <- predict(glm.sol, data.frame(x=d))2 d% y+ B6 A9 V) A# V# I5 s# d6 s
   
4. p <- exp(pre)/(1+exp(pre))9 q7 D9 W+ A) |: ?, [- t
   
5. norell$y <- norell$success/norell$n# ?6 v8 l& d# X
   
6. plot(norell$x, norell$y): Y$ f  d& [! r. E; u  l% X
   
7. lines(d, p)

复制代码

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

复制代码

1 S8 V% Q5 ]. {0 `5 R; s