# Example 2.3 (Table 2.5) in Davis, 2002

day <- c(1,8,15,22,29,36,43,44,50,57,64)
m <- matrix(scan("e:\\vhm\\vhm882\\data\\ch4.dat"),ncol=13,byrow=T)
group <- m[,1]
subject <- m[,2]
y <- m[,3:13]
slope <- c()
for (i in 1:16)
  slope=c(slope,coef(lm(y[i,]~day))[2])
names(slope) <- NULL
anova(lm(slope~as.factor(group))) # one-way ANOVA for slopes
kruskal.test(slope,group) # Kruskal-Wallis test for slopes

# code for graphics, using nlme library
library(nlme)
ch4 <- balancedGrouped( y ~ day|subject, matrix(m[,3:13],nrow=16,ncol=11,dimnames=list(subject,day)),
                        labels=list(y="Rat body weight"))
group <- rep(group,rep(11,16))
p1 <- plot(ch4, outer=~as.factor(group), 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\\01l\\d_eg2_3.eps")
p1 # file plot
dev.off()