Logistic回归实例2

2744557306

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

/ `$ @% @4 B6 p+ ]$ }2 H) pR软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。
3 B# \  \; u# e3 Q) w3 Y
4 R" Q) q. D! Q- f8 y5 c, u/ I
8 \/ X; R% w% z$ {+ k有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。
3 f: K; H* F2 G5 f: R
8 I- O# r3 L% S# ~4 H1）正态分布 正态分布族的使用方法： 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)有完全相同的结果，但效率却低得多。
8 B8 I  ?0 q  n
5 X( {7 A7 A+ P2）二项分布
8 \0 ]* F  l* e# g3 [( U

3 Q5 u& m' @% d; Slogistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。

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

1. #1、加载数据
   8 C7 W/ s: v9 g- ?( W- k/ d
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )
   % r4 }5 w* q: U$ V) J& G1 s% ~' \
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  
   9 l- }# T4 u$ _! a) d
4. ) P& i8 ~. ~$ }  O3 L8 p3 W
   
5. #2、建模
   2 O$ t' S8 {- }1 U' q4 p5 s
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)
   7 d7 P$ q: O' r1 E0 H
7. 
   / ?' z1 K  A6 d! {$ A5 N$ t, R( J' Z
8. #3、模型评估
     [/ q( {* \3 W; K; t0 b6 W3 B
9. summary(glm.sol)

复制代码

1. ##
   ! _  _/ q; ^: T# M2 k; [
2. ## Call:7 n9 b6 V9 o$ ~& D* q* u
   
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)
   + i* P5 A# X4 `& q6 O# T* ^( B
4. ## ! Z% N- ?( F0 \% a
   
5. ## Deviance Residuals:
   6 B3 ^3 J6 o5 t\" S/ s0 K
6. ##       1        2        3        4        5        6  . \% {* D/ t' U/ O% q/ s4 e% R* M
   
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
   , ?! _' Y5 |+ [+ b* a  b
8. ##
   / B9 E8 d3 h: @0 `
9. ## Coefficients:# u5 Y! s4 C- K! ]! r
   
10. ##             Estimate Std. Error z value Pr(>|z|)    \" A8 w% R0 k8 r0 a
    
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***% X0 E5 V7 G3 J0 n- t' E
    
12. ## x             1.2459     0.1119   11.13   <2e-16 ***
      G- J9 d/ d9 r. a
13. ## ---
    % R- w' K; C$ _9 E( p4 _
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 11 @1 y' a2 `% o
    
15. ##
    5 P6 T8 u9 `$ X5 }: _1 D6 c/ }% q
16. ## (Dispersion parameter for binomial family taken to be 1)! ^\" U3 G\" Z, @! t! K
    
17. ## - b% V+ ]% D4 k
    
18. ##     Null deviance: 250.4866  on 5  degrees of freedom
    8 _3 M* n2 P3 Q3 U
19. ## Residual deviance:   9.3526  on 4  degrees of freedom! G) h3 E9 E\" s  A, d
    
20. ## AIC: 34.093
      Y, n5 o) l4 {! i6 d2 l- d
21. ##
      t6 i- s/ T$ Q% `$ a
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率
   1 M+ R. o0 z. g  u
2. * u6 _& P7 f9 ~% m& A, C2 V1 s# k
   
3. #4、预测- F! T) H- I& t
   
4. pre <- predict(glm.sol, data.frame(x=3.5))
   8 I6 g0 H% J$ g
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1
   & `$ H! y4 K' j+ C& y& \1 a' X( _
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b11 F  k, u% Y5 t/ X7 k7 E7 N
   
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x
   : S1 ?; I! i2 y- ~\" P+ |. j5 J
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：
   \" d1 S$ f* r$ B
2. d <- seq(0, 5, length=100)- x- G9 X; i, ~6 G
   
3. pre <- predict(glm.sol, data.frame(x=d))
   ; M1 |; H. z* o
4. p <- exp(pre)/(1+exp(pre))' V9 R' L* k% r. z8 S5 W
   
5. norell$y <- norell$success/norell$n; @- ^5 C0 ^: `3 Z
   
6. plot(norell$x, norell$y)
   ' H3 V* h6 Z1 m. z
7. lines(d, p)

复制代码

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

复制代码

4 R, V) V- ]0 g' R3 O) v% g6 k  r