## Multivariate approach to repeated measures

########################################
# data: from table 3.3 of Davis, page 58
########################################
# A study was conducted at the U of N Carolina Dental School in two groups of
# children (16 boys and 11 girls). At ages 8, 10, 12, and 14, the distance (mm) 
# from the center of the pituitary gland to the pterygomaxillary fissure was measured. 
y8<-c(26,21.5,23,25.5,20,24.5,22,24,23,27.5,23,21.5,17.0,22.5,23,22,
   21,21,20.5,23.5,21.5,20,21.5,23,20,16.5,24.5)
y10<-c(25,22.5,22.5,27.5,23.5,25.5,22,21.5,20.5,28,23,23.5,24.5,25.5,24.5,21.5,
   20,21.5,24,24.5,23,21,22.5,23,21,19,25)
y12<-c(29,23,24,26.5,22.5,27,24.5,24.5,31,31,23.5,24,26,25.5,26,23.5,
   21.5,24,24.5,25,22.5,21,23,23.5,22,19,28)
y14<-c(31,26.5,27.5,27,26,28.5,26.5,25.5,26,31.5,25,28,29.5,26,30,25,
   23,25.5,26,26.5,23.5,22.5,25,24,21.5,19.5,28)
y<-cbind(y8,y10,y12,y14)
group<-rep(1:2,c(16,11))

> y
        y8  y10  y12  y14
 [1,] 26.0 25.0 29.0 31.0
 [2,] 21.5 22.5 23.0 26.5
 [3,] 23.0 22.5 24.0 27.5
 ....    ...    ...  ...
[27,] 24.5 25.0 28.0 28.0
> group
 [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2

######################
# exploratory analysis
######################
#long.plot: This creats the subject-specific trajectory plots over time.
# This function was written by Norm Breslow. 
# clus: the identity of the subject
# time: the time the observation was made
# x: the observation

long.plot<-function(clus,time,x,xlabel="Time",ylabel="Data",
  sub=" ",main=" ",cex=1)
{
ord<-sort.list(clus)
nn<-length(clus)
clus<-clus[ord]
time<-time[ord]
x<-x[ord]
n<-seq(nn)
cnj<-n[!duplicated(clus)]-1
m<-length(cnj)   # number of clusters
cnj<-c(cnj,nn)   # vector of cluster sizes
for (i in 1:m){
  tj<-time[(cnj[i]+1):cnj[i+1]]
  xj<-x[(cnj[i]+1):cnj[i+1]]
  sel<-!is.na(tj) & !is.na(xj)
  if(i==1) plot(tj[sel],xj[sel],xlim=c(min(time),max(time)),
     ylim=c(min(x[!is.na(x)]),max(x[!is.na(x)])),type='l',xlab=xlabel,
     ylab=ylabel,sub=sub,main=main,cex=cex)
  else lines(tj[sel],xj[sel],lty=i)
}
}

## boys (group 1)
id<-rep(1:16,rep(4,16))
age<-rep(c(8,10,12,14),16)
meas<-as.vector(t(y[1:16,]))
long.plot(id,age,meas,xlabel='Age (years)',ylabel='Distance (mm)',
  main='Dental measurements from 16 boys',cex=1.5)

## girls (group 2)
id<-rep(17:27,rep(4,11))
age<-rep(c(8,10,12,14),11)
meas<-as.vector(t(y[17:27,]))
long.plot(id,age,meas,xlabel='Age (years)',ylabel='Distance (mm)',
  main='Dental measurements from 11 girls',cex=1.5)

# means
plot(c(8,10,12,14),apply(y[1:16,],2,mean),ylim=c(15,30),xlab='Age (years)',
  ylab='Mean distance (mm)',type='b',cex=1.5,
  main='Sample means of dental measurements from 16 boys and 11 girls')
lines(c(8,10,12,14),apply(y[17:27,],2,mean),lty=2)
points(c(8,10,12,14),apply(y[17:27,],2,mean),cex=1.5)
legend(9,18,c('Boys','Girls'),lty=1:2)

###########################
# test for treatment effect
###########################
y1bar<-apply(y[1:16,],2,mean)
y2bar<-apply(y[17:27,],2,mean)
S1<-cov(y[1:16,])
#S1<-matrix(rep(0,16),ncol=4)
#for (i in 1:16){
#S1<-S1+(y[i,]-y1bar) %*% t(y[i,]-y1bar)
#}
#S1<-S1/(16-1)

S2<-cov(y[17:27,])
S<-((16-1)*S1+(11-1)*S2)/(27-2)
T2<-(16*11)/27 * t(y1bar-y2bar) %*% solve(S) %*% (y1bar-y2bar)
# T2=16.51

Ftrt<-(16+11-4-1)/((16+11-2)*4) *T2  
# Ftrt=3.63 with df1=4, df2=22
> 1-pf(3.63,4,22)
[1] 0.0203736              # The mean vectors for boys and girls are diff.

### using manova
fit1<-manova(y~group)
summary(fit1, test="Hotelling")
          Df Hotelling-Lawley approx F num Df den Df  Pr(>F)  
group      1           0.6603   3.6317      4     22 0.02034   #sig grp effect
Residuals 25             

########################################
# test for treatment by time interaction
########################################
C<-matrix(c(-1,1,0,0,0,-1,1,0,0,0,-1,1),byrow=T,ncol=4)
> C
     [,1] [,2] [,3] [,4]
[1,]   -1    1    0    0
[2,]    0   -1    1    0
[3,]    0    0   -1    1
SC<-C %*% S %*% t(C)   
diffyC<-C %*% (y1bar-y2bar)
T2a<-(16*11)/27 * t(diffyC) %*% solve(SC) %*% diffyC
# T2a=8.79
Fa<-(27-3-1)/(25*3)*T2a
# Fa=2.70 with df1=3, df2=23
> 1-pf(2.70,3,23)
[1] 0.06927386  # not enough evidence of a time by group interaction

#####################
# test for time trend
#####################
# since no sig. time by group interaction was found, 
# we will use the combined sample to test a time trend

newy1<-y8-y10
newy2<-y10-y12
newy3<-y12-y14
newy<-cbind(newy1,newy2,newy3)
Sstar<-cov(newy)
ybarstar<-apply(newy,2,mean)
T2b<-27*t(ybarstar) %*% solve(Sstar) %*% ybarstar 
# T2b=93.78
Fb<-(27-4+1)/(26*3) *T2b
#Fb=28.86 with df1=3, df2=24
> 1-pf(28.86,3,24)
[1] 3.936405e-08
# There is a very strong time trend: the average distance changes over age.