R version 2.4.0 (2006-10-03) Copyright (C) 2006 The R Foundation for Statistical Computing ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > # > # hivnet-CDA-geese.q > # > # ------------------------------------------------------------ > # > # PURPOSE: analysis of IC data > # > # AUTHOR: P. Heagerty > # > # DATE: 02Mar2007 > # > # ------------------------------------------------------------ > # > 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(" # > ################################################## > # Longitudinal data # > ################################################## > y0 <- vps.data$q4safe0 > y6 <- vps.data$q4safe6 > y12 <- vps.data$q4safe12 > nsubjects <- length(y0) > # > ##### stacked data for regression > # > vps.stacked <- data.frame( + y = as.vector( rbind( y0, y6, y12 ) ), + visit = rep( c(0,1,2), nsubjects ), + ICgroup = as.vector( rbind( vps.data$ICgroup, vps.data$ICgroup, + vps.data$ICgroup ) ), + id = rep( 1:nsubjects, rep(3,nsubjects) ) + ) > vps.stacked$visit6 <- as.integer( vps.stacked$visit==1 ) > vps.stacked$visit12 <- as.integer( vps.stacked$visit==2 ) > vps.stacked$post <- as.integer( vps.stacked$visit > 0 ) > # > ################################################## > # Pre/Post Analysis of month 0 and month 6 # > ################################################## > # > ##### gee analysis > # > library(geepack) > # > fit1a <- geese( y ~ visit6 + ICgroup + visit6*ICgroup, + family=binomial, + id=id, + corstr="independence", + subset=(visit<=1), data=vps.stacked ) > summary( fit1a ) Call: geese(formula = y ~ visit6 + ICgroup + visit6 * ICgroup, id = id, data = vps.stacked, subset = (visit <= 1), family = binomial, corstr = "independence") Mean Model: Mean Link: logit Variance to Mean Relation: binomial Coefficients: estimate san.se wald p (Intercept) 0.25741201 0.09018456 8.14691533 0.004313447 visit6 -0.06481890 0.09854411 0.43265525 0.510688927 ICgroup 0.01628382 0.12765263 0.01627246 0.898494378 visit6:ICgroup 0.36648722 0.14433417 6.44732200 0.011111971 Scale Model: Scale Link: identity Estimated Scale Parameters: estimate san.se wald p (Intercept) 1 0.009667174 10700.42 0 Correlation Model: Correlation Structure: independence Returned Error Value: 0 Number of clusters: 1000 Maximum cluster size: 2 > # > fit1b <- geese( y ~ visit6 + ICgroup + visit6*ICgroup, + family=binomial, + id=id, corstr="exchangeable", + subset=(visit<=1), data=vps.stacked ) > summary( fit1b ) Call: geese(formula = y ~ visit6 + ICgroup + visit6 * ICgroup, id = id, data = vps.stacked, subset = (visit <= 1), family = binomial, corstr = "exchangeable") Mean Model: Mean Link: logit Variance to Mean Relation: binomial Coefficients: estimate san.se wald p (Intercept) 0.25741201 0.09018456 8.14691533 0.004313447 visit6 -0.06481890 0.09854411 0.43265525 0.510688927 ICgroup 0.01628382 0.12765263 0.01627246 0.898494378 visit6:ICgroup 0.36648722 0.14433417 6.44732200 0.011111971 Scale Model: Scale Link: identity Estimated Scale Parameters: estimate san.se wald p (Intercept) 1 0.009667174 10700.42 0 Correlation Model: Correlation Structure: exchangeable Correlation Link: identity Estimated Correlation Parameters: estimate san.se wald p alpha 0.3696618 0.03007344 151.0925 0 Returned Error Value: 0 Number of clusters: 1000 Maximum cluster size: 2 > # > ################################################## > # Longitudinal Analysis of 0, 6, and 12 # > ################################################## > # > fit2a <- geese( y ~ visit6 + visit12 + ICgroup + + visit6*ICgroup + visit12*ICgroup, + family=binomial, + id=id, corstr="unstructured", + data=vps.stacked ) > summary( fit2a ) Call: geese(formula = y ~ visit6 + visit12 + ICgroup + visit6 * ICgroup + visit12 * ICgroup, id = id, data = vps.stacked, family = binomial, corstr = "unstructured") Mean Model: Mean Link: logit Variance to Mean Relation: binomial Coefficients: estimate san.se wald p (Intercept) 0.25741201 0.09018456 8.14691533 0.004313447 visit6 -0.06481890 0.09854411 0.43265525 0.510688927 visit12 0.08180371 0.11093456 0.54376634 0.460876011 ICgroup 0.01628382 0.12765263 0.01627246 0.898494378 visit6:ICgroup 0.36648722 0.14433417 6.44732200 0.011111971 visit12:ICgroup 0.24600305 0.15531956 2.50858397 0.113227632 Scale Model: Scale Link: identity Estimated Scale Parameters: estimate san.se wald p (Intercept) 1 0.01028674 9450.27 0 Correlation Model: Correlation Structure: unstructured Correlation Link: identity Estimated Correlation Parameters: estimate san.se wald p alpha.1:2 0.3696618 0.03002564 151.5739 0 alpha.1:3 0.2740313 0.03096106 78.3373 0 alpha.2:3 0.3901942 0.03029926 165.8431 0 Returned Error Value: 0 Number of clusters: 1000 Maximum cluster size: 3 > # > fit2b <- geese( y ~ post + visit12 + ICgroup + + post*ICgroup + visit12*ICgroup, + family=binomial, + id=id, corstr="unstructured", + data=vps.stacked ) > summary( fit2b ) Call: geese(formula = y ~ post + visit12 + ICgroup + post * ICgroup + visit12 * ICgroup, id = id, data = vps.stacked, family = binomial, corstr = "unstructured") Mean Model: Mean Link: logit Variance to Mean Relation: binomial Coefficients: estimate san.se wald p (Intercept) 0.25741201 0.09018456 8.14691533 0.004313447 post -0.06481890 0.09854411 0.43265525 0.510688927 visit12 0.14662262 0.10356299 2.00443857 0.156839352 ICgroup 0.01628382 0.12765263 0.01627246 0.898494378 post:ICgroup 0.36648722 0.14433417 6.44732200 0.011111971 visit12:ICgroup -0.12048417 0.14337561 0.70617026 0.400718132 Scale Model: Scale Link: identity Estimated Scale Parameters: estimate san.se wald p (Intercept) 1 0.01028674 9450.27 0 Correlation Model: Correlation Structure: unstructured Correlation Link: identity Estimated Correlation Parameters: estimate san.se wald p alpha.1:2 0.3696618 0.03002564 151.5739 0 alpha.1:3 0.2740313 0.03096106 78.3373 0 alpha.2:3 0.3901942 0.03029926 165.8431 0 Returned Error Value: 0 Number of clusters: 1000 Maximum cluster size: 3 > # > ##### 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(" vps.stacked$age <- rep( data[,4], rep(3,nsubjects) ) > vps.stacked$cohort <- rep( data[,5], rep(3,nsubjects) ) > # > fit2c <- geese( y ~ post + visit12 + ICgroup + + post*ICgroup + visit12*ICgroup + age + education + + cohort + risk.group, + family=binomial, + id=id, corstr="unstructured", + data=vps.stacked ) > summary( fit2c ) Call: geese(formula = y ~ post + visit12 + ICgroup + post * ICgroup + visit12 * ICgroup + age + education + cohort + risk.group, id = id, data = vps.stacked, family = binomial, corstr = "unstructured") Mean Model: Mean Link: logit Variance to Mean Relation: binomial Coefficients: estimate san.se wald p (Intercept) -0.653812235 0.303400434 4.6438020 3.116588e-02 post -0.072804833 0.110246507 0.4361047 5.090086e-01 visit12 0.164608064 0.115911699 2.0167296 1.555739e-01 ICgroup 0.057407302 0.135428242 0.1796866 6.716427e-01 age 0.002867246 0.006516176 0.1936174 6.599228e-01 educationHS 0.532181933 0.171119647 9.6720938 1.870880e-03 educationsome college 1.078366082 0.191577169 31.6843596 1.813770e-08 educationcollege 1.474141590 0.212775156 47.9994833 4.263367e-12 educationsome post 1.465868394 0.278698552 27.6643543 1.442973e-07 educationgrad/prof 1.719539978 0.240970578 50.9209522 9.615642e-13 cohort -0.176522778 0.104817255 2.8361940 9.216227e-02 risk.groupMaleIDU -0.358384469 0.170189288 4.4343903 3.522180e-02 risk.groupWAHR -0.280944024 0.199440810 1.9843192 1.589361e-01 risk.groupWAHR+IDU -0.579334066 0.206740864 7.8524554 5.075173e-03 post:ICgroup 0.413160363 0.162375803 6.4743277 1.094437e-02 visit12:ICgroup -0.135179205 0.161215903 0.7030788 4.017511e-01 Scale Model: Scale Link: identity Estimated Scale Parameters: estimate san.se wald p (Intercept) 1.001947 0.02805194 1275.745 0 Correlation Model: Correlation Structure: unstructured Correlation Link: identity Estimated Correlation Parameters: estimate san.se wald p alpha.1:2 0.2919784 0.03382777 74.49974 0.000000e+00 alpha.1:3 0.1854308 0.03343967 30.74962 2.935621e-08 alpha.2:3 0.3110622 0.03454029 81.10405 0.000000e+00 Returned Error Value: 0 Number of clusters: 1000 Maximum cluster size: 3 > # > # > # end-of-file... > > > > > > > > >