# R code for Exercise 9.11 in (Davis, 2002)

m <- matrix(scan("d:\\data.ext\\davis\\fl95.dat",what="character",na.strings="."),ncol=9,byrow=T)
#m <- matrix(scan("e:\\vhm\\vhm882\\data\\fl95.dat",what="character",na.strings="."),ncol=9,byrow=T)
week <- c(0,1,5,9,13)
subject <- type.convert(m[,1])
tx <- type.convert(m[,4])
sex <- type.convert(m[,2])
age <- type.convert(m[,3])

library(nlme)
arth <- balancedGrouped( arth ~ week|subject, matrix(type.convert(m[,5:9]),nrow=51,ncol=5,dimnames=list(subject,week)))
Sex <- rep(sex,rep(5,51))
Age <- rep(age,rep(5,51))
Tx <- rep(tx,rep(5,51))
Subject <- rep(subject,rep(5,51))
arth$activetx <- arth$week>1 & Tx=="A"

# (a)
# including week, sex, age, treatment and using unstructured working correlations
# gee (using gee library from gee package)
library(gee)
arth.gee1.uns <- gee(arth~as.factor(week)+as.factor(activetx)+as.factor(Sex)+Age, data=arth, family=binomial, id=subject, corstr="unstructured")
summary(arth.gee1.uns)
arth.gee1.exch <- update(arth.gee1.uns,corstr="exchangeable")
summary(arth.gee1.exch)
arth.gee1.ind <- update(arth.gee1.uns,corstr="independence")
summary(arth.gee1.ind)

# alr (using alr library from alr package)
library(alr)
int <- rep(1,255)
X1 <- matrix(c(int,as.numeric(arth$week==1),as.numeric(arth$week==5),as.numeric(arth$week==9),as.numeric(arth$week==13),as.numeric(arth$activetx),as.numeric(Sex=="M"),Age),nrow=255,ncol=8,byrow=F)
arth.alr1.exch <- alr(arth$arth~X1-1,id=Subject,depmodel="exchangeable",ainit=0.01)
summary(arth.alr1.exch)
ZIND <- rep(1:5,51)
ZMAST <- matrix(rep(0,200),nrow=20,ncol=10)
# form of ZMAST matrix corresponding to 5x5 correlation matrix:
# X 1 2 3 4
# 1 X 5 6 7
# 2 5 X 8 9
# 3 6 8 X 10
# 4 7 9 10 X
ZMAST[1,1] <- 1
ZMAST[2,2] <- 1
ZMAST[3,3] <- 1
ZMAST[4,4] <- 1
ZMAST[5,1] <- 1
ZMAST[6,5] <- 1
ZMAST[7,6] <- 1
ZMAST[8,7] <- 1
ZMAST[9,2] <- 1
ZMAST[10,5] <- 1
ZMAST[11,8] <- 1
ZMAST[12,9] <- 1
ZMAST[13,3] <- 1
ZMAST[14,6] <- 1
ZMAST[15,8] <- 1
ZMAST[16,10] <- 1
ZMAST[17,4] <- 1
ZMAST[18,7] <- 1
ZMAST[19,9] <- 1
ZMAST[20,10] <- 1
# does not give meaningful results
arth.alr1.uns <- alr(arth$arth~X1-1,id=Subject,depmodel="general",z=ZMAST,zmast=1,zloc=ZIND,ainit=rep(0.01,10))
summary(arth.alr1.uns)
# works (agreement with SAS)
arth.alr1.uns <- alr(arth$arth[!is.na(arth$arth)]~X1[!is.na(arth$arth),]-1,id=Subject[!is.na(arth$arth)],depmodel="general",z=ZMAST,zmast=1,zloc=ZIND[!is.na(arth$arth)],ainit=rep(0.5,10))
summary(arth.alr1.uns)

# (b)
# some differences are seen between estimates, depending on approach (GEE/ALR) and structure used
# a cautious approach is to use GEE with independence correlation structure,
# because correlations are variable with no obvious pattern, and the dataset may
# be too small to use the full correlation matrix
arth.gee2.ind <- update(arth.gee1.ind,arth~as.factor(week)+as.factor(activetx)+as.factor(Sex)+Age+I(Age^2))
summary(arth.gee2.ind)
arth.gee3.ind <- update(arth.gee1.ind,arth~as.factor(week)+as.factor(activetx)+as.factor(Sex)+Age+as.factor(activetx)*Age)
summary(arth.gee3.ind)
arth.gee4.ind <- update(arth.gee1.ind,arth~as.factor(week)+as.factor(activetx)+as.factor(Sex)+Age+as.factor(activetx)*as.factor(Sex))
summary(arth.gee4.ind)
# Wald test for week
bhat <- arth.gee1.ind$"coefficients"
rvar <- arth.gee1.ind$"robust.variance"
c.week <- matrix(c(0,1,rep(0,8),1,rep(0,8),1,rep(0,8),1,rep(0,3)),nrow=4,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,4,lower.tail=F)
# Wald test for comparison post-pre (post placebo only)
c.prepost <- matrix(c(0,0.5,-1/3,-1/3,-1/3,0,0,0),nrow=1,ncol=8,byrow=T)
library(MASS)
wald.prepost <- t(c.prepost %*% bhat) %*% ginv(c.prepost %*% rvar %*% t(c.prepost)) %*% c.prepost %*% bhat
wald.prepost
pchisq(wald.prepost,1,lower.tail=F)
# conclusion: all other effects than activetx seem non-significant
# strictly speaking, a backwards elimination should be performed
arth.gee5.ind <- update(arth.gee1.ind,arth~as.factor(activetx))
summary(arth.gee5.ind)
# a clear treatment effect, corresponding to OR=3.9