Transcript 非线性混合回归模型
非平衡数据非线性混合模型
数据
iq {spida}
创伤性脑损伤后智商恢复数据
(IQ after traumatic brain injuries)
数据描述
DAYSPC: 昏迷苏醒后天数
DCOMA:昏迷天数
SEX:性别
AgeIn: 受伤时的年龄
ID:身份
PIQ:行为IQ
VIQ:语言IQ
数据集
iq 为原始数据集
iqsim为针对一个病人在10个观测恢复天数
上的模拟数据集
tab(~ID,iq)
tab(~DAYSPC,iq)
tab(tab(~ID,iq))
tab(tab(~DAYSPC,iq))
xyplot( PIQ ~ DAYSPC , iq, groups = ID)
tab(~ID+SEX,iq)
iq.id=up(iq,~ID,all=T)
tab(~ID+SEX,iq.id)
非线性函数
多项式(polynomials)
样条(splines)
指数增长(衰减)exponential growth(decay)
周期(periodic)
多项式方程(Polynomials)
线性方程(linear)
二次方程(quadratic)
三次方程(cubic)
四次方程(quartic)
拟合程序
plot(iq~days,iqsim,xlim=c(0,800),ylim=c(85,12
0))
fit.lin<-lm(iq~days,iqsim)
pred<-expand.grid(days=seq(-20,800,1))
pred$iq.lin<-predict(fit.lin,pred)
lines(iq.lin~days, pred,lwd=2)
二次拟合
fit.quad<-lm(iq~days+I(days^2),iqsim)
pred$iq.quad<-predict(fit.quad,pred)
lines(iq.quad~days,pred,col="red",lwd=2)
三次拟合
fit.cub<-lm(iq~days+I(days^2)+I(days^3),iqsim)
summary(fit.cub)
pred$iq.cub<-predict(fit.cub,pred)
lines(iq.cub~days,pred,col="blue",lwd=2)
高次拟合
p8<-function(x) poly(x,8,raw=T)
#raw if true, use raw and not orthogonal
polynomials.
fit.8<-lm(iq~p8(days),iqsim)
summary(fit.8)
p9<-function(x) poly(x,9,raw=T)
fit.9<-lm(iq~p9(days),iqsim)
summary(fit.9)
pred$iq.8<-predict(fit.8,pred)
lines(iq.8~days,pred,col="green",lwd=2)
指数增长函数
指数衰减函数
指数渐进增长函数
IRRR(瞬间相对恢复率)
非线性增长曲线模型
公式 iq~b0+b1*exp(-alpha*days)
b0长期水平,b1在时间0点的相对差额
alpha恢复率
程序
fit.nl<-nls(iq~b0+b1*exp(alpha*days),iqsim,start=list(b0=100,b1=20,alpha=0.005))
summary(fit.nl)
pred$iq.nl<-predict(fit.nl,pred)
lines(iq.nl~days,pred,col="purple",lwd=2)
abline(h=coef(fit.nl)[1],col="gray",lwd=2)
转换时间
ttime<-function(x)exp(-0.0058*x)
fit.lin<-lm(iq~ttime(days),iqsim)
summary(fit.lin)
结果比较
IQ 数据
fit.nlme<-nlme(PIQ~b0+b1*exp(a*DAYSPC)+b.sex*(SEX=="Female"),
data=iq,
fixed=list(b0~1+sqrt(DCOMA),
b1~1+sqrt(DCOMA),a~1,b.sex~1),
random=list(ID=list(b0~1,b1~1)),
control=list(maxIter=200,returnObject=T),
start=list(fixed=c(100,-10,-20,1,0.05,0)),
verbose=T
)
IQ 数据
fit.nlme<-nlme(PIQ~b0+b1*exp(-a*DAYSPC),
data=iq,
fixed=list(b0~1+sqrt(DCOMA),
b1~1+sqrt(DCOMA),a~1),
random=list(ID=list(b0~1,b1~1)),
control=list(maxIter=200,returnObject=T),
start=list(fixed=c(100,-10,-20,1,0.05)),
verbose=T
)
NLME
Their algorithm proceeds in two alternating
steps, a penalized nonlinear least squares
(PNLS)
step and a linear mixed- effects (LME) step
混合模型
D=G 随机效应的协方差矩阵
最大循环100次,不收敛也返回值
verbose=T
显示PNLS和 LME步骤
plot(fit.nlme,resid(.,type="p")~fitted(.))
plot(fit.nlme,sqrt(abs(resid(.,type="p")))~fitted
(.))
plot(ranef(fit.nlme))
pairs(ranef(fit.nlme))
VIQ
fit.nlme.viq<-nlme(VIQ~b0+b1*exp(-a*DAYSPC),
data=iq,
fixed=list(b0~1+sqrt(DCOMA),
b1~1+sqrt(DCOMA),a~1),
random=list(ID=list(b0~1,b1~1)),
control=list(maxIter=100,returnObject=T),
start=list(fixed=c(100,-1,-10,-5,0.01)),
verbose=T
)
fit.nlme.viq<-nlme(VIQ~b0+b1*exp(-a*(DAYSPC30)),
data=iq,
fixed=list(b0~1+sqrt(DCOMA),
b1~1+sqrt(DCOMA),a~1),
random=list(ID=list(b0~1,b1~1)),
control=list(maxIter=100,returnObject=T),
start=list(fixed=c(100,-1,-10,-5,0.01)),
verbose=T
)
fit.nlme.viq<-nlme(VIQ~b0+b1*exp(-a*(DAYSPC-30)),
data=iq,
fixed=list(b0~1+sqrt(DCOMA),
b1~1+sqrt(DCOMA),a~1),
random=list(ID=list(b0~1,b1~1)),
control=list(maxIter=100,returnObject=T),
start=list(fixed=c(100,0,-10,-2,0.02)),
verbose=T,
subset=DCOMA<100
)
给定不同的初始值结果并不相同。
start=list(fixed=c(100,-0.3,-10,0,0.3))
b1的随机效应可去掉
fit.nlme.piq<-nlme(PIQ~b0+b1*exp(-a*(DAYSPC-30)),
data=iq,
fixed=list(b0~1+sqrt(DCOMA),
b1~1+sqrt(DCOMA),a~1),
random=list(ID=list(b0~1,b1~1)),
control=list(maxIter=100,returnObject=T),
start=list(fixed=c(100,-0.3,-10,0,0.1)),
verbose=T,
subset=DCOMA<100
)
summary(fit.nlme.piq)
比较
windows( height = 7, width = 8.5)
pred <- expand.grid( DCOMA=c(0,1,7,16,25,100),
DAYSPC = seq(30,365*2,5))
pred$piq <- predict( fit.nlme.piq, pred,level=0)
pred$viq <- predict( fit.nlme.viq, pred,level=0 )
zz <- factor( paste( 'DCOMA =', pred$DCOMA))
pred$DCOMA.lab <- reorder( zz, pred$DCOMA)
td( col = c('green','red'), lwd = 2)
xyplot( viq + piq ~ DAYSPC |DCOMA.lab, pred, type = 'l',
ylim = c(60,102),
lwd = 2, auto.key = list(columns = 2, points = F, lines = T))
pred2 <- expand.grid( DCOMA = 0:100,
DAYSPC = c(30,60,90.180,360,720))
pred2$piq <- predict( fit.nlme.piq, pred2, level = 0 )
pred2$viq <- predict( fit.nlme.viq, pred2, level = 0 )
zz <- factor( paste( 'DAYSPC =', pred2$DAYSPC))
pred2$dayspc.lab <- reorder( zz, pred2$DAYSPC)
xyplot( viq + piq ~ DCOMA | dayspc.lab, pred2, type = 'l',
ylim = c(60,102),
lwd = 2,
auto.key = list(columns = 2, points = F, lines = T))
predviq <- expand.grid( DCOMA = seq(0,100,10),
DAYSPC = seq(30,720,30))
predpiq <- predviq
predviq $ iq <- predict( fit.nlme.viq, predviq, level = 0)
predpiq $ iq <- predict( fit.nlme.piq , predpiq, level = 0)
predpiq$type <- factor( "PIQ" )
predviq$type <- factor( "VIQ" )
wireframe( iq ~ DAYSPC + DCOMA | type,
Rbind( predpiq, predviq))
wireframe( iq ~ DCOMA+DAYSPC | type,
Rbind( predpiq, predviq),scales = list( arrows =
F), col = 'blue')
wireframe( iq ~ DCOMA+DAYSPC | type,
Rbind(predpiq, predviq),
scales = list( arrows = F), col = 'blue',
xlab = 'Coma',
ylab = 'Days post coma',
screen = list( z = -65, x = -75 ))
wireframe( iq ~ DCOMA+DAYSPC , Rbind(predpiq,
predviq),
groups = type,
scales = list( arrows = F), col = 'blue',
xlab = 'Coma',
ylab = 'Days post coma',
screen = list( z = -65, x = -75 ),
auto.key = list(columns=2, lines = T, points = F),
alpha = .5)