# Example 6.3 in Davis, 2002

library(nlme)
jb <- groupedData(src~age|subject,data=read.table("e:/vhm/vhm882/data/jb.dat",col.names=c("subject","group","age","src")),
                  labels=list(y="Reciprocal serum creatinine level"))
jb$group <- as.factor(jb$group)

# code for graphics, using nlme and lattice libraries
# profile plot (subset of subjects from each group)
p1 <- plot(jb[1:200,,], outer=~as.factor(group[1:200]), aspect=2)
p1 # screen plot
ps.options(horizontal=FALSE,onefile=FALSE,print.it=FALSE) # p.75 in VR
library(lattice)
trellis.device("postscript",file="c:\\data.avc\\teaching\\vhm882\\06l\\d_eg6_3.eps")
p1 # file plot
dev.off()

# mean plot based on 10-year intervals
hour.group <- floor(jb$age/10)
jb.means <- aggregate(jb$src,list(hour.group,jb$group),mean,na.rm=TRUE)
names(jb.means)[1]<-"agegrp"
jb.means$agegrp <- as.numeric(levels(jb.means$agegrp))[jb.means$agegrp]
names(jb.means)[2]<-"group"
library(lattice)
xyplot(x~agegrp|group,jb.means,ylab="Mean Reciprocal serum creatinine level",aspect=2)
# plot not very useful

# model in Table 6.7 and model 4 in Figures 6.7-8
jb.m4 <- lme(src~group*age,random=~age|subject,correlation=corExp(form=~age|subject,nugget=TRUE),data=jb,method="REML")
summary(jb.m4)

# model 1 in Figures 6.7-8
jb.m1 <- gls(src~group*age,correlation=corExp(form=~age|subject,nugget=TRUE),data=jb,method="REML")
summary(jb.m1)

# model 2 in Figures 6.7-8
jb.m2 <- update(jb.m1,correlation=corExp(form=~age|subject,nugget=FALSE))
summary(jb.m2)
# jb.m2 <- update(jb.m1,correlation=corCAR1(form=~age|subject))
# same result: note that range=-1/ln(phi)

# model 3 in Figures 6.7-8
jb.m3 <- update(jb.m1,correlation=corCompSymm(form=~1|subject))
summary(jb.m3)

# model 5 in Figures 6.7-8
jb.m5 <- update(jb.m4,correlation=corExp(form=~age|subject,nugget=FALSE))
summary(jb.m5)
# jb.m5 <- update(jb.m4,correlation=corCAR1(form=~age|subject))
# same result: note that range=-1/ln(phi)

# model 6 in Figures 6.7-8
jb.m6 <- update(jb.m4,correlation=corIdent(form=~1))
summary(jb.m6)

# model 7 in Figures 6.7-8
jb.m7 <- update(jb.m4,random=~1|subject)
summary(jb.m7)

# model 8 in Figures 6.7-8
jb.m8 <- update(jb.m4,random=~1|subject,correlation=corExp(form=~age|subject,nugget=FALSE))
summary(jb.m8)
# jb.m8 <- update(jb.m4,random=~1|subject,correlation=corCAR1(form=~age|subject))
# same result: note that range=-1/ln(phi)

# model 9 in Figures 6.7-8
jb.m9 <- update(jb.m4,random=~1|subject,correlation=corIdent(form=~1))
summary(jb.m9)