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)