# Example 4.2 in Davis, 2002
# Continuation of Example 3.2

age <- c(8,8.5,9,9.5)
m <- matrix(scan("e:\\vhm\\vhm882\\data\\eg.dat"),ncol=5,byrow=T)
boy <- m[,1]
y <- m[,2:5]

# orthogonal T
orth <- poly(age,3)
t <- matrix(c(rep(0.5,4),orth),nrow=4,ncol=4,byrow=TRUE)
z <- y %*% t(t)
summary.aov(manova(z~1),intercept=TRUE)
t.test(z[,1])
t.test(z[,2])
t.test(z[,3])
t.test(z[,4])
# test for linear and quadratic terms combined in Example 3.2

library(MASS)
# linear T
t <- matrix(c(rep(1,4),age),nrow=2,ncol=4,byrow=TRUE)
# G=I(4)
g <- diag(4)
g.inv <- diag(4)
z <- y %*% g.inv %*% t(t) %*% ginv(t %*% g.inv %*% t(t))
summary.aov(manova(z~1),intercept=TRUE)
t.test(z[,1])
t.test(z[,2])
# G=S
x <- matrix(rep(1,20),nrow=20,ncol=1)
bhat <- ginv(t(x) %*% x) %*% crossprod(x,y) # t(x) %*% y
s <- crossprod(y - x %*% bhat)/(20-1)
s.inv <- ginv(s)
z <- y %*% s.inv %*% t(t) %*% ginv(t %*% s.inv %*% t(t))
summary.aov(manova(z~1),intercept=TRUE)
t.test(z[,1])
t.test(z[,2])