数学建模社区-数学中国
标题:
Logistic回归实例2
[打印本页]
作者:
2744557306
时间:
2023-11-30 17:34
标题:
Logistic回归实例2
# logistic回归
+ Z& L% E3 B6 I$ J2 \/ T% l/ ~
实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布,那么得到的是线性最小二乘回归,当随机变量服从伯努利分布,则得到的是Logistic回归。
3 L& f/ d2 `1 Z" P7 O
. `, g( d. o9 o1 g8 \; r2 Q
R软件提供了拟合计算广义线性模型的函数glm(),其命令格式如下:fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中,formula是拟合公式;family是分布族,即前面讲到的广义线性模型的种类,如正态分布、Poisson分布、二项分布等。
; t2 b) u9 ?3 D' X9 ^$ a
# }3 C: Y$ x8 Y' D
* P8 `/ \. {4 P; Q/ m7 E
有了上面这些分布族和连接函数,我们就可以完成相应的广义线性模型的拟合问题。
: i0 s2 J4 z( \* s/ n @
, j( U. I: G; [* f" W5 x9 C
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)有完全相同的结果,但效率却低得多。
# P$ u% T1 `$ U* }2 ~
# `9 Z& r' u# W$ W* b
2)二项分布
- {8 |$ t3 ^, ^
5 X% c4 t3 ~0 ]3 j' t# ]
( ^, v8 w- {% B. n
logistic回归模型是一个非线性回归模型,自变量可以是连续变量,也可以是分类变量,或哑变量。但可以使用线性回归模型对参数进行估计,所以Logistic回归模型属于广义线性模型。
3 s7 t2 y8 \8 y b5 u* g
2 P( R0 F5 E) O5 a
Logistic回归模型的公式为: fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中,link=logit可以不写,因为logit是二项分布族连接函数的缺省状态。
& G! a! x1 d7 ]
实例一、Norell实验,高压电线对牲畜的影响
#1、加载数据
a( S/ ?( W0 h( ^$ H- r7 Q
norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )
: `, ~; J5 O, x
norell$Ymat<-cbind(norell$success, norell$n-norell$success)
) {# p+ t! J$ N
5 e3 |! i# x8 l" j
#2、建模
0 q( p# A# a4 N# }
glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)
2 A; \; {- B* p" R/ e0 P z; y6 S7 L; J/ E
/ S* U% G _/ h' Q7 x9 w+ E
#3、模型评估
9 G/ K! ^9 v B7 I' P
summary(glm.sol)
复制代码
##
* k4 t' {, f0 Y. a) s
## Call:
3 K. X- z8 }6 N) K
## glm(formula = Ymat ~ x, family = binomial, data = norell)
: I; c/ L" `+ r, t% Q# S$ \$ R
##
1 f6 W1 ?2 }- T) h
## Deviance Residuals:
' ^" w. w( p0 X& K3 Z
## 1 2 3 4 5 6
1 d/ |( e8 j$ N2 E5 J' {! z1 m: n
## -2.2507 0.3892 -0.1466 1.1080 0.3234 -1.6679
2 b" {) q+ A5 D l! W0 T% u# w
##
: I4 w3 T4 ]) @; Z4 E
## Coefficients:
) ]+ b6 } `5 ?8 K4 a+ ^
## Estimate Std. Error z value Pr(>|z|)
: m6 x: X; m8 A( J8 P
## (Intercept) -3.3010 0.3238 -10.20 <2e-16 ***
) i' V4 x6 c4 i4 |8 z
## x 1.2459 0.1119 11.13 <2e-16 ***
3 P) f& r; O: \3 T4 Q6 s$ w8 Y; Y
## ---
5 B! M1 i; b5 m0 z
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
: |% P) L/ w* @8 H4 I( |. d
##
: {. ?- B7 J& u4 R6 Z$ ^( l" q6 W
## (Dispersion parameter for binomial family taken to be 1)
" q8 M7 A$ l+ z" ~+ z
##
% V+ X* D2 m2 k- `
## Null deviance: 250.4866 on 5 degrees of freedom
. f7 Y6 B6 w/ D% y' I& Q+ ]2 p
## Residual deviance: 9.3526 on 4 degrees of freedom
7 [+ z2 A8 e. w! j) ]5 \
## AIC: 34.093
- V4 d; G0 R( c3 f( M e
##
* l+ J. c5 I( p$ o
## Number of Fisher Scoring iterations: 4
复制代码
#与线性回归模型相同,在得到回归模型后,可以作预测:电流强度为3.5毫安时,有响应的牛的概率
6 Q9 O f' h8 Z' D2 s
7 |* p2 N8 }/ P$ F
#4、预测
5 f% Y5 ^( R4 ?4 ?. j, ]0 U2 ]
pre <- predict(glm.sol, data.frame(x=3.5))
# `/ b+ S: Y1 U9 z% R2 h) F
(p <- exp(pre)/(1+exp(pre)))
复制代码
## 1
# Q1 K: ^# y# G$ j' @% C) R. u N
## 0.742642
复制代码
#求有50%的牛响应时的电流强度:当P=0.5时,ln(P/(1-P))=0,所以X=-b0/b1
: }. g. t. f/ b2 U0 U, {
glm.sol$coefficients
复制代码
## (Intercept) x
3 d% H2 r5 S# C) y* [& y) _
## -3.301035 1.245937
复制代码
(X <- -glm.sol$coefficients[[1]]/glm.sol$coefficients[[2]])
复制代码
## [1] 2.649439
复制代码
#5、画出响应比例与logistic回归曲线:
/ t# O% v7 M v1 D/ O6 k
d <- seq(0, 5, length=100)
" Z* k+ \" `% y" k
pre <- predict(glm.sol, data.frame(x=d))
. n0 l r6 w5 e; W
p <- exp(pre)/(1+exp(pre))
/ G0 |0 P u! K0 J/ o1 ^
norell$y <- norell$success/norell$n
/ Q% Y4 s4 c7 p
plot(norell$x, norell$y)
: Z3 \2 ^) O* u- \) p' N6 U
lines(d, p)
复制代码
#其中,d是给出曲线横坐标的点,pre是计算预测值,p是相应的预测概率。用plot函数和lines给出散点图和对应的预测曲线。
复制代码
$ `, f$ \) P/ |# v9 ^
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5