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

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]

# test for parallel curves (already in Example 3.3)
diff <- matrix(c(y[,1]-y[,2],y[,2]-y[,3],y[,3]-y[,4]),nrow=27,ncol=3,byrow=FALSE)
summary(manova(diff~sex))

# test for equality between boys and girls, without assuming parallelism (already in Example 3.3)
summary(manova(y~sex))

# test for equality between boys and girls, assuming parallelism
sum <- y[,1]+y[,2]+y[,3]+y[,4]
t.test(sum~sex)
t.test(sum~sex,var.equal=TRUE)

# test for no time differences, not assuming parallelism
diff <- matrix(c(y[,1]-y[,2],y[,2]-y[,3],y[,3]-y[,4]),nrow=27,ncol=3,byrow=FALSE)
summary(manova(diff~as.factor(sex)-1),test="Wilk")

# test for no time differences for boys and girls separately
girl<- sex-1
boy<- 2-sex
summary(manova(diff~boy+girl-1))

# test for time differences, assuming parallelism, unequal weights
diff <- matrix(c(y[,1]-y[,2],y[,2]-y[,3],y[,3]-y[,4]),nrow=27,ncol=3,byrow=FALSE)
summary(manova(diff~sex),intercept=TRUE)

# test for no time differences, assuming parallelism, equal weights
intercept.corr <- 1/16*boy+1/11*girl
summary(manova(diff~intercept.corr+sex-1))


# Appendix: multivariate analysis using matrix formulae

# multivariate analysis for time differences
x <- matrix(c(2-sex,sex-1),nrow=27,ncol=2,byrow=FALSE)
library(MASS)
bhat <- ginv(t(x) %*% x) %*% crossprod(x,y) # t(x) %*% y

# test for no time differences, assuming parallelism, equal weights
a <- matrix(c(1,1),nrow=1,ncol=2)
c <- matrix(c(diag(3),-rep(1,3)),nrow=4,ncol=3,byrow=TRUE)
d <- matrix(0,nrow=1,ncol=3)
ssp.hyp <- t(a %*% bhat %*% c - d) %*% ginv(a %*% ginv(crossprod(x,x)) %*% t(a)) %*% (a %*% bhat %*% c - d)
ssp.res <- t(c) %*% (crossprod(y,y) - t(bhat) %*% crossprod(x,x) %*% bhat) %*% c
wilk.lambda <- det(ssp.res) / det(ssp.hyp+ssp.res)
# wilk.lambda ~ wilk (4-1,27-2,1)=(3,25,1)
f <- (1-wilk.lambda)/wilk.lambda * (25+1-3)/3
f
pf(f,3,23)

# test for time differences, assuming parallelism, unequal weights
a <- matrix(c(16/27,11/27),nrow=1,ncol=2)
c <- matrix(c(diag(3),-rep(1,3)),nrow=4,ncol=3,byrow=TRUE)
d <- matrix(0,nrow=1,ncol=3)
ssp.hyp <- t(a %*% bhat %*% c - d) %*% ginv(a %*% ginv(crossprod(x,x)) %*% t(a)) %*% (a %*% bhat %*% c - d)
ssp.res <- t(c) %*% (crossprod(y,y) - t(bhat) %*% crossprod(x,x) %*% bhat) %*% c
wilk.lambda <- det(ssp.res) / det(ssp.hyp+ssp.res)
# wilk.lambda ~ wilk (4-1,27-2,1)=(3,25,1)
f <- (1-wilk.lambda)/wilk.lambda * (25+1-3)/3
f
pf(f,3,23)
# same result from analysis of differences (test for intercept)
diff <- matrix(c(y[,1]-y[,2],y[,2]-y[,3],y[,3]-y[,4]),nrow=27,ncol=3,byrow=FALSE)
summary(manova(diff~sex),intercept=TRUE)

# test for no time differences, not assuming parallelism
a <- diag(2)
c <- matrix(c(diag(3),-rep(1,3)),nrow=4,ncol=3,byrow=TRUE)
d <- matrix(0,nrow=2,ncol=3)
ssp.hyp <- t(a %*% bhat %*% c - d) %*% ginv(a %*% ginv(crossprod(x,x)) %*% t(a)) %*% (a %*% bhat %*% c - d)
ssp.res <- t(c) %*% (crossprod(y,y) - t(bhat) %*% crossprod(x,x) %*% bhat) %*% c
wilk.lambda <- det(ssp.res) / det(ssp.hyp+ssp.res)
# wilk.lambda ~ wilk (4-1,27-2,2)=(3,25,2)
f <- (1-sqrt(wilk.lambda))/sqrt(wilk.lambda) * (25-3+1)/3
f
pf(f,2*3,2*(25-3+1))

# test for no time differences for boys
a <- matrix(c(1,0),nrow=1,ncol=2)
c <- matrix(c(diag(3),-rep(1,3)),nrow=4,ncol=3,byrow=TRUE)
d <- matrix(0,nrow=1,ncol=3)
ssp.hyp <- t(a %*% bhat %*% c - d) %*% ginv(a %*% ginv(crossprod(x,x)) %*% t(a)) %*% (a %*% bhat %*% c - d)
ssp.res <- t(c) %*% (crossprod(y,y) - t(bhat) %*% crossprod(x,x) %*% bhat) %*% c
wilk.lambda <- det(ssp.res) / det(ssp.hyp+ssp.res)
# wilk.lambda ~ wilk (4-1,27-2,1)=(3,25,1)
f <- (1-wilk.lambda)/wilk.lambda * (25+1-3)/3
f
pf(f,3,23)

# test for no time differences for girls
a <- matrix(c(0,1),nrow=1,ncol=2)
c <- matrix(c(diag(3),-rep(1,3)),nrow=4,ncol=3,byrow=TRUE)
d <- matrix(0,nrow=1,ncol=3)
ssp.hyp <- t(a %*% bhat %*% c - d) %*% ginv(a %*% ginv(crossprod(x,x)) %*% t(a)) %*% (a %*% bhat %*% c - d)
ssp.res <- t(c) %*% (crossprod(y,y) - t(bhat) %*% crossprod(x,x) %*% bhat) %*% c
wilk.lambda <- det(ssp.res) / det(ssp.hyp+ssp.res)
# wilk.lambda ~ wilk (4-1,27-2,1)=(3,25,1)
f <- (1-wilk.lambda)/wilk.lambda * (25+1-3)/3
f
pf(f,3,23)