# R code for Xerophthalmia data (Diggle et al., 2002)
# VHM 882, Winter 2005, lectures 11 and 12

xerop <- read.table("e:/vhm/vhm882/data/xerop.dat",header=F,col.names=c("id","resp","int","age","xero","cos","sin","sex","height","stunt","time","age0","seas","time2"),row.names=NULL)

# model 1
xerop.glm1 <- glm(resp~sex+height+xero+age+I(age^2), data=xerop, subset=(time==1), family=binomial)

# alr (using alr package; note that present version only works with R version <2.0)
library(alr)
attach(xerop)

# model 2
X2 <- matrix(c(int,sex,height,cos,sin,xero,age,age^2),ncol=8,byrow=F)
xerop.alr2 <- alr(resp~X1-1,id=id,depmodel="exchangeable",ainit=0.01)
summary(xerop.alr2)

# model 3
X3 <- matrix(c(int,sex,height,xero,age0,age0^2,age-age0,(age-age0)^2),ncol=8,byrow=F)
xerop.alr3 <- alr(resp~X3-1,id=id,depmodel="exchangeable",ainit=0.01)
summary(xerop.alr3)

# model 4
X4 <- matrix(c(int,sex,height,cos,sin,xero,age0,age0^2,age-age0,(age-age0)^2),ncol=10,byrow=F)
xerop.alr4 <- alr(resp~X4-1,id=id,depmodel="exchangeable",ainit=0.01)
summary(xerop.alr4)

# ordinary gee analysis
library(gee)
xerop.gee2 <- gee(resp~sex+height+cos+sin+xero+age+I(age^2), data=xerop, family=binomial, id=id, corstr="exchangeable")
summary(xerop.gee2)
xerop.gee3 <- gee(resp~sex+height+xero+age0+I(age0^2)+I(age-age0)+I((age-age0)^2), data=xerop, family=binomial, id=id, corstr="exchangeable")
summary(xerop.gee3)
xerop.gee4 <- gee(resp~sex+height+cos+sin+xero+age0+I(age0^2)+I(age-age0)+I((age-age0)^2), data=xerop, family=binomial, id=id, corstr="exchangeable")
summary(xerop.gee4)

xerop.gee4u <- gee(resp~sex+height+cos+sin+xero+age0+I(age0^2)+I(age-age0)+I((age-age0)^2), data=xerop, family=binomial, id=id, corstr="unstructured")
summary(xerop.gee4u)

# transitional models
xerop$prev <- xerop$resp-c(NA,diff(xerop$resp))
xerop$difftime <- c(NA,diff(xerop$time))
xerop$diffid <- c(NA,diff(xerop$id))
xerot <- xerop[xerop$diffid==0 & xerop$difftime==1,,]
# Tables 10.1-3, note a discrepancy of 1 additional obs in the dataset
table(xerot$prev,xerot$resp)
table(xerot$xero,xerot$resp)
table(xerot$xero,xerot$resp,xerot$prev)
# model 1 (Table 10.4)
xerop.tr1 <- gee(resp~xero, data=xerot, family=binomial, id=id, corstr="independence")
summary(xerop.tr1)
# model 2 (Table 10.4)
xerop.tr2 <- gee(resp~xero*prev, data=xerot, family=binomial, id=id, corstr="independence")
summary(xerop.tr2)
# model 3 (Table 10.4)
xerop.tr3 <- gee(resp~xero+prev, data=xerot, family=binomial, id=id, corstr="independence")
summary(xerop.tr3)
# model 4 (Table 10.4)
xerop.tr4 <- gee(resp~xero*prev+age*prev+I(seas==2)*prev, data=xerot, family=binomial, id=id, corstr="independence")
summary(xerop.tr4)
# model 5 (Table 10.4)
xerop.tr5 <- gee(resp~xero+prev+age+I(seas==2), data=xerot, family=binomial, id=id, corstr="independence")
summary(xerop.tr5)

# time-2 lagged models
xerop$prev2 <- xerop$resp-c(NA,NA,diff(xerop$resp,2))
xerop$difftime2 <- c(NA,NA,diff(xerop$time,2))
xerop$diffid2 <- c(NA,NA,diff(xerop$id,2))
xerot2 <- xerop[xerop$diffid==0 & xerop$difftime==1 & xerop$diffid2==0 & xerop$difftime2==2,,]
# model 6
xerop.tr6 <- gee(resp~xero+age+I(seas==2), data=xerot2, family=binomial, id=id, corstr="independence")
summary(xerop.tr6)
# model 7 
xerop.tr7 <- gee(resp~xero+age+I(seas==2)+prev, data=xerot2, family=binomial, id=id, corstr="independence")
summary(xerop.tr7)
# model 8 
xerop.tr8 <- gee(resp~xero+age+I(seas==2)+prev+prev2, data=xerot2, family=binomial, id=id, corstr="independence")
summary(xerop.tr8)
# models do not show the behaviour descripted in Section 10.3.1