数学建模社区-数学中国
标题:
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 I
R软件提供了拟合计算广义线性模型的函数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' J
Logistic回归模型的公式为: fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中,link=logit可以不写,因为logit是二项分布族连接函数的缺省状态。
! Z U& b( r! N( C* w5 Y1 T
实例一、Norell实验,高压电线对牲畜的影响
#1、加载数据
# t0 S8 V* [/ v) C- p4 a/ N9 E
norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )
- U2 p/ {& U ~ V7 k
norell$Ymat<-cbind(norell$success, norell$n-norell$success)
# X% V& h2 u: ~# b# m
( U1 u( I6 R" Z
#2、建模
' I2 ?) P9 D: p" f
glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)
E) {2 S R: F7 g. |9 p4 g
B7 {+ j6 s; T# Z
#3、模型评估
. T1 p l' o$ Y: o! n7 z
summary(glm.sol)
复制代码
##
9 e2 Q" h4 w/ w6 n, p! u) F
## Call:
: } Z& X5 Z- `/ I: x2 m# n
## glm(formula = Ymat ~ x, family = binomial, data = norell)
% d/ _3 Y. W! W+ _. i1 Z, C: Q
##
6 H* E: i" C8 n* r: w& [3 W& u
## Deviance Residuals:
$ R; C7 t/ c/ T; E$ ^( `
## 1 2 3 4 5 6
; x! f/ L7 I2 k% s+ E0 x5 j0 G
## -2.2507 0.3892 -0.1466 1.1080 0.3234 -1.6679
2 w* U' E" @# `# p
##
- N& k9 u+ {; j6 H. q; G
## Coefficients:
7 v J2 w N; d& |
## Estimate Std. Error z value Pr(>|z|)
. v' j* `1 i% ]5 Z" D8 @& g; L
## (Intercept) -3.3010 0.3238 -10.20 <2e-16 ***
9 t! j G7 w2 B! G) e& ]. H
## x 1.2459 0.1119 11.13 <2e-16 ***
5 E9 H& {* _! o8 R9 u' z$ N* K
## ---
3 m5 `1 U3 T. U8 G, k+ }' A
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
2 G7 O1 S, O( {( \. x3 }$ i
##
, s- {- o7 p |! V1 K1 C
## (Dispersion parameter for binomial family taken to be 1)
; f$ P) M, [2 k6 h4 Z! |
##
0 T9 n$ x& X" i5 h$ b3 j
## Null deviance: 250.4866 on 5 degrees of freedom
# f1 {' i: | [9 {" v! q. k* I
## Residual deviance: 9.3526 on 4 degrees of freedom
) F1 U/ ~5 I- Q' ]+ r, L" j" ^
## AIC: 34.093
4 Y* ?& n' _$ |
##
/ t' {3 [6 |! b' ~3 N3 g
## Number of Fisher Scoring iterations: 4
复制代码
#与线性回归模型相同,在得到回归模型后,可以作预测:电流强度为3.5毫安时,有响应的牛的概率
6 s6 d8 Q& Q+ ~$ l! C- P
2 @8 B' q8 ]* d$ n! g! o; v, F
#4、预测
. ~ i, B' C$ I
pre <- predict(glm.sol, data.frame(x=3.5))
3 R8 u4 n2 H) o& a$ i
(p <- exp(pre)/(1+exp(pre)))
复制代码
## 1
9 t! G8 |4 K6 X' F! g. D
## 0.742642
复制代码
#求有50%的牛响应时的电流强度:当P=0.5时,ln(P/(1-P))=0,所以X=-b0/b1
+ _/ @% s, P. W. R( \
glm.sol$coefficients
复制代码
## (Intercept) x
' @* g5 }- x1 T1 f
## -3.301035 1.245937
复制代码
(X <- -glm.sol$coefficients[[1]]/glm.sol$coefficients[[2]])
复制代码
## [1] 2.649439
复制代码
#5、画出响应比例与logistic回归曲线:
3 M" Y* }# g1 c# ]4 y8 l% `0 y5 G
d <- seq(0, 5, length=100)
' O. Z& j d+ {1 v( v7 S* H
pre <- predict(glm.sol, data.frame(x=d))
4 I6 \* h* B; i0 Q# J6 A
p <- exp(pre)/(1+exp(pre))
5 E3 X$ X5 [5 b- h
norell$y <- norell$success/norell$n
$ d) y( B6 ^9 l& D4 ?# _) i8 J5 |
plot(norell$x, norell$y)
7 s: j# h- n9 M- p) U
lines(d, p)
复制代码
#其中,d是给出曲线横坐标的点,pre是计算预测值,p是相应的预测概率。用plot函数和lines给出散点图和对应的预测曲线。
复制代码
2 j k( S& N, D0 Z. O8 ?. J
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5