# R code for Example 9.1 (Table 9.1) in (Davis, 2002)
# VHM 882, Winter 2005, lecture 10L
# note: gee implementation in gee library deals with missing values
#       in a different way than SAS for unstructured correlations

#m <- matrix(scan("d:\\data.avc\\teaching\\vhm882\\10l\\davis5r.dat",what="character",na.strings="."),ncol=7,byrow=T)
m <- matrix(scan("e:\\vhm\\vhm882\\data\\davis5r.dat",what="character",na.strings="."),ncol=7,byrow=T)
week <- c(2,4,8)
subject <- type.convert(m[,1])
sex <- type.convert(m[,4])
age <- type.convert(m[,3])
dose <- type.convert(m[,2])
#cgr <- type.convert(m[,5:7])

library(nlme)
cgr <- balancedGrouped( cgr ~ week|subject, matrix(type.convert(m[,5:7]),nrow=75,ncol=3,dimnames=list(subject,week)))
Dose <- rep(dose,rep(3,75))
Sex <- rep(sex,rep(3,75))
Age <- rep(age,rep(3,75))

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

# model 1
cgr.gee1 <- gee(cgr~as.factor(Dose)*as.factor(week)+as.factor(Sex)*as.factor(week)+Age*as.factor(week), data=cgr, family=binomial, id=subject, corstr="unstructured")
summary(cgr.gee1)

# model 2
cgr.gee2 <- update(cgr.gee1, cgr~as.factor(Dose)+as.factor(week)+as.factor(Sex)+Age)
summary(cgr.gee2)
# Wald test for week
bhat <- cgr.gee2[9]$coefficients
rvar <- cgr.gee2[20]$robust.variance
c.week <- matrix(c(rep(0,4),1,rep(0,8),1,rep(0,2)),nrow=2,ncol=8,byrow=T)
library(MASS)
wald.week <- t(c.week %*% bhat) %*% ginv(c.week %*% rvar %*% t(c.week)) %*% c.week %*% bhat
wald.week
pchisq(wald.week,2,lower.tail=F)

# model 3
cgr.gee3 <- update(cgr.gee1, cgr~as.factor(Dose)+as.factor(Sex)+Age)
summary(cgr.gee3)
# Wald test for dose
bhat <- cgr.gee3[9]$coefficients
rvar <- cgr.gee3[20]$robust.variance
c.dose <- matrix(c(0,1,rep(0,6),1,rep(0,6),1,rep(0,2)),nrow=3,ncol=6,byrow=T)
wald.dose <- t(c.dose %*% bhat) %*% ginv(c.dose %*% rvar %*% t(c.dose)) %*% c.dose %*% bhat
wald.dose
pchisq(wald.dose,3,lower.tail=F)

# exchangeable working correlation
cgr.gee3cs <- update(cgr.gee3, corstr="exchangeable")
summary(cgr.gee3cs)
# Wald test for dose
bhat <- cgr.gee3cs[9]$coefficients
rvar <- cgr.gee3cs[20]$robust.variance
c.dose <- matrix(c(0,1,rep(0,6),1,rep(0,6),1,rep(0,2)),nrow=3,ncol=6,byrow=T)
wald.dose.cs <- t(c.dose %*% bhat) %*% ginv(c.dose %*% rvar %*% t(c.dose)) %*% c.dose %*% bhat
wald.dose.cs
pchisq(wald.dose.cs,3,lower.tail=F)

# independence working correlation
cgr.gee3ind <- update(cgr.gee3, corstr="independence")
summary(cgr.gee3ind)

# model 4
cgr.gee4 <- update(cgr.gee1, cgr~I(Dose/250)+as.factor(Sex)+Age)
summary(cgr.gee4)