# # cfkids-CDA-NewLMM.q # # ------------------------------------------------------------ # # PURPOSE: Use linear mixed models to characterize longitudinal # change by gender and genotype. # # AUTHOR: P. Heagerty # # DATE: 00/07/10 Revised 14Feb2002 # # ------------------------------------------------------------ # # ##### ##### Read data ##### # source("cfkids-read.q") # library(nlme) # ##### ##### Trellis plots of individuals and groups ##### # # Create Grouped Data Set # ntotal <- cumsum( unlist( lapply( split( cfkids$id, cfkids$id) , length ) ) ) cf.subset <- groupedData( fev1 ~ age | id, outer = ~ factor(f508)*female, data = cfkids[ 1:ntotal[(8*4*1)], ] ) # cfkids <- groupedData( fev1 ~ age | id, outer = ~ factor(f508)*female, data = cfkids ) # # trellis plot, by id, first 1 pages, 8x4 # postscript( file="cfkids-trellis.1.ps", horiz=F ) plot( cf.subset, layout = c(4,8) ) graphics.off() postscript( file="cfkids-trellis.2.ps", horiz=T ) par( pch="." ) plot( cfkids, outer = ~ factor(f508)*factor(female), layout=c(3,2), aspect=1 ) graphics.off() # ##### ##### Linear Mixed Models ##### # options( contrasts=c("contr.treatment","contr.helmert") ) # ### Intercept only fit0 <- lme( fev1 ~ age0 + ageL + female*ageL + factor(f508)*ageL, method = "ML", random = reStruct( ~ 1 | id, pdClass="pdSymm", REML=F), data = cfkids ) summary( fit0 ) ### Intercept plus Slope fit1 <- lme( fev1 ~ age0 + ageL + female*ageL + factor(f508)*ageL, method = "ML", random = reStruct( ~ 1 + ageL | id, pdClass="pdSymm", REML=F), data = cfkids ) summary( fit1 ) ### EDA for serial correlation postscript( file="cfkids-Variogram.ps", horiz=T ) plot( Variogram( fit0, form = ~ age | id , resType="response" ) ) graphics.off() ### Intercept plus AR(1) fit2a <- lme( fev1 ~ age0 + ageL + female*ageL + factor(f508)*ageL, method = "ML", random = reStruct( ~ 1 | id, pdClass="pdSymm", REML=F), correlation = corAR1( form = ~ 1 | id ), data = cfkids ) summary( fit2a ) ### another way fit2b <- lme( fev1 ~ age0 + ageL + female*ageL + factor(f508)*ageL, method = "ML", random = reStruct( ~ 1 | id, pdClass="pdSymm", REML=F), correlation = corExp( form = ~ ageL | id, nugget=F), data = cfkids ) summary( fit2b ) ### another way fit2c <- lme( fev1 ~ age0 + ageL + female*ageL + factor(f508)*ageL, method = "ML", random = reStruct( ~ 1 | id, pdClass="pdSymm", REML=F), correlation = corCAR1( form = ~ ageL | id ), data = cfkids ) summary( fit2c ) fit2 <- fit2b ### Intercept plus AR(1) plus measurement error fit3 <- lme( fev1 ~ age0 + ageL + female*ageL + factor(f508)*ageL, method = "ML", random = reStruct( ~ 1 | id, pdClass="pdSymm", REML=F), correlation = corExp( form = ~ ageL | id, nugget=T), data = cfkids ) summary( fit3 ) # ##### compare these models # anova( fit0, fit1, fit2, fit3 ) # ##### ##### Residual Analysis -- using fit3 ##### # pop.res <- resid( fit3, level=0 ) cluster.res <- resid( fit3, level=1 ) print( var( pop.res ) ) print( var( cluster.res ) ) # postscript( file="cfkids-NewResiduals.ps", horiz=F ) par( mfrow=c(2,1) ) plot( cfkids$age0, pop.res, pch="." ) lines( smooth.spline( cfkids$age0, pop.res, df=5 ) ) title("Residuals (pop) vs Age0") abline( h=0, lty=2 ) plot( cfkids$ageL, pop.res, pch="." ) lines( smooth.spline( cfkids$ageL, pop.res, df=5 ) ) title("Residuals (pop) vs AgeL") abline( h=0, lty=2 ) graphics.off() # postscript( file="cfkids-NewResiduals2.ps", horiz=F ) par( mfrow=c(2,1) ) plot( cfkids$age0, cluster.res, pch="." ) lines( smooth.spline( cfkids$age0, cluster.res, df=5 ) ) abline( h=0, lty=2 ) title("Residuals (cluster) vs Age0") b0 <- unlist( fit2$coefficients$random ) age0 <- unlist( lapply( split( cfkids$age0, cfkids$id ), min ) ) plot( age0, b0 ) lines( smooth.spline( age0, b0, df=5 ) ) abline( h=0, lty=2 ) title("EB b0 versus Age0") graphics.off() # ##### Do we need a quadratic age0??? # fit4 <- lme( fev1 ~ age0 + age0^2 + ageL + female*ageL + factor(f508)*ageL, method = "ML", random = reStruct( ~ 1 | id, pdClass="pdSymm", REML=F), correlation = corExp( form = ~ ageL | id, nugget=T), data = cfkids ) summary( fit4 ) # anova( fit3, fit4 ) # # end-of-file...