1 / 52

非平衡数据非线性混合模型

非平衡数据非线性混合模型. 数据. 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 ))

pabla
Download Presentation

非平衡数据非线性混合模型

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. 非平衡数据非线性混合模型

  2. 数据 • iq {spida} • 创伤性脑损伤后智商恢复数据 • （IQ after traumatic brain injuries）

  3. 数据描述 • DAYSPC: 昏迷苏醒后天数 • DCOMA：昏迷天数 • SEX：性别 • AgeIn： 受伤时的年龄 • ID：身份 • PIQ：行为IQ • VIQ：语言IQ

  4. 数据集 • iq为原始数据集 • iqsim为针对一个病人在10个观测恢复天数上的模拟数据集

  5. 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)

  6. 非线性函数 • 多项式(polynomials) • 样条（splines） • 指数增长（衰减）exponential growth(decay) • 周期（periodic）

  7. 多项式方程（Polynomials） • 线性方程（linear） • 二次方程（quadratic） • 三次方程(cubic) • 四次方程(quartic)

  8. 拟合程序 • plot(iq~days,iqsim,xlim=c(0,800),ylim=c(85,120)) • 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)

  9. 二次拟合 • 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)

  10. 三次拟合 • 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)

  11. 高次拟合 • 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)

  12. 指数增长函数

  13. 指数衰减函数

  14. 指数渐进增长函数 • IRRR（瞬间相对恢复率）

  15. 非线性增长曲线模型 • 公式 iq~b0+b1*exp(-alpha*days) • b0长期水平，b1在时间0点的相对差额 alpha恢复率

  16. 程序 • 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)

  17. 转换时间 • ttime<-function(x)exp(-0.0058*x) • fit.lin<-lm(iq~ttime(days),iqsim) • summary(fit.lin)

  18. 结果比较

  19. 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 )

  20. 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 )

  21. NLME • Their algorithm proceeds in two alternating steps, a penalized nonlinear least squares (PNLS)step and a linear mixed- effects (LME) step

  22. 混合模型

  23. D=G 随机效应的协方差矩阵

  24. 最大循环100次，不收敛也返回值 • verbose=T 显示PNLS和 LME步骤

  25. plot(fit.nlme,resid(.,type="p")~fitted(.)) • plot(fit.nlme,sqrt(abs(resid(.,type="p")))~fitted(.)) • plot(ranef(fit.nlme)) • pairs(ranef(fit.nlme))

  26. 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 )

  27. 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,-1,-10,-5,0.01)), • verbose=T • )

  28. 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 )

  29. 给定不同的初始值结果并不相同。 • start=list(fixed=c(100,-0.3,-10,0,0.3)) • b1的随机效应可去掉

  30. 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)

  31. 比较

  32. 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))

  33. 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))

  34. 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))

  35. 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 ))

  36. 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)

More Related