数学建模社区-数学中国

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

作者: 2744557306    时间: 2023-11-30 17:34
标题: Logistic回归实例2
# logistic回归
2 r& T# x2 S0 W7 p! N7 @实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布,那么得到的是线性最小二乘回归,当随机变量服从伯努利分布,则得到的是Logistic回归。
3 q9 N8 k+ i$ T
& U( ^+ ?0 h# @* U" t5 vR软件提供了拟合计算广义线性模型的函数glm(),其命令格式如下:fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中,formula是拟合公式;family是分布族,即前面讲到的广义线性模型的种类,如正态分布、Poisson分布、二项分布等。6 t5 u( ~# r- I) j/ N9 T

" T7 b1 M. o, i. U/ H1 s
* n( W, m3 ]3 M1 l有了上面这些分布族和连接函数,我们就可以完成相应的广义线性模型的拟合问题。
# w# F& B9 n- E+ R
3 f; k! N/ l4 b5 r" k! c' t' w5 M1)正态分布 正态分布族的使用方法: 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)有完全相同的结果,但效率却低得多。1 p2 R5 O1 z; X3 T- j
% s# m0 n1 J# ~. n
2)二项分布' L% V) @% |/ t3 ~+ M4 g

; {( j8 e9 j5 r3 v) }2 c* v: j, p( f% R9 e" b  T
logistic回归模型是一个非线性回归模型,自变量可以是连续变量,也可以是分类变量,或哑变量。但可以使用线性回归模型对参数进行估计,所以Logistic回归模型属于广义线性模型。
0 E# {( c6 K' ~) n1 A! c4 m
% V$ |& z+ t* j, I( y/ WLogistic回归模型的公式为: fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中,link=logit可以不写,因为logit是二项分布族连接函数的缺省状态。: d- C2 M" }1 M: r$ y; D( D
实例一、Norell实验,高压电线对牲畜的影响
  1. #1、加载数据
    + i/ k" }! L  C$ U, Y" N
  2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )- T4 c: w" N9 x; l7 _4 V
  3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  6 q8 X) q0 p+ q2 O7 R

  4. 3 ?. n$ E" B, R
  5. #2、建模4 i& H! X2 R( s. [6 g7 k4 H' T% r
  6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)1 D, ]" a  `1 c
  7. 7 F! y* F9 n# ~9 m
  8. #3、模型评估; s0 d' ?  W3 H+ j+ c/ ]+ l* }* }
  9. summary(glm.sol)
复制代码
  1. ##
    " V5 @0 v  h1 {6 u) x3 h
  2. ## Call:) A$ l! w: M1 s" f$ L) t
  3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)9 D$ i8 W1 D1 z/ C9 j* t' s
  4. ##
    # A7 q& e( Z5 g7 d" R7 Z
  5. ## Deviance Residuals: 7 Y6 {# [( H9 t: ]5 z
  6. ##       1        2        3        4        5        6  9 E( x8 S, Y# E" V& z( Q
  7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
    $ L" L4 C% r8 D- Y
  8. ##
    # e( m. t4 W( \: y. z
  9. ## Coefficients:2 A8 ]3 Q! z/ C5 @9 ]
  10. ##             Estimate Std. Error z value Pr(>|z|)   
    # {; |8 g. j- P: U1 Z
  11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***; U8 j! G0 J1 t8 D  n1 {
  12. ## x             1.2459     0.1119   11.13   <2e-16 ***0 t% S# r' L) O  |, z0 ~$ O) c
  13. ## ---5 q7 ^, ^: F: A  A
  14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ) p7 ~8 X' r' M7 Y& h8 i
  15. ##   h* b( K2 k+ j2 I, B
  16. ## (Dispersion parameter for binomial family taken to be 1)0 U& M! S/ @5 C8 y! o4 e
  17. ## ( y0 O9 u8 f) u1 o
  18. ##     Null deviance: 250.4866  on 5  degrees of freedom
    ; v! i4 k) z; }7 H4 k+ E- W+ W
  19. ## Residual deviance:   9.3526  on 4  degrees of freedom
    + X* ~: k- J2 y' |! ?; a
  20. ## AIC: 34.093$ m- L# E; |+ S$ q9 z
  21. ##
    : G& b( n; W- E' ?
  22. ## Number of Fisher Scoring iterations: 4
复制代码
  1. #与线性回归模型相同,在得到回归模型后,可以作预测:电流强度为3.5毫安时,有响应的牛的概率% f4 P- L9 x& J: @, T! @
  2. 6 G+ ?" R. ?" j( }% p5 G4 Q" ?! h
  3. #4、预测
    1 E* P0 O3 w2 ]& \
  4. pre <- predict(glm.sol, data.frame(x=3.5))
    ) C1 }. T$ K. K" U6 L, M! A) b( r
  5. (p <- exp(pre)/(1+exp(pre)))
复制代码
  1. ##        1
    8 V/ I& d% D  H1 f- P3 J& {4 H& O& e4 x
  2. ## 0.742642
复制代码
  1. #求有50%的牛响应时的电流强度:当P=0.5时,ln(P/(1-P))=0,所以X=-b0/b1' U, V& U6 X, x- z7 X
  2. glm.sol$coefficients
复制代码
  1. ## (Intercept)           x ! P# H! g5 K) ?' f: @
  2. ##   -3.301035    1.245937
复制代码
  1. (X <- -glm.sol$coefficients[[1]]/glm.sol$coefficients[[2]])
复制代码
  1. ## [1] 2.649439
复制代码
  1. #5、画出响应比例与logistic回归曲线:& `, }0 Z$ Q! b/ _+ m/ D
  2. d <- seq(0, 5, length=100)' `% M1 G! a" S8 a) L- \
  3. pre <- predict(glm.sol, data.frame(x=d))) z% b3 _* L. p; V3 ]
  4. p <- exp(pre)/(1+exp(pre))
    " P- d) d( |3 S- m9 F
  5. norell$y <- norell$success/norell$n( _8 L; c' p: v( D9 Q6 h  P* c; B
  6. plot(norell$x, norell$y), i) m% {, t/ o! h/ S0 U
  7. lines(d, p)
复制代码
  1. #其中,d是给出曲线横坐标的点,pre是计算预测值,p是相应的预测概率。用plot函数和lines给出散点图和对应的预测曲线。
复制代码

% D6 @/ {. O4 W9 e+ F




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