# R code for 3x3 crossover trial (Diggle et al., 2002)
# VHM 882, Winter 2005, lecturea 10 and 11

xover3 <- read.table("e:/vhm/vhm882/data/xover3.dat",header=F,col.names=c("id","per","resp","int","x3","x4","x1","x2","prevA","x5","x6"),row.names=NULL)

# conditional logistic regression (in survival package)
library(splines)
library(survival)
xover3.cl <- clogit(resp~x1+x2+x3+x4+x5+x6+strata(id), data=xover3)
summary(xover3.cl)

# gee (using gee library from gee package)
library(gee)
xover3.gee1 <- gee(resp~x1+x2+x3+x4+x5+x6, data=xover3, family=binomial, id=id, corstr="exchangeable")
summary(xover3.gee1)
xover3.gee2 <- update(xover3.gee1, corstr="unstructured")
summary(xover3.gee2)

# alr (using alr package, seems to work only with versions <2.0.0)
library(alr)
attach(xover3)
X1 <- matrix(c(int,x1,x2,x3,x4,x5,x6),nrow=258,ncol=7,byrow=F)
xover3.alr1 <- alr(resp~X1-1,id=id,depmodel="exchangeable",ainit=0.01)
summary(xover3.alr1)
X2 <- matrix(c(int,x1,x2,x3,x4,x5,x6,x1*x3,x1*x4,x2*x3,x2*x4),nrow=258,ncol=11,byrow=F)
xover3.alr2 <- alr(resp~X2-1,id=id,depmodel="exchangeable",ainit=0.01)
summary(xover3.alr2)
ZIND <- rep(1:3,86)
ZMAST <- matrix(c(1,0,0,0,1,0,1,0,0,0,0,1,0,1,0,0,0,1),nrow=6,ncol=3,byrow=T)
xover3.alr3 <- alr(resp~X1-1,id=id,depmodel="general",z=ZMAST,zmast=1,zloc=ZIND,ainit=rep(0.01,3))
summary(xover3.alr3)
xover3.alr4 <- alr(resp~X2-1,id=id,depmodel="general",z=ZMAST,zmast=1,zloc=ZIND,ainit=rep(0.01,3))
summary(xover3.alr4)

# glmmML (crude implementation of numerical integration); works better for R-version>2.0
library(glmmML)
xover3.ML <- glmmML(resp~x1+x2+x3+x4+x5+x6, data=xover3, cluster=id, family=binomial,n.points=32)
summary(xover3.ML)