## Problem 1:

testdata <- read.table("http://www.stat.sc.edu/~hitchcock/testscoredata.txt", header=T)
attach(testdata)
testdata.noIDs <- testdata[,-1]  #to remove the ID numbers

summary(manova(cbind(math, reading ) ~   sex) ,test="Hotelling")

tapply( math,sex,mean)
#     boy     girl 
#84.09833 82.96250 

 tapply( reading,sex,mean)
#     boy     girl 
#81.78667 82.28000 


# Checking assumptions:
x1 <- testdata.noIDs[sex=='boy',-3]
x2 <- testdata.noIDs[sex=='girl',-3]

 chisplot(as.matrix(x1))
chisplot(as.matrix(x2))
 library(asbio)  # load the 'asbio' package once it is installed
Kullback(testdata.noIDs[,-3],sex)

 library(biotools)  # load the 'biotools' package once it is installed

boxM(testdata.noIDs[,-3],sex)


### PROBLEM 2:

hsb <- read.table("http://www.stat.sc.edu/~hitchcock/hsbdata.txt", header=T)
attach(hsb)

hsb.prob4 <- hsb[,c(5,8,9,10,11,12)]

hsb.numeric <- hsb[,c(8,9,10,11,12)]

my.manova <- manova(cbind(read, write, math, science, socst) ~ ses)

# Doing the MANOVA test using Wilks' Lambda:

summary(my.manova, test="Wilks")

###############
###
chisplot <- function(x) {
    if (!is.matrix(x)) stop("x is not a matrix")
    ### determine dimensions
    n <- nrow(x)
    p <- ncol(x)
    xbar <- apply(x, 2, mean)
    S <- var(x)
    S <- solve(S)
    index <- (1:n)/(n+1)
    xcent <- t(t(x) - xbar)
    di <- apply(xcent, 1, function(x,S) x %*% S %*% x,S)
    quant <- qchisq(index,p)
    plot(quant, sort(di), ylab = "Ordered distances",
         xlab = "Chi-square quantile", lwd=2,pch=1)
}
###
###############

chisplot(residuals(my.manova))


# Examining the sample covariance matrices for each group:

by(cbind(read, write, math, science, socst), ses, var)

# Optionally, doing the formal tests:

library(asbio)

Kullback(Y=cbind(read, write, math, science, socst),X=ses)

library(biotools)  
boxM(cbind(read, write, math, science, socst),ses)

##############################################

# Examining the sample mean vectors for each group:

by(cbind(read, write, math, science, socst), ses, colMeans)