# Exercise 2.13 in Davis, 2002

week <- 4*c(0:10)
#m <- matrix(scan("e:\\vhm\\vhm882\\data\\leppik.dat"),ncol=6,byrow=T)
m <- matrix(scan("d:\\data.ext\\davis\\amato.dat",what="character",na.strings="."),ncol=14,byrow=T)
subject <- type.convert(m[,1])
sex <- type.convert(m[,2])
group <- type.convert(m[,3])
ybase <- type.convert(m[,4])
y <- type.convert(m[,5:14])
yfull <- type.convert(m[,4:14])

# code for graphics, using nlme and lattice libraries
library(nlme)
amato <- balancedGrouped( y ~ week|subject, matrix(type.convert(m[,4:14]),nrow=227,ncol=11,dimnames=list(subject,week)),
                           labels=list(y="Time (sec) to walk 50 feet"))
Group <- rep(group,rep(11,227))
Sex <- rep(sex,rep(11,227))
# profile plots
plot(amato[1:200,,], outer=~Group[1:200], aspect=1)
plot(amato[1:200,,], outer=~Sex[1:200], aspect=1)
plot(amato[1:200,,], outer=~Group[1:200]+Sex[1:200], aspect=1)
plot(amato[1:200,,], aspect=1)
# mean plot
amato.means <- aggregate(amato$y,list(amato$week,Sex,Group),mean,na.rm=TRUE)
names(amato.means)[1]<-"Week"
amato.means$Week <- as.numeric(levels(amato.means$Week))[amato.means$Week]
names(amato.means)[2]<-"Sex"
names(amato.means)[3]<-"Group"
library(lattice)
xyplot(x~Week|Sex*Group,amato.means,ylab="Mean Time (sec) to walk 50 feet")

# summary statistic: slope
slope <- c()
for (i in 1:227) 
  slope=c(slope,coef(lm(yfull[i,]~week))[2])
names(slope) <- NULL
aggregate(slope,list(group),mean,na.rm=TRUE)
aggregate(slope,list(group),sd,na.rm=TRUE)
aggregate(slope,list(group),median,na.rm=TRUE)
t.test(slope ~ group)
wilcox.test(slope ~ group)

# summary statistic: ymean-ybase
ymean <- rowMeans(y,na.rm=TRUE)
drop <- ybase-ymean
aggregate(drop,list(group),mean,na.rm=TRUE)
aggregate(drop,list(group),sd,na.rm=TRUE)
aggregate(drop,list(group),median,na.rm=TRUE)
t.test(drop ~ group)
wilcox.test(drop ~ group)

# separate analyses by sex
# summary statistic: slope
aggregate(slope,list(group,sex),mean,na.rm=TRUE)
aggregate(slope,list(group,sex),sd,na.rm=TRUE)
aggregate(slope,list(group,sex),median,na.rm=TRUE)
t.test(slope ~ group,subset=(sex=="F"))
t.test(slope ~ group,subset=(sex=="M"))
wilcox.test(slope ~ group,subset=(sex=="F"))
wilcox.test(slope ~ group,subset=(sex=="M"))

# summary statistic: ymean-ybase
aggregate(drop,list(group,sex),mean,na.rm=TRUE)
aggregate(drop,list(group,sex),sd,na.rm=TRUE)
aggregate(drop,list(group,sex),median,na.rm=TRUE)
t.test(drop ~ group,subset=(sex=="F"))
t.test(drop ~ group,subset=(sex=="M"))
wilcox.test(drop ~ group,subset=(sex=="F"))
wilcox.test(drop ~ group,subset=(sex=="M"))