# # cfkids-CDA-gee.q # # ------------------------------------------------------------ # # PURPOSE: Use GEE to characterize longitudinal # change by gender and genotype. # # AUTHOR: P. Heagerty # # DATE: 00/07/10 # # ------------------------------------------------------------ # # ##### ##### Read data ##### # source("cfkids-read.q") # # ##### ##### GEE ##### # options( contrasts="contr.treatment" ) library( gee ) # fit1a <- gee( fev1 ~ age0 + ageL + female + factor(f508) + female * ageL + factor(f508) * ageL, id = id, corstr="independence", data = cfkids ) summary( fit1a ) # fit1b <- gee( fev1 ~ age0 + ageL + female + factor(f508) + female * ageL + factor(f508) * ageL, id = id, corstr="exchangeable", data = cfkids ) summary( fit1b ) # ##### ##### Residual Analysis ##### # resids <- residuals( fit1b ) # attach( cfkids ) # ##### Residuals versus Age0, ageL # postscript( file="cfkids-CDA-gee-resids1.ps", horiz=F ) par( mfrow=c(2,1) ) # plot( age0, resids, pch="." ) abline( h=0, lty=3 ) lines( smooth.spline( age0, resids, df=7 ) ) title("Residuals versus Age0") # plot( ageL, resids, pch="." ) abline( h=0, lty=3 ) lines( smooth.spline( ageL, resids, df=7 ) ) title("Residuals versus AgeL") # graphics.off() # ##### Residuals versus AgeL for male / female # postscript( file="cfkids-CDA-gee-resids2.ps", horiz=F ) par( mfrow=c(2,1) ) # females <- female==1 plot( ageL[females], resids[females], pch="." ) abline( h=0, lty=3 ) lines( smooth.spline( ageL[females], resids[females], df=7 ) ) title("Residuals versus AgeL - Females") # males <- female==0 plot( ageL[males], resids[males], pch="." ) abline( h=0, lty=3 ) lines( smooth.spline( ageL[males], resids[males], df=7 ) ) title("Residuals versus AgeL - Males") # graphics.off() # postscript( file="cfkids-CDA-gee-resids3.ps", horiz=F ) par( mfrow=c(3,1) ) # f508.is.0 <- f508==0 plot( ageL[f508.is.0], resids[f508.is.0], pch="." ) abline( h=0, lty=3 ) lines( smooth.spline( ageL[f508.is.0], resids[f508.is.0], df=7 ) ) title("Residuals versus AgeL - F508=0") # f508.is.1 <- f508==1 plot( ageL[f508.is.1], resids[f508.is.1], pch="." ) abline( h=0, lty=3 ) lines( smooth.spline( ageL[f508.is.1], resids[f508.is.1], df=7 ) ) title("Residuals versus AgeL - F508=1") # f508.is.2 <- f508==2 plot( ageL[f508.is.2], resids[f508.is.2], pch="." ) abline( h=0, lty=3 ) lines( smooth.spline( ageL[f508.is.2], resids[f508.is.2], df=7 ) ) title("Residuals versus AgeL - F508=2") # graphics.off() # # # end-of-file...