nfl2014 <- read.csv("http://www.stat.ufl.edu/~winner/sta4702/nfl_combine_2014.csv") attach(nfl2014); names(nfl2014) Height_OL <- Height[Position=="OL"] Weight_OL <- Weight[Position=="OL"] HandLen_OL <- HandLen[Position=="OL"] ArmLen_OL <- ArmLen[Position=="OL"] Height_DL <- Height[Position=="DL"] Weight_DL <- Weight[Position=="DL"] HandLen_DL <- HandLen[Position=="DL"] ArmLen_DL <- ArmLen[Position=="DL"] X1 <- cbind(Height_OL, Weight_OL, HandLen_OL, ArmLen_OL) X2 <- cbind(Height_DL, Weight_DL, HandLen_DL, ArmLen_DL) n1 <- nrow(X1); n2 <- nrow(X2); p <- ncol(X1) I_n1 <- diag(n1); J_n1 <- matrix(rep(1,n1^2),n1,n1) I_n2 <- diag(n2); J_n2 <- matrix(rep(1,n2^2),n2,n2) (xbar1 <- (1/n1) * (t(X1) %*% rep(1,n1))) (xbar2 <- (1/n2) * (t(X2) %*% rep(1,n2))) (S1 <- (1/(n1-1)) * (t(X1) %*% (I_n1 - (1/n1) * J_n1) %*% X1)) (S2 <- (1/(n2-1)) * (t(X2) %*% (I_n2 - (1/n2) * J_n2) %*% X2)) (Sp <- ((n1-1) * S1 + (n2-1) * S2) / (n1+n2-2)) ### Test mu1 = mu2 based on Hotelling's T^2 (Tsq <- t(xbar1-xbar2) %*% solve(((1/n1)+(1/n2)) * Sp) %*% (xbar1-xbar2)) (CV.wt <- (n1+n2-2)*p/(n1+n2-p-1)) (Tsq.CV <- CV.wt * qf(.95,p,n1+n2-p-1)) (Tsq.PV <- 1-pf(Tsq/CV.wt, p, n1+n2-p-1)) ### Simultaneous CI's for differences of mean characteristics diff <- xbar1 - xbar2 sediff <- sqrt(diag(Sp)*(1/n1+1/n2)) simCILo <- diff - sqrt(Tsq.CV) * sediff simCIHi <- diff + sqrt(Tsq.CV) * sediff bonCILo <- diff - qt(1-.05/(2*4),n1+n2-2) * sediff bonCIHi <- diff + qt(1-.05/(2*4),n1+n2-2) * sediff oldlci <- cbind(diff,simCILo,simCIHi,bonCILo,bonCIHi) colnames(oldlci) <- c("Diff","Sim Lo","Sim Hi", "Bon Lo", "Bon Hi") round(oldlci, 3) ### Linear combination with maximum population differerence (ahat <- solve(Sp) %*% (xbar1-xbar2)) ### Case When Sigma1 ^= Sigma2 (Tsq_un <- t(xbar1-xbar2) %*% solve((1/n1)*S1+(1/n2)*S2) %*% (xbar1-xbar2)) (Tsq_un.CV <- qchisq(.95,p)) (Tsq_un.PV <- 1-pchisq(Tsq_un,p)) S12_mean <- (1/n1)*S1+(1/n2)*S2 diff <- xbar1 - xbar2 sediff <- sqrt(diag(S12_mean)) simCILo <- diff - sqrt(Tsq_un.CV) * sediff simCIHi <- diff + sqrt(Tsq_un.CV) * sediff bonCILo <- diff - qt(1-.05/(2*4),n1+n2-2) * sediff bonCIHi <- diff + qt(1-.05/(2*4),n1+n2-2) * sediff oldlci_un <- cbind(diff,simCILo,simCIHi) colnames(oldlci_un) <- c("Diff","Sim Lo","Sim Hi") round(oldlci_un,3) ### Test for Equality of Covariance Matrices g <- 2 ### number of groups/treatments Lambda <- ((det(S1)/det(Sp))^((n1-1)/2)) * ((det(S2)/det(Sp))^((n2-1)/2)) M <- -2*log(Lambda) u1 <- (1/(n1-1)) + (1/(n2-1)) - (1/(n1+n2-2)) u2 <- (2*p^2+3*p-1) / (6*(p+1)*(g-1)) u <- u1*u2 (Box.C <- (1-u)*M) (df.X2 <- (1/2)*p*(p+1)*(g-1)) (Box.C.CV <- qchisq(.95,df.X2)) (Box.C.PV <- 1-pchisq(Box.C,df.X2))