数学建模社区-数学中国
标题:
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 v
R软件提供了拟合计算广义线性模型的函数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 M
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)有完全相同的结果,但效率却低得多。
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/ W
Logistic回归模型的公式为: fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中,link=logit可以不写,因为logit是二项分布族连接函数的缺省状态。
: d- C2 M" }1 M: r$ y; D( D
实例一、Norell实验,高压电线对牲畜的影响
#1、加载数据
+ i/ k" }! L C$ U, Y" N
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
norell$Ymat<-cbind(norell$success, norell$n-norell$success)
6 q8 X) q0 p+ q2 O7 R
3 ?. n$ E" B, R
#2、建模
4 i& H! X2 R( s. [6 g7 k4 H' T% r
glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)
1 D, ]" a `1 c
7 F! y* F9 n# ~9 m
#3、模型评估
; s0 d' ? W3 H+ j+ c/ ]+ l* }* }
summary(glm.sol)
复制代码
##
" V5 @0 v h1 {6 u) x3 h
## Call:
) A$ l! w: M1 s" f$ L) t
## glm(formula = Ymat ~ x, family = binomial, data = norell)
9 D$ i8 W1 D1 z/ C9 j* t' s
##
# A7 q& e( Z5 g7 d" R7 Z
## Deviance Residuals:
7 Y6 {# [( H9 t: ]5 z
## 1 2 3 4 5 6
9 E( x8 S, Y# E" V& z( Q
## -2.2507 0.3892 -0.1466 1.1080 0.3234 -1.6679
$ L" L4 C% r8 D- Y
##
# e( m. t4 W( \: y. z
## Coefficients:
2 A8 ]3 Q! z/ C5 @9 ]
## Estimate Std. Error z value Pr(>|z|)
# {; |8 g. j- P: U1 Z
## (Intercept) -3.3010 0.3238 -10.20 <2e-16 ***
; U8 j! G0 J1 t8 D n1 {
## x 1.2459 0.1119 11.13 <2e-16 ***
0 t% S# r' L) O |, z0 ~$ O) c
## ---
5 q7 ^, ^: F: A A
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
) p7 ~8 X' r' M7 Y& h8 i
##
h* b( K2 k+ j2 I, B
## (Dispersion parameter for binomial family taken to be 1)
0 U& M! S/ @5 C8 y! o4 e
##
( y0 O9 u8 f) u1 o
## Null deviance: 250.4866 on 5 degrees of freedom
; v! i4 k) z; }7 H4 k+ E- W+ W
## Residual deviance: 9.3526 on 4 degrees of freedom
+ X* ~: k- J2 y' |! ?; a
## AIC: 34.093
$ m- L# E; |+ S$ q9 z
##
: G& b( n; W- E' ?
## Number of Fisher Scoring iterations: 4
复制代码
#与线性回归模型相同,在得到回归模型后,可以作预测:电流强度为3.5毫安时,有响应的牛的概率
% f4 P- L9 x& J: @, T! @
6 G+ ?" R. ?" j( }% p5 G4 Q" ?! h
#4、预测
1 E* P0 O3 w2 ]& \
pre <- predict(glm.sol, data.frame(x=3.5))
) C1 }. T$ K. K" U6 L, M! A) b( r
(p <- exp(pre)/(1+exp(pre)))
复制代码
## 1
8 V/ I& d% D H1 f- P3 J& {4 H& O& e4 x
## 0.742642
复制代码
#求有50%的牛响应时的电流强度:当P=0.5时,ln(P/(1-P))=0,所以X=-b0/b1
' U, V& U6 X, x- z7 X
glm.sol$coefficients
复制代码
## (Intercept) x
! P# H! g5 K) ?' f: @
## -3.301035 1.245937
复制代码
(X <- -glm.sol$coefficients[[1]]/glm.sol$coefficients[[2]])
复制代码
## [1] 2.649439
复制代码
#5、画出响应比例与logistic回归曲线:
& `, }0 Z$ Q! b/ _+ m/ D
d <- seq(0, 5, length=100)
' `% M1 G! a" S8 a) L- \
pre <- predict(glm.sol, data.frame(x=d))
) z% b3 _* L. p; V3 ]
p <- exp(pre)/(1+exp(pre))
" P- d) d( |3 S- m9 F
norell$y <- norell$success/norell$n
( _8 L; c' p: v( D9 Q6 h P* c; B
plot(norell$x, norell$y)
, i) m% {, t/ o! h/ S0 U
lines(d, p)
复制代码
#其中,d是给出曲线横坐标的点,pre是计算预测值,p是相应的预测概率。用plot函数和lines给出散点图和对应的预测曲线。
复制代码
% D6 @/ {. O4 W9 e+ F
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5