# R code for Exercise 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)))
Site <- rep(site,rep(4,111))
Sex <- rep(sex,rep(4,111))
Age <- rep(age,rep(4,111))
Tx <- rep(tx,rep(4,111))
Base <- rep(base,rep(4,111))

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

# including visit, unstructured working correlations
resp.gee1 <- gee(I(resp>2)~as.factor(Site)+as.factor(Tx)+as.factor(Sex)+Age+as.factor(visit)+I(Base>2), data=resp, family=binomial, id=subject, corstr="unstructured")
summary(resp.gee1)
# Wald test for visit
bhat <- resp.gee1[9]$coefficients
rvar <- resp.gee1[20]$robust.variance
c.visit <- matrix(c(rep(0,5),1,rep(0,9),1,rep(0,9),1,0),nrow=3,ncol=9,byrow=T)
library(MASS)
wald.visit <- t(c.visit %*% bhat) %*% ginv(c.visit %*% rvar %*% t(c.visit)) %*% c.visit %*% bhat
wald.visit
pchisq(wald.visit,3,lower.tail=F)

# excluding visit, unstructured working correlations
resp.gee2 <- 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.gee2)

# stationary working correlations
resp.gee2.stat <- update(resp.gee2, corstr="stat_M_dep",Mv=3)
summary(resp.gee2.stat)

# exchangeable working correlations
resp.gee2.cs <- update(resp.gee2, corstr="exchangeable")
summary(resp.gee2.cs)

# independent working correlations
resp.gee2.ind <- update(resp.gee2, corstr="independence")
summary(resp.gee2.ind)

# comments:
# minor differences in estimates between working correlation structures,
# most sensible seems to be the stationary structure
# OR's: site (2 vs 1) 1.95, tx (A vs P) 3.49, sex (M vs F) 0.87 NS, 
#       age (10 yrs) 0.83, baseline (good vs poor) 6.37