bac3 <- read.csv("http://www.stat.ufl.edu/~winner/data/bac_3methods_multi.csv") attach(bac3); names(bac3) (X <- cbind(BAC1,BAC2,BAC3)) (q <- ncol(X)) (n <- nrow(X)) (xbar <- (1/n) * (t(X) %*% rep(1,n))) I_n <- diag(n); J_n <- matrix(rep(1,n^2),n,n) (S <- (1/(n-1)) * (t(X) %*% (I_n - (1/n)*J_n) %*% X)) C <- matrix(c(1, -1, 0, 1, 0, -1),byrow=T,ncol=3) T2 <- n * (t(C%*%xbar) %*% solve(C %*% S %*% t(C)) %*% (C %*% xbar)) T2s <- ((n-(q-1))/((n-1)*(q-1))) * T2 CV <- qf(.95,q-1,n-q+1) PVal <- 1 - pf(T2s, q-1, n-q+1) T2.out <- cbind(T2, T2s, q-1, n-q+1, CV, PVal) colnames(T2.out) <- c("Hotelling's T2", "Scaled T2", "df1", "df2","Critical Value", "P-value") round(T2.out, 4) ## Test using HotellingsT2 function in ICSNP package # install.packages("ICSNP") library(ICSNP) C.mu1 <- BAC1 - BAC2 ## Contrast 1 (within subjects) C.mu2 <- BAC1 - BAC3 ## Contrast 2 (within subjects) muH0 <- c(0,0) ## mu0 NOT as column vector X.HT2 <- cbind(C.mu1, C.mu2) HotellingsT2(X.HT2,mu=muH0) ## 95% Confidence Region for Cm = [mu1-mu2 mu1-mu3]' alpha <- 0.05 crit.dist <- sqrt(((q-1)*(n-1)/(n*(n-(q-1))))*qf(1-alpha,q-1,n-(q-1))) A <- C %*% S %*% t(C) mu.test <- c(0.0,0.0) ctr <- c(mean(C.mu1), mean(C.mu2)) angles <- seq(0, 2*pi, length.out=200) eigVal <- eigen(A)$values eigVec <- eigen(A)$vectors # scale eigenvectors to length = square-root lambda eigScl <- eigVec %*% diag(sqrt(eigVal)) xMat <- rbind(ctr[1] + eigScl[1, ]*crit.dist, ctr[1] - eigScl[1, ]*crit.dist) yMat <- rbind(ctr[2] + eigScl[2, ]*crit.dist, ctr[2] - eigScl[2, ]*crit.dist) ellBase <- cbind(sqrt(eigVal[1])*crit.dist*cos(angles), sqrt(eigVal[2])*crit.dist*sin(angles)) # normal ellipse ellRot <- eigVec %*% t(ellBase) # rotated ellipse plot((ellRot+ctr)[1, ], (ellRot+ctr)[2, ], asp=1, type="l", lwd=2, main="100(1-a)% Confidence Ellipsoid", xlab="mu1-mu2", ylab="mu1-mu3") matlines(xMat, yMat, lty=1, lwd=2, col="blue") points(ctr[1], ctr[2], pch=4, col="orange", lwd=3) points(mu.test[1],mu.test[2], pch=8, col="purple") bon11 <- ctr[1] - qt(1-alpha/(2*(q-1)),n-1)*sqrt(A[1,1]/n) bon12 <- ctr[1] + qt(1-alpha/(2*(q-1)),n-1)*sqrt(A[1,1]/n) bon21 <- ctr[2] - qt(1-alpha/(2*(q-1)),n-1)*sqrt(A[2,2]/n) bon22 <- ctr[2] + qt(1-alpha/(2*(q-1)),n-1)*sqrt(A[2,2]/n) rect(bon11,bon21,bon12,bon22,border="violetred3")