Logistic回归实例2

2744557306

# logistic回归
% I+ x. K2 h9 U( v# h6 I) z实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。

2 J5 b1 W, }3 R. \* n" |, xR软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。
3 d* s# o0 P6 T% a) j  H- C" m

7 Y: Q8 b9 p; p4 R3 i1 W# q8 N+ |有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。

) Q; f9 u( t, o5 o1）正态分布 正态分布族的使用方法： 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)有完全相同的结果，但效率却低得多。

! f3 V, o% ~; b: v% n2）二项分布
' n+ J" Q, F1 F" C2 m! @1 [
) c& f8 M0 f& i5 r, M, _: U
: E! o- P) i) llogistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。
6 |& E; b$ t  u0 f. l1 h
/ d* n* x& ~7 @  G# y* vLogistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据- G0 ?( q. y\" u1 G
   
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )
   & [! H3 u( i8 d% I
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  / ?  u) M4 n5 z5 S( |) Q; G
   
4. 
   * Z& k$ M  v\" i4 _+ _
5. #2、建模2 C: O% H1 z5 c8 \4 p
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)0 }& J# j% Y; i0 O
   
7. 2 e* `8 \0 g3 n- F7 s3 H, f9 D3 d* s
   
8. #3、模型评估! O6 p$ ^: [3 l* y8 I# T
   
9. summary(glm.sol)

复制代码

1. ##
   * Y# z) {\" {' {5 D8 ?( N. x1 {
2. ## Call:& \/ V. ^, Q4 A6 A& L: X
   
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)\" @, W/ ^3 P1 m. b: s7 b# d
   
4. ##
   6 |# d/ A' c5 G4 B: `2 }
5. ## Deviance Residuals:
   0 y5 l$ q( k2 X) O
6. ##       1        2        3        4        5        6  
   \" ~9 j/ d2 x, V' B* {
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
   & J, D+ M' |1 ]8 G2 m- G* n
8. ##
   1 ~- ^4 z3 B' K- v% X/ q9 v
9. ## Coefficients:
   , s- N2 _3 u% e2 s8 M3 b
10. ##             Estimate Std. Error z value Pr(>|z|)    0 U\" l, f7 R! r7 v- L; n  q
    
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***
    # Y% y\" s7 @0 _
12. ## x             1.2459     0.1119   11.13   <2e-16 ***+ _0 ~' _% K+ q1 J: [
    
13. ## ---+ k4 b6 n) t+ s  L9 ~& Q
    
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    0 k0 c9 C9 d3 E5 [
15. ## . t) [  n' X/ ^( u, b5 n! D
    
16. ## (Dispersion parameter for binomial family taken to be 1)
    1 ^* E# d( t6 I
17. ## 1 Q. I' B) Z# V6 G
    
18. ##     Null deviance: 250.4866  on 5  degrees of freedom  \' i4 U3 z8 Q\" V, {* G
    
19. ## Residual deviance:   9.3526  on 4  degrees of freedom
    6 h* v& ^! W9 M7 a+ ^
20. ## AIC: 34.093
    4 D0 y% |- d: s0 j. S
21. ##
    \" k! j5 ~# r) {- k& b
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率9 b* l) g6 z) b$ M: F+ ^, D
   
2. 
   & Q0 }: x- G: L; b! r2 b
3. #4、预测, }+ |  n6 z( m! m( I
   
4. pre <- predict(glm.sol, data.frame(x=3.5))6 {. F- j0 F& C3 {$ J( j; }
   
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1 ( n& f) _( s' Z' Y& h
   
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1. p) X. e. X. ?' q9 n6 c
   
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x
   0 l; f. C- s  s& z. Y' |\" J
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：
   \" v+ l9 m  L/ F8 ?2 u
2. d <- seq(0, 5, length=100), ^# c( R2 V2 n$ Q, V' M) u
   
3. pre <- predict(glm.sol, data.frame(x=d))# }/ f3 L9 c( Q/ k8 W: W, u\" K
   
4. p <- exp(pre)/(1+exp(pre))
   & Y# J& B0 R: Y/ r# ^' r2 s  [
5. norell$y <- norell$success/norell$n  }# j7 l\" Z8 A  b) e5 |7 p0 Y
   
6. plot(norell$x, norell$y)
   : H6 J3 z. ]- m* N3 Y% n' V8 O
7. lines(d, p)

复制代码

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

复制代码

8 k5 W' P. U9 K! F) R. a