# Example 6.2 in Davis, 2002

hour <- c(1:24)
m <- matrix(scan("e:\\vhm\\vhm882\\data\\davis4.dat",what="character",na.strings="."),ncol=26,byrow=T)
#m <- matrix(scan("d:\\data.ext\\davis\\davis4.dat",what="character",na.strings="."),ncol=26,byrow=T)
group <- as.factor(m[,1])
subject <- type.convert(m[,2])
y <- type.convert(m[,3:26])

# code for graphics, using nlme and lattice libraries
library(nlme)
rescue <- balancedGrouped( y ~ hour|subject, matrix(type.convert(m[,3:26]),nrow=122,ncol=24,dimnames=list(subject,hour)),
                           labels=list(y="Amount of rescue medication used"))
Group <- rep(group,rep(24,122))
# profile plot (5 subjects from each group)
index<-c(1:120,1105:1224,2209:2328)
p1 <- plot(rescue[index,,], outer=~as.factor(Group[index]), 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_2.eps")
p1 # file plot
dev.off()

# mean plot
rescue.means <- aggregate(rescue$y,list(rescue$hour,Group),mean,na.rm=TRUE)
names(rescue.means)[1]<-"hour"
rescue.means$hour <- as.numeric(levels(rescue.means$hour))[rescue.means$hour]
names(rescue.means)[2]<-"Group"
library(lattice)
xyplot(x~hour|Group,rescue.means,ylab="Mean Amount of rescue medication used",aspect=2)

# simple AR(1) structure
rescue.ar1 <- gls(y~Group+I(hour-12)+I((hour-12)^2),correlation=corAR1(form=~1|subject),data=rescue,method="REML")
summary(rescue.ar1)

# ar(1) structure + random intercepts
rescue.ar1re <- lme(y~Group+I(hour-12)+I((hour-12)^2),random=~1|subject,correlation=corAR1(form=~1|subject),data=rescue,method="REML")
summary(rescue.ar1re)

# random slopes model, no AR(1) structure
rescue.rslp <- lme(y~Group+I(hour-12)+I((hour-12)^2),random=~I(hour-12)|subject,data=rescue,method="REML")
summary(rescue.rslp)

# extra models
# arma(1,1) errors
rescue.arma11 <- update(rescue.ar1,correlation=corARMA(c(0.2,0.2),form=~1|subject,p=1,q=1))
summary(rescue.arma11)

# arma(1,1) errors + random intercept
rescue.arma11re <- update(rescue.ar1re,correlation=corARMA(c(0.2,0.2),form=~1|subject,p=1,q=1))
summary(rescue.arma11re)

# ar(1) structure + random intercepts + heterogeneous variances
rescue.ar1rehet <- update(rescue.ar1re,weights=varIdent(form=~1|hour))
summary(rescue.ar1rehet)

# ar(1) structure + random slopes
rescue.rslp.ar1 <- update(rescue.ar1re,random=~hour|subject)
summary(rescue.rslp.ar1)

# arma(1,1) structure + random slopes
rescue.rslp.arma11 <- update(rescue.rslpar1,correlation=corARMA(c(0.2,0.2),form=~1|subject,p=1,q=1))
summary(rescue.rslp.arma11)