## R-code for assignment 5 solutions: ## AHW3 ## Nov 24, 2007 ###################################################################### ## 5.1 y <- matrix(c(2,12,8,9,6,9,8,10), 4,2,byrow=T) y ## n <- nrow(y) p <- ncol(y) y.bar <- apply(y, 2, mean) y.S <- cov(y) ## null hypothesis alpha <- 0.05 mu.0 <- matrix(c(7,11),2,1) mu.0 ## (a) T.sq.obs <- t(y.bar - mu.0)%*%solve(y.S/n)%*%(y.bar - mu.0) T.sq.obs ## (b) T.sq.crit <- (((n-1)*p)/(n-p))*qf(1-alpha, p, n-p) T.sq.crit ## (c) T.sq.obs > T.sq.crit #################################################################### ## 5.2 ## First run the test from example 5.1 y <- matrix(c(6,9,10,6,8,3), 3,2,byrow=T) y ## n <- nrow(y) p <- ncol(y) y.bar <- apply(y, 2, mean) y.S <- cov(y) ## null hypothesis alpha <- 0.05 mu.0 <- matrix(c(9,5),2,1) mu.0 T.sq.obs <- t(y.bar - mu.0)%*%solve(y.S/n)%*%(y.bar - mu.0) T.sq.obs T.sq.crit <- (((n-1)*p)/(n-p))*qf(1-alpha, p, n-p) T.sq.crit T.sq.obs > T.sq.crit ## Now re-run the test C <-matrix(c(1,-1,1,1),2,2, byrow=T) y.new <- matrix(0, n, p) for(i in 1:n){ y.new[i,] <- C%*%y[i,] } y.new ## n <- nrow(y.new) p <- ncol(y.new) y.bar <- apply(y.new, 2, mean) y.S <- cov(y.new) ## null hypothesis alpha <- 0.05 mu.0 <- C%*%matrix(c(9,5),2,1) mu.0 T.sq.obs <- t(y.bar - mu.0)%*%solve(y.S/n)%*%(y.bar - mu.0) T.sq.obs y.new <- matrix(0, n, p) for(i in 1:n){ y.new[i,] <- C%*%y[i,] } y.new ## n <- nrow(y.new) p <- ncol(y.new) y.bar <- apply(y.new, 2, mean) y.S <- cov(y.new) ## null hypothesis alpha <- 0.05 mu.0 <- C%*%matrix(c(9,5),2,1) mu.0 T.sq.obs <- t(y.bar - mu.0)%*%solve(y.S/n)%*%(y.bar - mu.0) T.sq.obs T.sq.crit <- (((n-1)*p)/(n-p))*qf(1-alpha, p, n-p) T.sq.crit T.sq.obs > T.sq.crit ####################################################################### ## 5.3 y <- matrix(c(2,12,8,9,6,9,8,10), 4,2,byrow=T) y ## n <- nrow(y) p <- ncol(y) y.bar <- apply(y, 2, mean) y.S <- cov(y) ## null hypothesis alpha <- 0.05 mu.0 <- matrix(c(7,11),2,1) mu.0 ## (a) SS.mu.0 <- 0 for(i in 1:n){ SS.mu.0 <- SS.mu.0 + (y[i,]-mu.0)%*%t(y[i,]-mu.0) } SS.y.bar <- 0 for(i in 1:n){ SS.y.bar <- SS.y.bar + (y[i,]-y.bar)%*%t(y[i,]-y.bar) } T.sq.obs <- ((n-1)*det(SS.mu.0))/(det(SS.y.bar)) - (n-1) T.sq.obs ## (b) Lambda <- (det(SS.y.bar)/det(SS.mu.0))^(n/2) Lambda ########################################################################### ## 5.4 y <- read.table("sweat.txt", header=T) y ## n <- nrow(y) p <- ncol(y) y.bar <- apply(y, 2, mean) y.S <- cov(y) ## (a) ## Construct the T.sq 90% ellipse axes ## calculate the eigen values and vectors eig <- eigen(y.S) eig.val.1 <- eig$values[1] eig.val.2 <- eig$values[2] eig.val.3 <- eig$values[3] eig.val.1 eig.val.2 eig.val.3 eig.vec.1 <- eig$vectors[,1] eig.vec.2 <- eig$vectors[,2] eig.vec.3 <- eig$vectors[,3] ## calculate the F-distribution modified by the appropriate constant alpha <- 0.10 mod.F <- (((n-1)*p)/((n-p)))*qf(1-alpha, p,n-p) mod.F ## determine the half lengths hl.1 <- sqrt(eig.val.1)*sqrt(mod.F/n) hl.2 <- sqrt(eig.val.2)*sqrt(mod.F/n) hl.3 <- sqrt(eig.val.3)*sqrt(mod.F/n) hl.1 hl.2 hl.3 ## (b) par(mfrow=c(3,1)) qqnorm(y[,1], pch=16, main="Normal Q-Q Plot: Sweat.Rate") qqnorm(y[,2], pch=16, main="Normal Q-Q Plot: Sodium") qqnorm(y[,3], pch=16, main="Normal Q-Q Plot: Potassium") pairs(y, pch=16) ######################################################################################## ## 5.9 y.bar <- matrix(c(95.52,164.38,55.69,93.39,17.98,31.13),6,1) y.S <- matrix(c(3266.46, 1343.97, 731.54, 1175.50, 162.68, 238.37, 1343.97, 721.91, 324.25, 537.35, 80.17, 117.73, 731.54, 324.25, 179.28, 281.17, 39.15, 56.80, 1175.50, 537.35, 281.17, 474.98, 63.73, 94.85, 162.68, 80.17, 39.15, 63.73, 9.95, 13.88, 238.37, 117.73, 56.80, 94.85, 13.88, 21.26),6,6) p <- 6 n <- 61 y.names <- c("weight", "body.length", "neck", "girth", "head.length", "head.width") ## (a) alpha <- 0.05 for(i in 1:6){ out.i <- c(y.bar[i,] - sqrt(qchisq(1-alpha, p))*sqrt(y.S[i,i]/n), y.bar[i,] + sqrt(qchisq(1-alpha, p))*sqrt(y.S[i,i]/n)) cat(y.names[i], round(out.i,3), "\n") } ## (b) y.bar.sub <- y.bar[c(1,4),] y.S.sub <- y.S[c(1,4), c(1,4)] alpha <- 0.05 eig.val <- eigen(y.S.sub)$values eig.vec <- eigen(y.S.sub)$vectors p.cen <- t(y.bar.sub) p.1 <- sqrt(qchisq(1-alpha, p)/n)*sqrt(eig.val[1])*eig.vec[,1] p.2 <- sqrt(qchisq(1-alpha, p)/n)*sqrt(eig.val[2])*eig.vec[,2] p.1 p.2 Ell <- rbind(p.cen, p.cen+p.1, p.cen+p.2, p.cen+(-p.1), p.cen+(-p.2)) Ell library(ellipse) plot(ellipse(y.S.sub, centre = c(y.bar.sub), t = sqrt(qchisq(1-alpha, p)/n)), lwd=3, xlab=paste("mu", y.names[1]), ylab=paste("mu", y.names[4])) points(Ell[,1], Ell[,2], pch=16, cex=3, col="blue") ## (c) alpha <- 0.05 for(i in 1:6){ out.i <- c(y.bar[i,] - qt(1-(alpha/(2*p)), n-1)*sqrt(y.S[i,i]/n), y.bar[i,] + qt(1-(alpha/(2*p)), n-1)*sqrt(y.S[i,i]/n)) cat(y.names[i], round(out.i,3), "\n") } ## (d) library(ellipse) plot(ellipse(y.S.sub, centre = c(y.bar.sub), t = sqrt(qchisq(1-alpha, p)/n)), lwd=3, xlab=paste("mu", y.names[1]), ylab=paste("mu", y.names[4])) points(Ell[,1], Ell[,2], pch=16, cex=3, col="blue") out.1 <- c(y.bar[1,] - qt(1-(alpha/(2*p)), n-1)*sqrt(y.S[1,1]/n), y.bar[1,] + qt(1-(alpha/(2*p)), n-1)*sqrt(y.S[1,1]/n)) out.4 <- c(y.bar[4,] - qt(1-(alpha/(2*p)), n-1)*sqrt(y.S[4,4]/n), y.bar[4,] + qt(1-(alpha/(2*p)), n-1)*sqrt(y.S[4,4]/n)) abline(v=out.1, col="green", lwd=2) abline(h=out.4, col="green", lwd=2) ## (e) y.bar.min <- y.bar[6,]-y.bar[5,] s.star <- y.S[6,6]+y.S[5,5]-2*y.S[6,5] round(c(y.bar.min - qt(1-(alpha/(2*(p+1))), n-1)*sqrt(s.star/n), y.bar.min + qt(1-(alpha/(2*(p+1))), n-1)*sqrt(s.star/n)),3) ####################################################################### ## 5.10 y <- matrix(c(141, 157, 168, 183, 140, 168, 174, 170, 145, 162, 172, 177, 146, 159, 176, 171, 150, 158, 168, 175, 142, 140, 178, 189, 139, 171, 176, 175), 7, 4, byrow=T) y.bar <- matrix(apply(y,2,mean),4,1) y.S <- cov(y) n <- nrow(y) p <- ncol(y) ## (a) alpha <- 0.05 mod.F <- ((p*(n-1))/(n-p))*qf(1-alpha, p, n-p) for(i in 1:4){ out.i <- c(y.bar[i,] - sqrt(mod.F)*sqrt(y.S[i,i]/n), y.bar[i,] + sqrt(mod.F)*sqrt(y.S[i,i]/n)) cat("Length",i+1, round(out.i,3), "\n") } ## (b) A <- matrix(c(-1,1,0,0, 0,-1,1,0, 0,0,-1,1),3,4,byrow=T) y.bar.diff <- A%*%y.bar y.S.diff <- A%*%y.S%*%t(A) for(i in 1:3){ out.i <- c(y.bar.diff[i,] - sqrt(mod.F)*sqrt(y.S.diff[i,i]/n), y.bar.diff[i,] + sqrt(mod.F)*sqrt(y.S.diff[i,i]/n)) cat("Length",i+1, "-","Length", i, ":",round(out.i,3), "\n") } ## (c) A <- matrix(c(-1,1,0,0, 0,0,-1,1),2,4,byrow=T) y.bar.diff <- A%*%y.bar y.S.diff <- A%*%y.S%*%t(A) alpha <- 0.05 eig.val <- eigen(y.S.diff)$values eig.vec <- eigen(y.S.diff)$vectors p.cen <- t(y.bar.diff) p.1 <- sqrt(mod.F/n)*sqrt(eig.val[1])*eig.vec[,1] p.2 <- sqrt(mod.F/n)*sqrt(eig.val[2])*eig.vec[,2] p.1 p.2 Ell <- rbind(p.cen, p.cen+p.1, p.cen+p.2, p.cen+(-p.1), p.cen+(-p.2)) Ell mod.F library(ellipse) plot(ellipse(y.S.diff/n, centre = c(y.bar.diff), t = sqrt(mod.F)), lwd=3, xlab=paste("mu 3-2"), ylab=paste("mu 5-4")) points(Ell[,1], Ell[,2], pch=16, cex=3, col="blue") ## (d) A <- matrix(c(1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1, -1,1,0,0, 0,-1,1,0, 0,0,-1,1),7,4,byrow=T) y.comb.names <- c("L2:", "L3:", "L4:", "L5:", "L3-L2:", "L4-L3:", "L5-L4:") y.bar.comb <- A%*%y.bar y.S.comb <- A%*%y.S%*%t(A) ## note this matrix is not a proper covariance matrix b/c it is not invertible!! for(i in 1:7){ out.i <- c(y.bar.comb[i,] - qt(1-(alpha/(2*(4+3))), n-1)*sqrt(y.S.comb[i,i]/n), y.bar.comb[i,] + qt(1-(alpha/(2*(4+3))), n-1)*sqrt(y.S.comb[i,i]/n)) cat(y.comb.names[i], round(out.i,3), "\n") } ## (e) A <- matrix(c(-1,1,0,0, 0,0,-1,1),2,4,byrow=T) y.bar.diff <- A%*%y.bar y.S.diff <- A%*%y.S%*%t(A) alpha <- 0.05 eig.val <- eigen(y.S.diff)$values eig.vec <- eigen(y.S.diff)$vectors p.cen <- t(y.bar.diff) p.1 <- sqrt(mod.F/n)*sqrt(eig.val[1])*eig.vec[,1] p.2 <- sqrt(mod.F/n)*sqrt(eig.val[2])*eig.vec[,2] p.1 p.2 Ell <- rbind(p.cen, p.cen+p.1, p.cen+p.2, p.cen+(-p.1), p.cen+(-p.2)) Ell mod.F library(ellipse) plot(ellipse(y.S.diff/n, centre = c(y.bar.diff), t = sqrt(mod.F)), lwd=3, xlab=paste("mu 3-2"), ylab=paste("mu 5-4"), ylim=c(-25,60), xlim=c(-25,60)) points(Ell[,1], Ell[,2], pch=16, cex=3, col="blue") ## add the Bonferroni intervals abline(v=c(-1.423, 33.423), lwd=2, col="green") abline(h=c(-7.539, 15.539), lwd=2, col="green") ## add the T^2 intervals abline(v=c(-21.226, 53.226), lwd=2, col="blue") abline(h=c(-20.654, 28.654), lwd=2, col="blue") ############################################################## ## 5.11 y1 <- c(0.48, 40.53, 2.19, 0.55, 0.74, 0.66, 0.93, 0.37, 0.22) y2 <- c(12.57, 73.68, 11.13, 20.03, 20.29, 0.78, 4.64, 0.43, 1.08) y <- cbind(y1,y2) y.bar <- matrix(apply(y, 2, mean), 2, 1) y.S <- cov(y) n <- nrow(y) p <- ncol(y) ## (a) alpha <- 0.10 mod.F <- ((p*(n-1))/(n-p))*qf(1-alpha, p, n-p) eig.val <- eigen(y.S)$values eig.vec <- eigen(y.S)$vectors p.cen <- t(y.bar) p.1 <- sqrt(mod.F/n)*sqrt(eig.val[1])*eig.vec[,1] p.2 <- sqrt(mod.F/n)*sqrt(eig.val[2])*eig.vec[,2] p.1 p.2 Ell <- rbind(p.cen, p.cen+p.1, p.cen+p.2, p.cen+(-p.1), p.cen+(-p.2)) Ell mod.F library(ellipse) plot(ellipse(y.S/n, centre = c(y.bar), t = sqrt(mod.F)), lwd=3, xlab=paste("mu 1"), ylab=paste("mu 2")) points(Ell[,1], Ell[,2], pch=16, cex=3, col="blue") ## (b) out <- NULL for(i in 1:2){ out.i <- c(y.bar[i,] - sqrt(mod.F)*sqrt(y.S[i,i]/n), y.bar[i,] + sqrt(mod.F)*sqrt(y.S[i,i]/n)) out <- rbind(out, out.i) } ## The confidence intervals out ## add the T^2 intervals abline(v=c(out[1,]), lwd=2, col="blue") abline(h=c(out[2,]), lwd=2, col="blue") ## add the specific point points(0.30, 10, col="red", cex=2, pch=16) ## (c) par(mfrow=c(2,1)) qqnorm(y[,1], pch=16) qqnorm(y[,2], pch=16) pairs(y, pch=16) ## (d) ## remove the row from the data y <- y[-2,] y.bar <- matrix(apply(y, 2, mean), 2, 1) y.S <- cov(y) n <- nrow(y) p <- ncol(y) alpha <- 0.10 mod.F <- ((p*(n-1))/(n-p))*qf(1-alpha, p, n-p) eig.val <- eigen(y.S)$values eig.vec <- eigen(y.S)$vectors p.cen <- t(y.bar) p.1 <- sqrt(mod.F/n)*sqrt(eig.val[1])*eig.vec[,1] p.2 <- sqrt(mod.F/n)*sqrt(eig.val[2])*eig.vec[,2] p.1 p.2 Ell <- rbind(p.cen, p.cen+p.1, p.cen+p.2, p.cen+(-p.1), p.cen+(-p.2)) Ell mod.F par(mfrow=c(1,1)) library(ellipse) plot(ellipse(y.S/n, centre = c(y.bar), t = sqrt(mod.F)), lwd=3, xlab=paste("mu 1"), ylab=paste("mu 2")) points(Ell[,1], Ell[,2], pch=16, cex=3, col="blue") out <- NULL for(i in 1:2){ out.i <- c(y.bar[i,] - sqrt(mod.F)*sqrt(y.S[i,i]/n), y.bar[i,] + sqrt(mod.F)*sqrt(y.S[i,i]/n)) out <- rbind(out, out.i) } ## The confidence intervals out ## add the T^2 intervals abline(v=c(out[1,]), lwd=2, col="blue") abline(h=c(out[2,]), lwd=2, col="blue") ## add the specific point points(0.30, 10, col="red", cex=2, pch=16) ######################################################### ## 6.1 d.bar <- matrix(c(-9.36, 13.27),2,1) S.d <- matrix(c(199.26, 88.38, 88.38, 418.61),2,2) n <- 11 p <- 2 alpha <- 0.05 mod.F <- ((p*(n-1))/(n-p))*qf(1-alpha, p, n-p) eig.val <- eigen(S.d)$values eig.vec <- eigen(S.d)$vectors p.cen <- t(d.bar) p.1 <- sqrt(mod.F/n)*sqrt(eig.val[1])*eig.vec[,1] p.2 <- sqrt(mod.F/n)*sqrt(eig.val[2])*eig.vec[,2] p.1 p.2 Ell <- rbind(p.cen, p.cen+p.1, p.cen+p.2, p.cen+(-p.1), p.cen+(-p.2)) Ell mod.F par(mfrow=c(1,1)) library(ellipse) plot(ellipse(S.d/n, centre = c(d.bar), t = sqrt(mod.F)), lwd=3, xlab=paste("mu diff 1"), ylab=paste("mu diff 2")) points(Ell[,1], Ell[,2], pch=16, cex=3, col="blue") out <- NULL for(i in 1:2){ out.i <- c(d.bar[i,] - sqrt(mod.F)*sqrt(S.d[i,i]/n), d.bar[i,] + sqrt(mod.F)*sqrt(S.d[i,i]/n)) out <- rbind(out, out.i) } ## add the T^2 intervals abline(v=c(out[1,]), lwd=2, col="blue") abline(h=c(out[2,]), lwd=2, col="blue") ## add the specific point points(0, 0, col="red", cex=2, pch=16)