Logistic回归实例2

2744557306

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

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

有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。
. r3 K$ U' l0 v0 O0 H6 o4 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)有完全相同的结果，但效率却低得多。

5 T7 e5 q! `9 T9 w% {& |2）二项分布
. w+ Y# K3 I3 G1 Y+ }
5 v! G; P) q9 d* [
" e( u" p. O7 x5 X9 x6 o1 llogistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。
7 w7 Z3 h# i7 g; s, k1 s
Logistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据& J+ U7 O+ [\" r! g; T# i; g
   
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )2 O$ V9 y7 }( ?/ J9 _
   
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  
   + }6 v2 h+ M2 ?1 j
4. # D3 A% _- {/ \\" L! c
   
5. #2、建模
   / {\" a1 Q! A% n* K! h
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)  E7 F, Z/ i. T, u9 [
   
7. - X; J. x: l) m3 g* Q2 I; L: u3 D
   
8. #3、模型评估
   0 s0 K2 T- I$ t- t
9. summary(glm.sol)

复制代码

1. ##
   7 p/ i2 Y1 q) k7 n0 \: l, v1 d
2. ## Call:
   * U3 Y5 N6 J: R
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)4 M& W2 x' f) E  q( w\" l
   
4. ##
   9 j& J# W, e. d
5. ## Deviance Residuals:
   9 ?3 v  U6 k7 F8 Z$ }
6. ##       1        2        3        4        5        6  
   8 K6 u9 u& g0 X1 d# J7 z0 F, \
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
   4 v: w7 W) y+ j
8. ##
   $ A. S7 }' C- e! G. j6 T
9. ## Coefficients:4 E1 W3 ?/ i4 _* I0 h
   
10. ##             Estimate Std. Error z value Pr(>|z|)   
    3 Z, x5 D' V# M! H\" N3 F/ Y  U, C\" Q
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***
    / z, P; ]& n7 B- c% G) a9 ~; {
12. ## x             1.2459     0.1119   11.13   <2e-16 ***2 Y8 B$ w; M2 o' T6 c
    
13. ## ---
    2 p2 p0 ?2 V6 e4 \) Z$ i# T
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    + Y6 H9 R- l6 P\" C; r2 `5 h3 M6 O5 G
15. ## + I( b\" H# l\" z) Y  g$ c
    
16. ## (Dispersion parameter for binomial family taken to be 1)/ @4 G) C; b5 R
    
17. ##
    % j0 D  }7 I0 w- K  l8 S
18. ##     Null deviance: 250.4866  on 5  degrees of freedom- `$ A6 J4 p8 i8 G- J! S
    
19. ## Residual deviance:   9.3526  on 4  degrees of freedom
    3 p/ Q9 z* h5 J. O
20. ## AIC: 34.0939 u+ `. u- B5 C$ f0 J
    
21. ##
    + U' y+ Y- K* c! T1 ^: Q
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率5 e3 |/ d& Y8 G6 e4 v7 s
   
2. 
   5 c. K  O4 q3 m. A3 l: L
3. #4、预测
   6 ~  F8 Y: S* R
4. pre <- predict(glm.sol, data.frame(x=3.5))
   & `) z- ~! v: {% x6 W
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1 5 U4 B, F: y- u- w! V2 G2 r
   
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1
   3 Z$ s\" N4 ^0 @0 y
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x ' H! x& `4 x3 s6 K
   
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：
   * d; T) k' d& Y/ y) H0 _8 ~* i
2. d <- seq(0, 5, length=100)+ ^* y& M+ c% R3 v+ e; [\" e
   
3. pre <- predict(glm.sol, data.frame(x=d))
   * [, G7 @9 y, N+ {\" W
4. p <- exp(pre)/(1+exp(pre))2 `0 o8 d) ~. Z5 G. [$ d
   
5. norell$y <- norell$success/norell$n$ J# d# j- d\" P4 n7 Y1 P
   
6. plot(norell$x, norell$y)
   4 x, d. n( p# ^\" x+ |
7. lines(d, p)

复制代码

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

复制代码

3 I. m# O6 z; ?