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

#m <- matrix(scan("d:\\data.ext\\davis\\kz.dat",what="character",na.strings="."),ncol=8,byrow=T)
m <- matrix(scan("e:\\vhm\\vhm882\\data\\kz.dat",what="character",na.strings="."),ncol=8,byrow=T)
quarter <- c(1:4)
subject <- type.convert(m[,1])
age <- type.convert(m[,2])
sex <- type.convert(m[,3])
smoking <- type.convert(m[,4])

library(nlme)
visits <- balancedGrouped( visits ~ quarter|subject, matrix(type.convert(m[,5:8]),nrow=73,ncol=4,dimnames=list(subject,quarter)))
visits$Sex <- rep(sex,rep(4,73))
visits$Age <- rep(age,rep(4,73))
visits$Smoking <- rep(smoking,rep(4,73))
summary(visits)

# profile plot
plot(visits, outer=~as.factor(Smoking), aspect=1)
plot(visits, outer=~as.factor(Sex), aspect=1)
plot(visits, outer=~as.factor(Age %/% 12), aspect=1)
# mean plot
visits.means <- aggregate(visits$visit,list(visits$quarter,visits$Smoking),mean,na.rm=TRUE)
names(visits.means)[1]<-"Quarter"
visits.means$Quarter <- as.numeric(levels(visits.means$Quarter))[visits.means$Quarter]
names(visits.means)[2]<-"Smoking"
library(lattice)
xyplot(x~Quarter|Smoking,visits.means,ylab="Mean No. of visits in 4 quarters")
visits.means2 <- aggregate(visits$visit,list(visits$quarter,visits$Age %/% 12),mean,na.rm=TRUE)
names(visits.means2)[1]<-"Quarter"
visits.means2$Quarter <- as.numeric(levels(visits.means2$Quarter))[visits.means2$Quarter]
names(visits.means2)[2]<-"Age"
library(lattice)
xyplot(x~Quarter|Age,visits.means2,ylab="Mean No. of visits in 4 quarters")

# gee (using gee library from gee package)
library(gee)

# including quarter, unstructured working correlations
visits.gee1 <- gee(visits~as.factor(Smoking)+as.factor(Sex)+Age+as.factor(quarter), data=visits, family=poisson, id=subject, corstr="unstructured")
summary(visits.gee1)

# linear effect of quarter, quadratic of Age
visits.gee2 <- gee(visits~as.factor(Smoking)+as.factor(Sex)+Age+I(Age^2)+quarter, data=visits, family=poisson, id=subject, corstr="unstructured")
summary(visits.gee2)

# added interaction with age
visits.gee3 <- gee(visits~as.factor(Smoking)+as.factor(Sex)+I(Age-36)*quarter, data=visits, family=poisson, id=subject, corstr="unstructured")
summary(visits.gee3)

# stationary correlation structure
visits.gee3.stat <- update(visits.gee3, corstr="stat_M_dep",Mv=3)
summary(visits.gee3.stat)

# exchangeable correlation structure
visits.gee3.exch <- update(visits.gee3, corstr="exchangeable",Mv=3)
summary(visits.gee3.exch)

# exchangeable correlation structure
visits.gee3.ind <- update(visits.gee3, corstr="independence",Mv=3)
summary(visits.gee3.ind)
# results quite robust to correlation structure
# however, data still overdispersed (scale=1.8)

# random effects models fit by glmmPQL
library(MASS)
visits.pql3 <- glmmPQL(visits~as.factor(Smoking)+as.factor(Sex)+I(Age-36)*quarter, random=~1|subject,data=visits, family=poisson)
summary(visits.pql3)
# similar results but residual variance only slightly larger than 1

# main significant effect is of quarter: visits become less frequent over time
# no significance of sex and smoking whatsoever
# small age effect manifested as an interaction with quarter: younger children decline less rapidly in their visits over time
# quarter and age may be longitudinal and cross-sectional age effects.