rm(list = ls())  # clear the memory
###########################

data = read.table("http://www.stat.ualberta.ca/~wiens/stat575/datasets/T6-13.dat")
X = as.matrix(data[,1:4])
time = data[,5]

n1 = sum(time==1)
n2 = sum(time==2)
n3 = sum(time==3)
n = n1 + n2 + n3
ens = c(n1,n2,n3)
p=4
g = 3

S1 = cov(X[time==1,])
S2 = cov(X[time==2,])
S3 = cov(X[time==3,])
W = (n1-1)*S1 + (n2-1)*S2 + (n3-1)*S3
Spooled = W/(n-g)

# Test equality of the covariances
u = (sum(1/(ens-1)) - 1/sum(ens-1))*(2*p*p + 3*p - 1)/(6*(p+1)*(g-1))
M = (n-g)*log(det(Spooled)) - 
	((n1-1)*log(det(S1)) + (n2-1)*log(det(S2)) + (n3-1)*log(det(S3)))
C = (1-u)*M
dfC = p*(p+1)*(g-1)/2
pval = pchisq(C, dfC, lower.tail = 0)
pval # .39


xbar = colMeans(X)
tauhat1 = colMeans(X[time==1,]) - xbar
tauhat2 = colMeans(X[time==2,]) - xbar
tauhat3 = colMeans(X[time==3,]) - xbar
tauhat = cbind(tauhat1, tauhat2, tauhat3)

#Test overall treatment effects using Bartlett's approximation:
B = tauhat%*%diag(ens)%*%t(tauhat)
Wilks = det(W)/det(B+W) #  0.8301027
bart = -(n-1-((p+g)/2))*log(Wilks)
asym.pval = pchisq(bart, p*(g-1), lower.tail = 0)
asym.pval # .0435

# 'matplot' plots all the columns of one matrix against all the columns of another.  
# Here we plot the 4 columns of tauhat against the one column c(1,2,3,4):

matplot(1:4, tauhat, type = 'l', xaxp = c(1,4,3), 
	lty = c(1,2,4), col = c(1,2,1), xlab = "component", ylab = "treatment effects")
legend("topright", c("tau1", "tau2", "tau3"), lty = c(1,2,4), col = c(1,2,1))


fit = manova(X ~ as.factor(time)) #Note: time is a vector of labels - 'factors'
summary(fit)
summary(fit, test = "Wilks") # p = .044
summary(fit, test = "Roy") # p = .005
summary(fit, test = "Hotelling-Lawley") # p = .039
summary.aov(fit)

#Bonferroni intervals on all differences
alpha = .25 # Some are significant at alpha = .25
m = p*g*(g-1)/2
q = qt(alpha/(2*m), n-g, lower.tail = 0)
confint = array(dim=c(g,g,p,3),
	dimnames = list(NULL, NULL, NULL, c("lower", "point", "upper")))
for(k in 1:g) {
	for(l in 1:g) {
		halfwidth = q*sqrt((diag(W)/(n-g))*(1/ens[k] + 1/ens[l]))
		mid = tauhat[,k] - tauhat[,l]
		lower = mid - halfwidth
		upper = mid + halfwidth
		confint[k,l, , ] = cbind(lower,mid,upper)
	}
}

for(k in 2:g) { for(l in 1:(k-1)) {
mat = round(confint[k,l, ,],3)
sig = sign(mat[,1]*mat[,3])
sig[sig==-1] = ' '
sig[sig==1] = '  *'
mat = cbind(mat, sig)
cat(paste(100*(1-alpha),"% Bonferroni confidence intervals on tau", k, "minus tau", l, "are:","\n"))
prmatrix(mat, quote = 0)
cat("\n")
}}