dryeff <- read.table("http://www.stat.ufl.edu/~winner/data/dryer_efficiency.dat", header=F,col.names=c("clothcat","dryertype","energy.eff","cycle.time")) attach(dryeff) X <- cbind(energy.eff, cycle.time) clothcat <- factor(clothcat) dryertype <- factor(dryertype) eff.mod <- manova(X ~ clothcat*dryertype) summary(eff.mod, test="Wilks") ### Direct Computations (N <- nrow(X)) (p <- ncol(X)) (g <- length(unique(clothcat))) (b <- length(unique(dryertype))) (gb <- g*b) n.grp <- rep(0,gb) n.A <- rep(0,g) n.B <- rep(0,b) trt.hold <- 0 for (i1 in 1:g) { for (i2 in 1:b) { trt.hold <- trt.hold + 1 n.grp[trt.hold] <- length(energy.eff[(clothcat==i1)&(dryertype==i2)]) n.A[i1] <- n.A[i1] + n.grp[trt.hold] n.B[i2] <- n.B[i2] + n.grp[trt.hold] } } n.grp n.A n.B (xbar.all <- (1/N) * t(X) %*% rep(1,N)) X.AB <- matrix(rep(0,N*gb), ncol=gb) cume.grp <- 0 for (i1 in 1:gb) { for (i2 in 1:n.grp[i1]) { X.AB[cume.grp + i2, i1] <- 1 } cume.grp <- cume.grp + n.grp[i1] } X.AB (xbar.AB <- (1/n.grp) * t(X) %*% X.AB) X.A <- matrix(rep(0,N*g), ncol=g) cume.A <- 0 cume.grp <- 0 for (i1 in 1:g) { for (i2 in 1:b) { cume.grp <- cume.grp + 1 for (i3 in 1:n.grp[cume.grp]) { X.A[cume.A+i3,i1] <- 1 } cume.A <- cume.A + n.grp[cume.grp] } } X.A (xbar.A <- (1/n.A) * t(X) %*% X.A) X.B <- matrix(rep(0,N*b), ncol=b) cume.B <- 0 cume.grp <- 0 for (i1 in 1:g) { for (i2 in 1:b) { cume.grp <- cume.grp + 1 for (i3 in 1:n.grp[cume.grp]) { X.B[cume.B+i3,i2] <- 1 } cume.B <- cume.B + n.grp[cume.grp] } } X.B (xbar.B <- (1/n.B) * t(X) %*% X.B) Xbar.mat.AB <- matrix(rep(0,N*p), ncol=p) Xbar.mat.A <- matrix(rep(0,N*p), ncol=p) Xbar.mat.B <- matrix(rep(0,N*p), ncol=p) Xbar.mat.all <- matrix(rep(0,N*p), ncol=p) for (i0 in 1:p) { cume.grp <- 0 num.grp <- 0 for (i1 in 1:g) { for (i2 in 1:b) { num.grp <- num.grp + 1 for (i3 in 1:n.grp[num.grp]) { Xbar.mat.AB[cume.grp + i3, i0] <- xbar.AB[i0,num.grp] Xbar.mat.A[cume.grp + i3, i0] <- xbar.A[i0,i1] Xbar.mat.B[cume.grp + i3, i0] <- xbar.B[i0,i2] Xbar.mat.all[cume.grp + i3, i0] <- xbar.all[i0,1] } cume.grp <- cume.grp + n.grp[num.grp] } } } Xbar.mat.AB Xbar.mat.A Xbar.mat.B Xbar.mat.all (SSCP.ERR <- t(X - Xbar.mat.AB) %*% (X - Xbar.mat.AB)) (SSCP.AB <- t(Xbar.mat.AB-Xbar.mat.A-Xbar.mat.B+Xbar.mat.all) %*% (Xbar.mat.AB-Xbar.mat.A-Xbar.mat.B+Xbar.mat.all)) (SSCP.A <- t(Xbar.mat.A-Xbar.mat.all) %*% (Xbar.mat.A-Xbar.mat.all)) (SSCP.B <- t(Xbar.mat.B-Xbar.mat.all) %*% (Xbar.mat.B-Xbar.mat.all)) ## Wilks' Lambda -- F-test -- AB Interaction Lambda.AB <- det(SSCP.ERR)/det(SSCP.ERR + SSCP.AB) q.AB <- (g-1) * (b-1) r.AB <- (N-gb) - (p-q.AB+1)/2 u.AB <- (p*q.AB-2)/4 t.AB <- ifelse(p^2+q.AB^2-5>0,sqrt((p^2*q.AB^2-4)/(p^2+q.AB^2-5)),1) F.AB.W <- ((1-Lambda.AB^(1/t.AB))/(Lambda.AB^(1/t.AB))) * ((r.AB*t.AB - 2*u.AB)/(p*q.AB)) df1.AB.W <- p*q.AB; df2.AB.W <- r.AB*t.AB - 2*u.AB CV.AB.W <- qf(.95, df1.AB.W, df2.AB.W) PV.AB.W <- 1 - pf(F.AB.W, df1.AB.W, df2.AB.W) W.AB.out <- cbind(Lambda.AB, F.AB.W, df1.AB.W, df2.AB.W, CV.AB.W, PV.AB.W) colnames(W.AB.out) <- c("Lambda AB", "F", "df1", "df2", "F(.05)", "P-value") round(W.AB.out, 4) ## Chi-square Test -- AB Interaction n <- n.grp[1] ## Must be balanced sample sizes X2.AB <- -(g*b*(n-1) - (p+1-(g-1)*(b-1))/2) * log(Lambda.AB) df.AB.X2 <- p*(g-1)*(b-1) CV.AB.X2 <- qchisq(.95, df.AB.X2) PV.AB.X2 <- 1 - pchisq(X2.AB, df.AB.X2) X2.AB.out <- cbind(Lambda.AB, X2.AB, df.AB.X2, CV.AB.X2, PV.AB.X2) colnames(X2.AB.out) <- c("Lambda AB", "X2-stat", "df", "X2(.05)", "P-value") round(X2.AB.out, 4) ## Wilks' Lambda -- F-test -- Factor A Lambda.A <- det(SSCP.ERR)/det(SSCP.ERR + SSCP.A) q.A <- g-1 r.A <- (N-gb) - (p-q.A+1)/2 u.A <- (p*q.A-2)/4 t.A <- ifelse(p^2+q.A^2-5>0,sqrt((p^2*q.A^2-4)/(p^2+q.A^2-5)),1) F.A.W <- ((1-Lambda.A^(1/t.A))/(Lambda.A^(1/t.A))) * ((r.A*t.A - 2*u.A)/(p*q.A)) df1.A.W <- p*q.A; df2.A.W <- r.A*t.A - 2*u.A CV.A.W <- qf(.95, df1.A.W, df2.A.W) PV.A.W <- 1 - pf(F.A.W, df1.A.W, df2.A.W) W.A.out <- cbind(Lambda.A, F.A.W, df1.A.W, df2.A.W, CV.A.W, PV.A.W) colnames(W.A.out) <- c("Lambda A", "F", "df1", "df2", "F(.05)", "P-value") round(W.A.out, 4) ## Chi-square Test -- Factor A n <- n.grp[1] ## Must be balanced sample sizes X2.A <- -(g*b*(n-1) - (p+1-(g-1))/2) * log(Lambda.A) df.A.X2 <- p*(g-1) CV.A.X2 <- qchisq(.95, df.A.X2) PV.A.X2 <- 1 - pchisq(X2.A, df.A.X2) X2.A.out <- cbind(Lambda.A, X2.A, df.A.X2, CV.A.X2, PV.A.X2) colnames(X2.AB.out) <- c("Lambda A", "X2-stat", "df", "X2(.05)", "P-value") round(X2.A.out, 4) ## Wilks' Lambda -- F-test -- Factor B Lambda.B <- det(SSCP.ERR)/det(SSCP.ERR + SSCP.B) q.B <- b-1 r.B <- (N-gb) - (p-q.B+1)/2 u.B <- (p*q.B-2)/4 t.B <- ifelse(p^2+q.B^2-5>0,sqrt((p^2*q.B^2-4)/(p^2+q.B^2-5)),1) F.B.W <- ((1-Lambda.B^(1/t.B))/(Lambda.B^(1/t.B))) * ((r.B*t.B - 2*u.B)/(p*q.B)) df1.B.W <- p*q.B; df2.B.W <- r.B*t.B - 2*u.B CV.B.W <- qf(.95, df1.B.W, df2.B.W) PV.B.W <- 1 - pf(F.B.W, df1.B.W, df2.B.W) W.B.out <- cbind(Lambda.B, F.B.W, df1.B.W, df2.B.W, CV.B.W, PV.B.W) colnames(W.B.out) <- c("Lambda B", "F", "df1", "df2", "F(.05)", "P-value") round(W.B.out, 4) ## Chi-square Test -- Factor B n <- n.grp[1] ## Must be balanced sample sizes X2.B <- -(g*b*(n-1) - (p+1-(b-1))/2) * log(Lambda.B) df.B.X2 <- p*(b-1) CV.B.X2 <- qchisq(.95, df.B.X2) PV.B.X2 <- 1 - pchisq(X2.B, df.B.X2) X2.B.out <- cbind(Lambda.B, X2.B, df.B.X2, CV.B.X2, PV.B.X2) colnames(X2.B.out) <- c("Lambda B", "X2-stat", "df", "X2(.05)", "P-value") round(X2.B.out, 4) #### Bonferroni Comparisons among Factor A m <- p*g*(g-1)/2 ### Number of Comparisons Bon.t <- qt(1-.05/(2*m),N-gb) Mean.21 <- matrix(rep(0,2*m),ncol=2) SE.diff <- rep(0,m) Bon.CI <- matrix(rep(0,m*3),ncol=3) Bon.MSD <- rep(0,m) Bon.count <- 0 Trt.Lab <- matrix(rep(0,2*m),ncol=2) for (i1 in 1:p) { SE.diff1 <- sqrt((SSCP.ERR[i1,i1]/(N-gb))*(2/(b*n))) for (i2 in 1:(g-1)) { for (i3 in (i2+1):g) { Bon.count <- Bon.count + 1 SE.diff[Bon.count] <- SE.diff1 Bon.MSD[Bon.count] <- Bon.t * SE.diff[Bon.count] Bon.CI[Bon.count,1] <- (xbar.A[i1,i3] - xbar.A[i1,i2]) Bon.CI[Bon.count,2] <- (xbar.A[i1,i3] - xbar.A[i1,i2]) - Bon.MSD[Bon.count] Bon.CI[Bon.count,3] <- (xbar.A[i1,i3] - xbar.A[i1,i2]) + Bon.MSD[Bon.count] Trt.Lab[Bon.count,1] <- i3 Trt.Lab[Bon.count,2] <- i2 Mean.21[Bon.count,1] <- xbar.A[i1,i3] Mean.21[Bon.count,2] <- xbar.A[i1,i2] } } } Bon.out <- cbind(Trt.Lab[,1],Trt.Lab[,2],Mean.21[,1], Mean.21[,2],Bon.CI[,1], SE.diff,Bon.MSD, Bon.CI[,2], Bon.CI[,3]) colnames(Bon.out) <- c("Trt 1", "Trt 2", "Mean 1", "Mean 2", "Diff", "Std Err", "MSD", "LB", "UB") round(Bon.out,3)