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

# Profile analysis of spouses rating each other; data in Table 6.14

data = read.table("http://www.stat.ualberta.ca/~wiens/stat575/datasets/T6-14.dat")
X = as.matrix(data[,1:4])
spouse = c(rep(1,30), rep(2,30)) # 1 = Husband rating wife, 2 = Wife rating husband
spouse = as.factor(spouse)

xbar1 = colMeans(X[spouse==1,])
xbar2 = colMeans(X[spouse==2,])

# Plot of estimated profiles:
matplot(1:4, cbind(xbar1,xbar2), type = 'l', xaxp = c(1,4,3), 
	lty = 1:2, col= 1:2, xlab = "i", ylab = "average reponse to question 'i'", main = "Estimated profiles")
legend("topleft", c("husbands", "wives"), lty = 1:2, col= 1:2)

##################################
## Test for parallel profiles, two ways
C = rbind(c(-1, 1, 0, 0), c(0, -1, 1, 0), c(0, 0, -1, 1))

Y = X%*%t(C)

# Calc Tsqd:
ybar1 = colMeans(Y[spouse==1,])
ybar2 = colMeans(Y[spouse==2,])
n1 = nrow(Y[spouse==1,])
n2 = nrow(Y[spouse==2,])

dbar = ybar1 - ybar2
S1 = cov(Y[spouse==1,])
S2 = cov(Y[spouse==2,])
dfS = (n1-1) + (n2-1)
Spooled = ((n1-1)*S1 + (n2-1)*S2)/dfS
cov.est = (1/n1 + 1/n2)*Spooled
Tsqd = t(dbar)%*%solve(cov.est, dbar)
df1 = length(dbar)
df2 = dfS - df1 + 1
F = df2*Tsqd/(df1*dfS)
pval = pf(F, df1, df2, lower.tail = 0)
pval 	# .062559

# Using the R function manova():
fit = manova(Y ~ spouse)
summary(fit)

summary.aov(fit) # Third difference in slopes is most significant

##################################
## Test for coincident profiles

x1 = rowSums(X[spouse==1,])
x2 = rowSums(X[spouse==2,])
t.test(x1,x2, var.equal = T) # pval =  .2207
t.test(x1,x2) # Welch test; still pval =  .2207

# Using the R function aov():
fit = aov(rowSums(X) ~ spouse)
summary(fit)

##################################
## Test for level profiles, using the combined sample of Y's

dbar = colMeans(Y)
df1 = length(dbar)
S = cov(Y)
cov.est = S/nrow(Y)
dfS = nrow(Y) - 1
Tsqd = t(dbar)%*%solve(cov.est, dbar)
df2 = dfS - df1 + 1
Fcalc = df2*Tsqd/(df1*dfS)
pval = pf(Fcalc, df1, df2, lower.tail = 0)

pval 	# .00015


##################################
# Test for parallel profiles WITHOUT the assumption of equal covariance:
# Much of this is just copied and pasted from above

C = rbind(c(-1, 1, 0, 0), c(0, -1, 1, 0), c(0, 0, -1, 1))
Y = X%*%t(C)

# Calc Tsqd.tilde, start out exxactly as before:
ybar1 = colMeans(Y[spouse==1,])
ybar2 = colMeans(Y[spouse==2,])
n1 = nrow(Y[spouse==1,])
n2 = nrow(Y[spouse==2,])
dbar = ybar1 - ybar2
S1 = cov(Y[spouse==1,])
S2 = cov(Y[spouse==2,])

#Here things change. Compute a different cov.est and compute nu, then proceed as before:
cov.est = S1/n1 + S2/n2 #done formally; no change since n1=n2.
A1 = t(solve(cov.est,S1/n1))
A2 = t(solve(cov.est,S2/n2))
tr = function(M) sum(diag(M)) # define ones own "trace" function
denominator = (tr(A1%*%A1) + tr(A1)^2)/n1 + (tr(A2%*%A2) + tr(A2)^2)/n2
nu = (3 + 9)/denominator # p=3

# Now proceed exactly as before, with dfS = nu:
Tsqd.tilde = t(dbar)%*%solve(cov.est, dbar)
dfS = nu 		# 57.26954, previously 58
df1 = length(dbar) 
df2 = dfS - df1 + 1 # 55.26954, previously 56
F = df2*Tsqd.tilde/(df1*dfS)
pval = pf(F, df1, df2, lower.tail = 0)
pval 			# 0.06279063, previously .062559