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