# Exercise 2.4 in Davis, 2002

week <- c(0:4)
#m <- matrix(scan("e:\\vhm\\vhm882\\data\\leppik.dat"),ncol=6,byrow=T)
m <- matrix(scan("c:\\data.ext\\davis\\box.dat"),ncol=7,byrow=T)
group <- m[,1]
rat <- m[,2]
y <- m[,3:7]

# code for graphics, using nlme and lattice libraries
library(nlme)
box <- balancedGrouped( y ~ week|rat, matrix(m[,3:7],nrow=27,ncol=5,dimnames=list(rat,week)),
                           labels=list(y="Weight of rats in 4 weeks"))
Group <- rep(group,rep(5,27))
# profile plot
plot(box, outer=~as.factor(Group), aspect=2)
# mean plot
box.means <- aggregate(box$y,list(box$week,Group),mean,na.rm=TRUE)
names(box.means)[1]<-"Week"
box.means$Week <- as.numeric(levels(box.means$Week))[box.means$Week]
names(box.means)[2]<-"Group"
library(lattice)
xyplot(x~Week|Group,box.means,ylab="Mean Weight of rats in 4 weeks",aspect=2)

# summary statistic: slope
slope <- c()
for (i in 1:27) 
  slope=c(slope,coef(lm(y[i,]~week))[2])
names(slope) <- NULL
aggregate(slope,list(group),mean)
aggregate(slope,list(group),sd)
aggregate(slope,list(group),median)
anova(lm(slope~as.factor(group))) # one-way ANOVA for slopes
shapiro.test(residuals(lm(slope~as.factor(group))))
pairwise.t.test(slope,group,p.adj="bonf")
kruskal.test(slope,group) # Kruskal-Wallis test for slopes
pairwise.wilcox.test(slope,group,p.adj="bonf")