pf1 <- read.csv("http://www.stat.ufl.edu/~winner/data/perfume_fakeorig.csv") attach(pf1); names(pf1) Cl.o <- ClMn[Original == 1] Na.o <- NaMn[Original == 1] K.o <- KMn[Original == 1] Ethnol.o <- EthnolMn[Original == 1] Cl.f <- ClMn[Original == 0] Na.f <- NaMn[Original == 0] K.f <- KMn[Original == 0] Ethnol.f <- EthnolMn[Original == 0] X.o <- cbind(Cl.o, Na.o, K.o, Ethnol.o) X.f <- cbind(Cl.f, Na.f, K.f, Ethnol.f) X.o <- as.matrix(X.o); X.f <- as.matrix(X.f) p <- ncol(X.o) n.o <- length(Cl.o); n.f <- length(Cl.f) I.o <- diag(n.o); I.f <- diag(n.f) J.o <- matrix(rep(1, n.o^2), n.o,n.o) J.f <- matrix(rep(1, n.f^2), n.f,n.f) one.o <- rep(1, n.o); one.f <- rep(1, n.f) (xbar.o <- (1/n.o) * t(X.o) %*% one.o) (S.o <- (1/(n.o-1)) * t(X.o) %*% (I.o - (1/n.o) * J.o) %*% X.o) (xbar.f <- (1/n.f) * t(X.f) %*% one.f) (S.f <- (1/(n.f-1)) * t(X.f) %*% (I.f - (1/n.f) * J.f) %*% X.f) (S.pooled <- ((n.o-1) * S.o + (n.f-1) * S.f) / (n.o + n.f -2)) delta.0 <- matrix(rep(0,length(xbar.o)), ncol=1) T2 <- t((xbar.o-xbar.f)-delta.0) %*% solve((1/n.o+1/n.f)*S.pooled) %*% ((xbar.o-xbar.f)-delta.0) T2s <- ((n.o+n.f-p-1)/(p*(n.o+n.f-2))) * T2 df1 <- p; df2 <- n.o+n.f-p-1 CV <- qf(.95,df1,df2) ## Critical Value (alpha=0.05) PVal <- 1-pf(T2s,df1,df2) ## P-value T2.out <- cbind(T2, T2s, df1, df2, 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) HotellingsT2(X.o, X.f) ### Simultaneous CI's for differences of mean characteristics # Critical Values fro Simultaneous and Bonferroni Intervals CV.sim <- sqrt(((n.o+n.f-2)*p/(n.o+n.f-p-1)) * qf(.95,df1,df2)) CV.bon <- qt(1-0.05/(2*p), n.o+n.f-2) diff <- xbar.o - xbar.f sediff <- sqrt(diag(S.pooled)*(1/n.o+1/n.f)) simCILo <- diff - sqrt(CV.sim) * sediff simCIHi <- diff + sqrt(CV.sim) * sediff bonCILo <- diff - CV.bon * sediff bonCIHi <- diff + CV.bon * sediff perfci <- cbind(diff,simCILo,simCIHi,bonCILo,bonCIHi) colnames(perfci) <- c("Diff","Sim Lo","Sim Hi", "Bon Lo", "Bon Hi") rownames(perfci) <- c("Cl", "Na", "K", "Ethnol") round(perfci, 3) ### Linear combination with maximum population differerence (ahat <- solve(S.pooled) %*% (xbar.o-xbar.f)) ### Case When Sigma1 ^= Sigma2 Tsq_un <- t((xbar.o-xbar.f)-delta.0) %*% solve((1/n.o)*S.o+(1/n.f)*S.f) %*% ((xbar.o-xbar.f)-delta.0) Tsq_un.CV <- qchisq(.95,p) Tsq_un.PV <- 1-pchisq(Tsq_un,p) Tsq_un.out <- cbind(Tsq_un, p, Tsq_un.CV, Tsq_un.PV) colnames(Tsq_un.out) <- c("Chi-Square Stat", "df", "Critical-Value", "P-value") round(Tsq_un.out,4) S12_mean <- (1/n.o)*S.o+(1/n.f)*S.f diff <- xbar.o - xbar.f sediff <- sqrt(diag(S12_mean)) simCILo <- diff - sqrt(Tsq_un.CV) * sediff simCIHi <- diff + sqrt(Tsq_un.CV) * sediff perfci_un <- cbind(diff,simCILo,simCIHi) colnames(perfci_un) <- c("Diff","Sim Lo","Sim Hi") rownames(perfci_un) <- c("Cl", "Na", "K", "Ethnol") round(perfci_un,3)