# Exercise 2.3 in Davis, 2002

week <- c(0,2,4,6,8,10,12,14,16,18,19)
#m <- matrix(scan("e:\\vhm\\vhm882\\data\\leppik.dat"),ncol=6,byrow=T)
m <- matrix(scan("c:\\data.ext\\davis\\kenward.dat"),ncol=13,byrow=T)
group <- m[,1]
calf <- m[,2]
y <- m[,3:13]

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

# summary statistic: slope
slope <- c()
for (i in 1:60) 
  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)
t.test(slope ~ group)
wilcox.test(slope ~ group)

# summary statistic: curvature
curv <- c()
for (i in 1:60) 
  curv=c(curv,coef(lm(y[i,]~week+I(week^2)))[3])
names(curv) <- NULL
aggregate(curv,list(group),mean)
aggregate(curv,list(group),sd)
aggregate(curv,list(group),median)
t.test(curv ~ group)
wilcox.test(curv ~ group)