# R code for 2x2 crossover trial (Diggle et al, 2002)
# VHM 882, Winter 2005, lectures 10 and 11

xover2 <- read.table("e:/vhm/vhm882/data/xover2.dat",header=F,col.names=c("id","int1","resp","int2","tx","per","txper","groupBA"),row.names=NULL)
xover2$tx <- 1-xover2$tx             # switch to indicator for A
xover2$txper <- xover2$tx*xover2$per # switch to indicator for A

# conditional logistic regression (in survival package)
library(splines)
library(survival)
xover2.cl1 <- clogit(resp~tx+strata(id), data=xover2)
summary(xover2.cl1)
xover2.cl2 <- clogit(resp~tx+per+strata(id), data=xover2)
summary(xover2.cl2)
xover2.cl3 <- clogit(resp~tx+per+txper+strata(id), data=xover2)
summary(xover2.cl3)

# gee (using gee library from gee package)
library(gee)
xover2.gee1 <- gee(resp~tx, data=xover2, family=binomial, id=id, corstr="exchangeable")
summary(xover2.gee1)
xover2.gee2 <- update(xover2.gee1, resp~tx+per)
summary(xover2.gee2)
xover2.gee3 <- update(xover2.gee1, resp~tx+per+txper)
summary(xover2.gee3)
xover2.gee4 <- update(xover2.gee2, corstr="independence")
summary(xover2.gee4)

# alr (using alr package, seems to work only with versions <2.0.0)
library(alr)
attach(xover2)
X1 <- matrix(c(int1,tx,per),nrow=134,ncol=3,byrow=F)
xover2.alr1 <- alr(resp~X1-1,id=id,depmodel="exchangeable",ainit=0.01)
summary(xover2.alr1)
X2 <- matrix(c(int1,tx,per,txper),nrow=134,ncol=4,byrow=F)
xover2.alr2 <- alr(resp~X2-1,id=id,depmodel="exchangeable",ainit=0.01)
summary(xover2.alr2)