#  Fit gee with linear estimating equation for beta and
#  moment methods as in Liang and Zeger (1986) for correlation
#  parameter using the R function gee()

#  load the gee library

library(gee)

#  output data set

outfile <- "gee_seizure.Rout"

#  read in the data set

thedata <- matrix(scan("seize.dat"),ncol=6,byrow=T)

subj <- thedata[,1]
seize <- thedata[,2]
visit <- thedata[,3]
trt <- thedata[,4]
base <- thedata[,5]
age <- thedata[,6]

#  log tranform baseline and age

logbase <- log(base/4)
logage <- log(age)

#  other variables that could be used in more complicated models

basetrt <- logbase*trt
visit4 <- as.numeric(visit==4)

#  calls to gee() function

unfit <- gee(seize ~ logbase+logage+trt+visit4+basetrt,id=subj,family=poisson,
  corstr="unstructured")

cat("UNSTRUCTURED WORKING CORRELATION","\n","\n",file=outfile,append=F)
sink(file=outfile,append=T)
print(summary(unfit))
sink()
cat("\n","\n",file=outfile,append=T)

csfit <- gee(seize ~ logbase+logage+trt+visit4+basetrt,id=subj,family=poisson,
  corstr="exchangeable")

cat("EXCHANGEABLE WORKING CORRELATION","\n","\n",file=outfile,append=T)
sink(file=outfile,append=T)
print(summary(csfit))
sink()
cat("\n","\n",file=outfile,append=T)

#  note demonstration of use of user-supplied starting values here

ar1fit <- gee(seize ~ logbase+logage+trt+visit4+basetrt,id=subj,family=poisson,
   corstr="AR-M",M=1,b=c(-2.0,1,1,-1,-0.1,0.5) )

cat("AR(1) WORKING CORRELATION","\n","\n",file=outfile,append=T)
sink(file=outfile,append=T)
print(summary(ar1fit))
sink()