# # hivnet-CDA-gee.q # # ------------------------------------------------------------ # # PURPOSE: analysis of IC data # # AUTHOR: P. Heagerty # # DATE: 00/05/09 # # ------------------------------------------------------------ # data <- read.table( "HivnetWide.dat", header=F ) # vps.data <- data.frame( id = data[,1], risk.group = factor( data[,2], levels=1:4, labels=c("MSM","MaleIDU","WAHR","WAHR+IDU") ), education = factor( data[,3], levels=1:6, labels=c(" 0 ) # ################################################## # Pre/Post Analysis of month 0 and month 6 # ################################################## # ##### gee analysis # library(gee) options( contrasts="contr.treatment" ) # fit1a <- gee( y ~ visit6 + ICgroup + visit6*ICgroup, family=binomial, id=id, corstr="independence", subset=(visit<=1), data=vps.stacked ) summary( fit1a ) # fit1b <- gee( y ~ visit6 + ICgroup + visit6*ICgroup, family=binomial, id=id, corstr="exchangeable", subset=(visit<=1), data=vps.stacked ) summary( fit1b ) # ################################################## # Longitudinal Analysis of 0, 6, and 12 # ################################################## # fit2a <- gee( y ~ visit6 + visit12 + ICgroup + visit6*ICgroup + visit12*ICgroup, family=binomial, id=id, corstr="unstructured", data=vps.stacked ) summary( fit2a ) # fit2b <- gee( y ~ post + visit12 + ICgroup + post*ICgroup + visit12*ICgroup, family=binomial, id=id, corstr="unstructured", data=vps.stacked ) summary( fit2b ) # ##### Adjusted for demographics / risk # vps.stacked$risk.group <- factor( rep(data[,2],rep(3,nsubjects)), levels=1:4, labels=c("MSM","MaleIDU","WAHR","WAHR+IDU") ) vps.stacked$education <- factor( rep(data[,3],rep(3,nsubjects)), levels=1:6, labels=c("