数学建模社区-数学中国
标题:
Logistic回归实例2
[打印本页]
作者:
2744557306
时间:
2023-11-30 17:34
标题:
Logistic回归实例2
# logistic回归
" O2 K3 I* s4 V t
实际上线性最小二乘回归和Logistic回归都是广义线性模型的一个特例。当随机变量Y服从高斯分布,那么得到的是线性最小二乘回归,当随机变量服从伯努利分布,则得到的是Logistic回归。
: n r9 V# X" N" T
% ^ A. G6 n6 Q; \2 i: R4 d5 o: I
R软件提供了拟合计算广义线性模型的函数glm(),其命令格式如下:fitted.model <- glm(formula, family=family.generator, data=data.frame) 其中,formula是拟合公式;family是分布族,即前面讲到的广义线性模型的种类,如正态分布、Poisson分布、二项分布等。
$ `1 p, T& u" r4 u* a
* F& s% r; M$ u' w4 A5 Z; o
5 D+ {; p5 ]4 p, `. d' Z
有了上面这些分布族和连接函数,我们就可以完成相应的广义线性模型的拟合问题。
, B$ `; [- ]: Z
4 \) D* [- y0 |
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)有完全相同的结果,但效率却低得多。
& n7 `( J2 H! K: ?
) E5 X% {) k1 Z8 c' R5 x5 n
2)二项分布
$ f5 h' B5 T- `4 @8 a/ t: \
# I! X+ o' ?1 X# Y
% Z e: ~, @" s# |5 I
logistic回归模型是一个非线性回归模型,自变量可以是连续变量,也可以是分类变量,或哑变量。但可以使用线性回归模型对参数进行估计,所以Logistic回归模型属于广义线性模型。
/ K' d d, E) Z# Y
) K! m) ~ E3 w7 [
Logistic回归模型的公式为: fm <- glm(formula, family=binomial(link=logit), data=data.frame) 其中,link=logit可以不写,因为logit是二项分布族连接函数的缺省状态。
8 ]- Y- {* U. p5 U8 j
实例一、Norell实验,高压电线对牲畜的影响
#1、加载数据
! A. c" @) m# u6 n4 q
norell<-data.frame( x=0:5, n=rep(70,6), success=c(0,9,21,47,60,63) )
( d1 y- l+ n/ H3 U5 C
norell$Ymat<-cbind(norell$success, norell$n-norell$success)
3 ^0 y/ u. F( l
! ]6 ]6 I; Y4 j! [
#2、建模
3 a- t- a6 T* S
glm.sol <- glm(Ymat ~ x, family=binomial, data=norell)
3 g$ ~9 W" ~6 I5 |8 \
* Q% Q$ `* N* G$ V+ D, B
#3、模型评估
5 P ? q7 {4 l, g& F4 A& L
summary(glm.sol)
复制代码
##
# h ~, }' t- z" k
## Call:
( @ i* `0 g- V e- S: t. E
## glm(formula = Ymat ~ x, family = binomial, data = norell)
; v% M1 e) u; f2 R0 a; e$ ? ~
##
8 W' Y) ]- o, U3 F2 [
## Deviance Residuals:
) i6 Q8 P" |4 q2 {
## 1 2 3 4 5 6
7 v, V) G6 O. e+ O u
## -2.2507 0.3892 -0.1466 1.1080 0.3234 -1.6679
+ S& `, |9 o; t5 A: B
##
1 h! W# B+ [/ F% v# V1 A% ~8 I
## Coefficients:
9 X c6 b u& a7 v* _ e* X
## Estimate Std. Error z value Pr(>|z|)
~6 {7 }5 c( ~3 ]( t
## (Intercept) -3.3010 0.3238 -10.20 <2e-16 ***
" E6 y' d/ C: {# n0 F4 e
## x 1.2459 0.1119 11.13 <2e-16 ***
6 ~. J z7 n. ?- j
## ---
" J2 ~) @" v$ { N
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
0 u7 |8 G+ {. K3 O+ `, S
##
8 s9 Y/ z/ R, v) z. O- X$ W$ ]% R
## (Dispersion parameter for binomial family taken to be 1)
7 ~6 I* I; _& t# m* P, _
##
' c n0 }8 N; r6 Q, _1 G
## Null deviance: 250.4866 on 5 degrees of freedom
5 ]: S. C2 w) F
## Residual deviance: 9.3526 on 4 degrees of freedom
0 e" K7 j3 d: B5 H w; G) s
## AIC: 34.093
( }/ N3 U! c. d. z; C
##
! q7 b) V7 F1 y& D
## Number of Fisher Scoring iterations: 4
复制代码
#与线性回归模型相同,在得到回归模型后,可以作预测:电流强度为3.5毫安时,有响应的牛的概率
1 \' R n2 p! J* z
. g* ^: Y! ?( X, c9 }
#4、预测
P m+ R9 Q' w, e2 u2 S" }' q( K) _
pre <- predict(glm.sol, data.frame(x=3.5))
7 q! y8 @( }# ` x1 O
(p <- exp(pre)/(1+exp(pre)))
复制代码
## 1
: W0 U" f, u! @( j
## 0.742642
复制代码
#求有50%的牛响应时的电流强度:当P=0.5时,ln(P/(1-P))=0,所以X=-b0/b1
1 Q8 T# d n( F& j0 m
glm.sol$coefficients
复制代码
## (Intercept) x
4 @, v. g( R( {$ y$ L& `5 }8 D. k
## -3.301035 1.245937
复制代码
(X <- -glm.sol$coefficients[[1]]/glm.sol$coefficients[[2]])
复制代码
## [1] 2.649439
复制代码
#5、画出响应比例与logistic回归曲线:
9 z( Z1 R8 G- k
d <- seq(0, 5, length=100)
7 Y1 }) L/ d8 w. X! M
pre <- predict(glm.sol, data.frame(x=d))
# e+ H* z V; ^) R) j- t
p <- exp(pre)/(1+exp(pre))
c" {6 s' k. f2 \
norell$y <- norell$success/norell$n
" r6 ?. M' s, ?1 z
plot(norell$x, norell$y)
6 I- s$ q* D. ~# J
lines(d, p)
复制代码
#其中,d是给出曲线横坐标的点,pre是计算预测值,p是相应的预测概率。用plot函数和lines给出散点图和对应的预测曲线。
复制代码
. C) h% R* G# t& l' m$ y0 v
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5