# Example 4.1 (Table 3.3) in Davis, 2002; Orthodont data in Pinheiro & Bates
# Continuation of Examples 3.3 and 4.1

age <- c(8,10,12,14)
m <- matrix(scan("e:\\vhm\\vhm882\\data\\pr.dat"),ncol=6,byrow=T)
sex <- m[,1]
subject <- m[,2]
y <- m[,3:6]

# 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~as.factor(sex)-1))
t.test(z[,1]~sex)
t.test(z[,2]~sex)
t.test(z[,3]~sex)
t.test(z[,4]~sex)
# test of combined quadratic and cubic effects
summary(manova(z[,3:4]~as.factor(sex)-1),test="Wilk")

library(MASS)
# linear T (non-orthogonal)
t <- matrix(c(rep(1,4),age-8),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~sex),intercept=TRUE)
t.test(z[,1]~sex,var.equal=TRUE)
t.test(z[,2]~sex,var.equal=TRUE)

# G=S
x <- matrix(c(2-sex,sex-1),nrow=27,ncol=2,byrow=FALSE)
bhat <- ginv(t(x) %*% x) %*% crossprod(x,y) # t(x) %*% y
s <- crossprod(y - x %*% bhat)/(27-2)
s.inv <- ginv(s)
z <- y %*% s.inv %*% t(t) %*% ginv(t %*% s.inv %*% t(t))
summary.aov(manova(z~sex),intercept=TRUE)
summary(lm(z[,1]~as.factor(sex)-1))
summary(lm(z[,2]~as.factor(sex)-1))
t.test(z[,1]~sex,var.equal=TRUE)
t.test(z[,2]~sex,var.equal=TRUE)