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"] nfl2014_OL <- data.frame(Height_OL, Weight_OL, HandLen_OL, ArmLen_OL) detach(nfl2014) attach(nfl2014_OL) #### Data now only contains players who are offensive linemen (OL) (n_OL <- length(Height_OL)) ### n=sample size (number of players) one_OL <- rep(1,n_OL) ### nx1 vecor of 1s ### create X matrix by column binding the p=2 variables Ht and Wt X_OL <- cbind(Height_OL, Weight_OL) ### Obtain px1 vector xbar = (1/n) X'1_n * is scalar multiply, %*% is matrix (xbar_OL <- (1/n_OL) * t(X_OL) %*% one_OL) ### J_n is nxn matrix of 1s, I_n is nxn Identity matrix J_n <- matrix(rep(1,n_OL^2),ncol=n_OL) I_n <- diag(n_OL) ### Deviation matrix = (I - (1/n)J)X dmat_OL <- (I_n - (1/n_OL)*J_n) %*% X_OL ### Sum of Squares and Crossproducts matrix: SSCP=X'(I - (1/n)J)X and S (SS_dmat_OL <- t(dmat_OL) %*% dmat_OL) (S_dmat_OL <- (1/(n_OL-1)) * SS_dmat_OL) ### Test H0: mu'=m0'=[76 300] mu0 <- c(76, 300) (T2 <- n_OL * t(xbar_OL - mu0) %*% solve(S_dmat_OL) %*% (xbar_OL - mu0)) (T2_RR <- (2*(n_OL-1)/(n_OL-2)) * qf(.95,2,n_OL-2)) ### Obtain eigenvalues (lambda) and eigenvectors (as matrix P) of S eigen_OL <- eigen(S_dmat_OL) (lambda_OL <- eigen_OL$val) (P_OL <- eigen_OL$vec) ### P is [e1 ... ep] (pxp) ### Obtain diagonal matrix of eigenvalues and sqrt(eigvals) (Lambda_OL <- diag(lambda_OL)) (Lambda12_OL <- diag(sqrt(lambda_OL))) ### Simultaneous CI for Mean HT,WT (SimHtLo <- xbar_OL[1] - sqrt(((2*49)/(50*48))*qf(.95,2,48)*S_dmat_OL[1,1])) (SimHtHi <- xbar_OL[1] + sqrt(((2*49)/(50*48))*qf(.95,2,48)*S_dmat_OL[1,1])) (SimWtLo <- xbar_OL[2] - sqrt(((2*49)/(50*48))*qf(.95,2,48)*S_dmat_OL[2,2])) (SimWtHi <- xbar_OL[2] + sqrt(((2*49)/(50*48))*qf(.95,2,48)*S_dmat_OL[2,2])) ## Bonferroni CIs for Mean HT,WT (BonHtLo <- xbar_OL[1] - qt(1-.05/4,49)*sqrt(S_dmat_OL[1,1]/50)) (BonHtHi <- xbar_OL[1] + qt(1-.05/4,49)*sqrt(S_dmat_OL[1,1]/50)) (BonWtLo <- xbar_OL[2] - qt(1-.05/4,49)*sqrt(S_dmat_OL[2,2]/50)) (BonWtHi <- xbar_OL[2] + qt(1-.05/4,49)*sqrt(S_dmat_OL[2,2]/50)) #### Points a constant distance from the origin: drive distance/fairway pct x1 <- Height_OL; x2 <- Weight_OL n <- length(x1); p <- 2 alpha <- 0.05 crit.dist <- sqrt((p*(n-1)/(n*(n-p)))*qf(1-alpha,p,n-p)) A <- rbind(cbind(var(x1),cov(x1,x2)),cbind(cov(x1,x2),var(x2))) mu.test <- c(76,300) ctr <- c(mean(x1),mean(x2)) angles <- seq(0, 2*pi, length.out=200) eigVal <- eigen(A)$values eigVec <- eigen(A)$vectors eigScl <- eigVec %*% diag(sqrt(eigVal)) # scale eigenvectors to length = square-root 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="x1", ylab="x2") 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="green") bon11 <- ctr[1] - qt(1-alpha/(2*p),n-1)*sqrt(A[1,1]/n) bon12 <- ctr[1] + qt(1-alpha/(2*p),n-1)*sqrt(A[1,1]/n) bon21 <- ctr[2] - qt(1-alpha/(2*p),n-1)*sqrt(A[2,2]/n) bon22 <- ctr[2] + qt(1-alpha/(2*p),n-1)*sqrt(A[2,2]/n) rect(bon11,bon21,bon12,bon22,border="violetred3")