数学建模社区-数学中国

标题: Logistic回归实例2 [打印本页]

作者: 2744557306    时间: 2023-11-30 17:34
标题: Logistic回归实例2
# logistic回归: _- Z% N7 i  x2 y
实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布,那么得到的是线性最小二乘回归,当随机变量服从伯努利分布,则得到的是Logistic回归。- z; G' b# i5 S

3 N* a$ J, D- S( N/ A1 IR软件提供了拟合计算广义线性模型的函数glm(),其命令格式如下:fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中,formula是拟合公式;family是分布族,即前面讲到的广义线性模型的种类,如正态分布、Poisson分布、二项分布等。* k% z, B, E0 q+ D  L. `
4 m5 w3 O1 V. q8 S0 i

$ B, s4 t' U) ?) v0 t$ K' G有了上面这些分布族和连接函数,我们就可以完成相应的广义线性模型的拟合问题。3 |! T2 n7 y, C) I3 ?* i) [+ [, _- N

( Y  X4 m% o4 {) ~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)有完全相同的结果,但效率却低得多。/ c+ a3 @$ d2 a) |( J
# O0 s8 _6 r- k3 B  l  t
2)二项分布! K+ ~3 P; A; @+ @; J

8 z7 Y1 W% o* G3 W0 C! X( ~6 ]
7 v6 i5 m7 I/ a1 q4 \logistic回归模型是一个非线性回归模型,自变量可以是连续变量,也可以是分类变量,或哑变量。但可以使用线性回归模型对参数进行估计,所以Logistic回归模型属于广义线性模型。' t7 _- i2 e* h3 r% I! @

: A* K" T# |) y* R9 Q2 z% q# I' JLogistic回归模型的公式为: fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中,link=logit可以不写,因为logit是二项分布族连接函数的缺省状态。! Z  U& b( r! N( C* w5 Y1 T
实例一、Norell实验,高压电线对牲畜的影响
  1. #1、加载数据
    # t0 S8 V* [/ v) C- p4 a/ N9 E
  2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )
    - U2 p/ {& U  ~  V7 k
  3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  # X% V& h2 u: ~# b# m

  4. ( U1 u( I6 R" Z
  5. #2、建模
    ' I2 ?) P9 D: p" f
  6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)
      E) {2 S  R: F7 g. |9 p4 g
  7.   B7 {+ j6 s; T# Z
  8. #3、模型评估
    . T1 p  l' o$ Y: o! n7 z
  9. summary(glm.sol)
复制代码
  1. ##
    9 e2 Q" h4 w/ w6 n, p! u) F
  2. ## Call:
    : }  Z& X5 Z- `/ I: x2 m# n
  3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)
    % d/ _3 Y. W! W+ _. i1 Z, C: Q
  4. ## 6 H* E: i" C8 n* r: w& [3 W& u
  5. ## Deviance Residuals: $ R; C7 t/ c/ T; E$ ^( `
  6. ##       1        2        3        4        5        6  
    ; x! f/ L7 I2 k% s+ E0 x5 j0 G
  7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  2 w* U' E" @# `# p
  8. ##
    - N& k9 u+ {; j6 H. q; G
  9. ## Coefficients:7 v  J2 w  N; d& |
  10. ##             Estimate Std. Error z value Pr(>|z|)   
    . v' j* `1 i% ]5 Z" D8 @& g; L
  11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***9 t! j  G7 w2 B! G) e& ]. H
  12. ## x             1.2459     0.1119   11.13   <2e-16 ***5 E9 H& {* _! o8 R9 u' z$ N* K
  13. ## ---3 m5 `1 U3 T. U8 G, k+ }' A
  14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12 G7 O1 S, O( {( \. x3 }$ i
  15. ##
    , s- {- o7 p  |! V1 K1 C
  16. ## (Dispersion parameter for binomial family taken to be 1); f$ P) M, [2 k6 h4 Z! |
  17. ##
    0 T9 n$ x& X" i5 h$ b3 j
  18. ##     Null deviance: 250.4866  on 5  degrees of freedom# f1 {' i: |  [9 {" v! q. k* I
  19. ## Residual deviance:   9.3526  on 4  degrees of freedom) F1 U/ ~5 I- Q' ]+ r, L" j" ^
  20. ## AIC: 34.0934 Y* ?& n' _$ |
  21. ## / t' {3 [6 |! b' ~3 N3 g
  22. ## Number of Fisher Scoring iterations: 4
复制代码
  1. #与线性回归模型相同,在得到回归模型后,可以作预测:电流强度为3.5毫安时,有响应的牛的概率
    6 s6 d8 Q& Q+ ~$ l! C- P
  2. 2 @8 B' q8 ]* d$ n! g! o; v, F
  3. #4、预测. ~  i, B' C$ I
  4. pre <- predict(glm.sol, data.frame(x=3.5))3 R8 u4 n2 H) o& a$ i
  5. (p <- exp(pre)/(1+exp(pre)))
复制代码
  1. ##        1
    9 t! G8 |4 K6 X' F! g. D
  2. ## 0.742642
复制代码
  1. #求有50%的牛响应时的电流强度:当P=0.5时,ln(P/(1-P))=0,所以X=-b0/b1+ _/ @% s, P. W. R( \
  2. glm.sol$coefficients
复制代码
  1. ## (Intercept)           x
    ' @* g5 }- x1 T1 f
  2. ##   -3.301035    1.245937
复制代码
  1. (X <- -glm.sol$coefficients[[1]]/glm.sol$coefficients[[2]])
复制代码
  1. ## [1] 2.649439
复制代码
  1. #5、画出响应比例与logistic回归曲线:3 M" Y* }# g1 c# ]4 y8 l% `0 y5 G
  2. d <- seq(0, 5, length=100)' O. Z& j  d+ {1 v( v7 S* H
  3. pre <- predict(glm.sol, data.frame(x=d))
    4 I6 \* h* B; i0 Q# J6 A
  4. p <- exp(pre)/(1+exp(pre))5 E3 X$ X5 [5 b- h
  5. norell$y <- norell$success/norell$n
    $ d) y( B6 ^9 l& D4 ?# _) i8 J5 |
  6. plot(norell$x, norell$y)7 s: j# h- n9 M- p) U
  7. lines(d, p)
复制代码
  1. #其中,d是给出曲线横坐标的点,pre是计算预测值,p是相应的预测概率。用plot函数和lines给出散点图和对应的预测曲线。
复制代码

2 j  k( S& N, D0 Z. O8 ?. J




欢迎光临 数学建模社区-数学中国 (http://www.madio.net/) Powered by Discuz! X2.5