数学建模社区-数学中国

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

---

作者: 2744557306    时间: 2023-11-30 17:34
标题: Logistic回归实例2
# logistic回归
实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布，那么得到的是线性最小二乘回归，当随机变量服从伯努利分布，则得到的是Logistic回归。

9 i; r1 t; g' |, t. o5 AR软件提供了拟合计算广义线性模型的函数glm()，其命令格式如下：fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中，formula是拟合公式；family是分布族，即前面讲到的广义线性模型的种类，如正态分布、Poisson分布、二项分布等。
8 g* b8 j3 G: F% k7 P/ n  W( s/ v5 C

+ X3 P2 [  N$ a' ^* h有了上面这些分布族和连接函数，我们就可以完成相应的广义线性模型的拟合问题。

8 c$ }2 C  C' V& A* R- }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)有完全相同的结果，但效率却低得多。

2）二项分布
# O2 [: f/ R# Q% s3 ^

0 }7 X+ n( u1 F: P* llogistic回归模型是一个非线性回归模型，自变量可以是连续变量，也可以是分类变量，或哑变量。但可以使用线性回归模型对参数进行估计，所以Logistic回归模型属于广义线性模型。

( M7 C! d+ f" @' _Logistic回归模型的公式为： fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中，link=logit可以不写，因为logit是二项分布族连接函数的缺省状态。
实例一、Norell实验，高压电线对牲畜的影响

1. #1、加载数据
   
2. norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )
   ) Q( r( h. U! h% _( _% H
3. norell$Ymat<-cbind(norell$success, norell$n-norell$success)  
   
4. 
   
5. #2、建模
   
6. glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)
   : Y# g' ^* T- ?3 Y& Q8 H' g8 v
7. 
   
8. #3、模型评估
   
9. summary(glm.sol)

复制代码

1. ##
   % ~$ P- s9 g. @  y1 l3 S* C; P
2. ## Call:
   5 y9 w# p" h7 B9 ~
3. ## glm(formula = Ymat ~ x, family = binomial, data = norell)
   # y1 R# J) d5 B
4. ##
   
5. ## Deviance Residuals:
   & e- M( D( ~2 i4 c; C. t
6. ##       1        2        3        4        5        6  
   1 ^# X. V2 w1 c$ s: j4 B( `4 ?; }0 [2 c
7. ## -2.2507   0.3892  -0.1466   1.1080   0.3234  -1.6679  
   7 U9 E7 u- u- c
8. ##
   ) R  n; b8 O5 S/ o4 Y
9. ## Coefficients:
   / A7 ^1 I2 X* O( w  `1 i
10. ##             Estimate Std. Error z value Pr(>|z|)   
    
11. ## (Intercept)  -3.3010     0.3238  -10.20   <2e-16 ***
    
12. ## x             1.2459     0.1119   11.13   <2e-16 ***
    
13. ## ---
    
14. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
15. ##
    
16. ## (Dispersion parameter for binomial family taken to be 1)
    
17. ##
    
18. ##     Null deviance: 250.4866  on 5  degrees of freedom
    
19. ## Residual deviance:   9.3526  on 4  degrees of freedom
    4 g3 B4 x3 L6 X0 |, Z
20. ## AIC: 34.093
    " \* k# }: V: d# k8 k2 M, i
21. ##
    # ^. T( e! |8 _" c3 h6 k
22. ## Number of Fisher Scoring iterations: 4

复制代码

1. #与线性回归模型相同，在得到回归模型后，可以作预测：电流强度为3.5毫安时，有响应的牛的概率
   7 K5 x9 k! X5 w2 q0 V( F4 P; W
2. 
   5 \9 j3 b% V9 W5 g* t5 ^
3. #4、预测
   7 U6 h' J$ j- G3 d/ A4 {" z
4. pre <- predict(glm.sol, data.frame(x=3.5))
   
5. (p <- exp(pre)/(1+exp(pre)))

复制代码

1. ##        1
   $ L# _9 b/ S, G" J6 M1 f
2. ## 0.742642

复制代码

1. #求有50%的牛响应时的电流强度：当P=0.5时，ln(P/(1-P))=0,所以X=-b0/b1
   5 b: \8 H, b7 S$ H9 v2 W' @
2. glm.sol$coefficients

复制代码

1. ## (Intercept)           x
   
2. ##   -3.301035    1.245937

复制代码

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

复制代码

1. ## [1] 2.649439

复制代码

1. #5、画出响应比例与logistic回归曲线：
   ) x  G/ u4 m9 F$ J/ A4 u" }) ]
2. d <- seq(0, 5, length=100)
   " r: r7 y) W9 Y$ G+ z! l; I
3. pre <- predict(glm.sol, data.frame(x=d))
   
4. p <- exp(pre)/(1+exp(pre))
   
5. norell$y <- norell$success/norell$n
   + h! A1 t2 L8 C4 r: L+ f7 F
6. plot(norell$x, norell$y)
   % X. x: c: y# g
7. lines(d, p)

复制代码

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

复制代码

---

欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)

Powered by Discuz! X2.5