# R code for Additional exercise 12.1 (addition to Exericse 9.5 in Davis,2002)

#m <- matrix(scan("d:\\data.ext\\davis\\davis1.dat",what="character",na.strings="."),ncol=10,byrow=T)
m <- matrix(scan("e:\\vhm\\vhm882\\data\\davis1.dat",what="character",na.strings="."),ncol=10,byrow=T)
visit <- c(1:4)
site <- type.convert(m[,1])
subject <- type.convert(m[,2])+(site-1)*100
tx <- type.convert(m[,3])
sex <- type.convert(m[,4])
age <- type.convert(m[,5])
base <- type.convert(m[,6])

library(nlme)
resp <- balancedGrouped( resp ~ visit|subject, matrix(type.convert(m[,7:10]),nrow=111,ncol=4,dimnames=list(subject,visit)))
resp$Site <- rep(site,rep(4,111))
resp$Sex <- rep(sex,rep(4,111))
resp$Age <- rep(age,rep(4,111))
resp$Tx <- rep(tx,rep(4,111))
resp$Base <- rep(base,rep(4,111))
resp$prev <- resp$resp-c(NA,diff(resp$resp))
for (i in 1:444) { if (resp$visit[i]==1) resp$prev[i]<- resp$Base[i] }

# gee (using gee library from gee package)
library(gee)

# marginal model
resp.gee0 <- gee(I(resp>2)~as.factor(Site)+as.factor(Tx)+as.factor(Sex)+Age+I(Base>2), data=resp, family=binomial, id=subject, corstr="unstructured")
summary(resp.gee0)

# transitional model, no interactions
resp.gee1 <- gee(I(resp>2)~as.factor(Site)+as.factor(Tx)+as.factor(Sex)+Age+I(prev>2), data=resp, family=binomial, id=subject, corstr="unstructured")
summary(resp.gee1)
# note: smaller working correlations, smaller Tx effect

# transitional model, interaction with Tx - gives an error
resp.gee2 <- gee(I(resp>2)~as.factor(Site)+as.factor(Tx)*I(prev>2)+as.factor(Sex)+Age, data=resp, family=binomial, id=subject, corstr="unstructured")
summary(resp.gee2)

# transitional model, interaction with Tx - stationary working correlations
resp.gee2.stat <- update(resp.gee2, corstr="stat_M_dep",Mv=3)
summary(resp.gee2.stat)
# no evidence of Tx by prev interaction

# transitional model, interaction with Site - stationary working correlations
resp.gee3.stat <- update(resp.gee2.stat, I(resp>2)~as.factor(Site)*I(prev>2)+as.factor(Tx)+as.factor(Sex)+Age)
summary(resp.gee3.stat)
# no evidence of Site by prev interaction

# transitional model, no interactions, stationary working correlation
resp.gee1.stat <- update(resp.gee1, corstr="stat_M_dep",Mv=3)
summary(resp.gee1.stat)

# transitional model, no interactions, independence working correlation
resp.gee1.ind <- update(resp.gee1, corstr="independence")
summary(resp.gee1.ind)
# very similar results of all three working correlation structures
# OR's (stationary working correlation): 
#     site (2 vs 1) 2.07, tx (A vs P) 2.35, sex (M vs F) 0.95 NS, 
#     age (10 yrs) 0.83, prev (good vs poor) 9.35
# about 2/3 for Tx, a bit larger for site, and value for prev 3/2 of baseline