rm(list = ls()) # clear the memory ########################### X = as.matrix(read.table("http://www.stat.ualberta.ca/~wiens/stat575/datasets/T5-2.DAT")) n = nrow(X) p = ncol(X) colnames(X) = c("SS&H", "Vb", "Sc") xbar = colMeans(X) S = cov(X) S ######################################## # Test marginal normality. par(mfrow = c(2,2)) # A 2 by 2 panel of plots for(i in 1:3) { y = X[,i] v=qqnorm(y, ylab = colnames(X)[i]) text(0, .9*max(v\$y), paste("p = ", round(shapiro.test(y)[[2]],3))) qqline(y) } # Trivariate normality: dsqd = vector(length = n) qsqd = vector(length = n) for(i in 1:n) { dsqd[i] = t(X[i,] - xbar)%*%solve(S,(X[i,] - xbar)) qsqd[i] = qchisq((i-.5)/n, p, lower.tail = T) } dsqd = sort(dsqd) qqplot(qsqd, dsqd, main = "Chisquare Q-Q Plot", xlab = "Chisquare quantiles", ylab = "sample quantiles") abline(0,1) text(6, max(qsqd-2), paste("corr = ", round(cor(qsqd,dsqd),3))) ######################################## # Bonferroni intervals on the three marginal means; 95% alpha = .05 m = 3 alpham = alpha/m int = matrix(nrow = m, ncol = 2, data = NA) for(i in 1:m) { psihat = xbar[i] sd_psihat = sqrt(S[i,i]/n) qT = qt(alpham/2, n-1, lower.tail = 0) int[i,] = c(psihat - qT*sd_psihat, psihat + qT*sd_psihat) } int_Bon = cbind(int[,1], xbar, int[,2]) colnames(int_Bon) = c("lower", "xbar", "upper") cat("Bonferroni intervals on the three marginal means","\n") int_Bon ######################################## # Simultaneous intervals on these AND ALL OTHERS alpha = .05 m = 3 int = matrix(nrow = m, ncol = 2, data = NA) for(i in 1:m) { I = diag(m) a = I[,i] # Any other vectors a could be used without lowering the overall confidence psihat = t(a)%*%xbar sd_psihat = sqrt(t(a)%*%S%*%a/n) qF = sqrt(((n-1)*p/(n-p))*qf(alpha, p, n-p, lower.tail = 0)) int[i,] = c(psihat - qF*sd_psihat, psihat + qF*sd_psihat) } int_Scheffe = cbind(int[,1], xbar, int[,2]) colnames(int_Scheffe) = c("lower", "xbar", "upper") int_Scheffe # Values in text not very accurate - lots of rounding ######################################## # 95% confidence ellipsoid on means of VB and Sc: X23 = X[, 2:3] n = nrow(X23) p = ncol(X23) alpha = .05 xbar23 = colMeans(X23) S23 = cov(X23) U = chol(S23) # S23 = U'U c = sqrt(((n-1)*p/(n-p))*qf(alpha, p, n-p, lower.tail = 0)) phi = 2*pi*seq(from = 0, to = 1, length = 201) # 2*pi*(0,1/200, 2/200, ... 199/200, 1) z = rbind(cos(phi), sin(phi)) # 201 columns, each of norm 1 mu = xbar23 + (c/sqrt(n))*t(U)%*%z dev.new() plot(mu[1,], mu[2,], xlab = "mean Vb", ylab = "mean Sc", type = 'l') ######################################## # Simultaneous intervals on these two means, and all other comparisons of these two means only; # should agree with the extremes of the confidence ellipsoid: alpha = .05 m = 2 int = matrix(nrow = m, ncol = 2, data = NA) for(i in 1:m) { I = diag(m) a = I[,i] psihat = t(a)%*%xbar23 sd_psihat = sqrt(t(a)%*%S23%*%a/n) qF = sqrt(((n-1)*p/(n-p))*qf(alpha, p, n-p, lower.tail = 0)) int[i,] = c(psihat - qF*sd_psihat, psihat + qF*sd_psihat) } int2_Scheffe = cbind(int[,1], xbar23, int[,2]) colnames(int2_Scheffe) = c("lower", "xbar", "upper") int2_Scheffe abline(v=int[1,]) # vertical lines at the endpoints of the interval on mean(Vb) abline(h=int[2,]) # horizontal lines at the endpoints of the interval on mean(Sc) extremes_of_ellipsoid = rbind(c(min(mu[1,]), max(mu[1,])),c(min(mu[2,]), max(mu[2,]))) extremes_of_ellipsoid ######################################## # 95% prediction region for a new (VB,Sc): c = sqrt(((n-1)*p/(n-p))*qf(alpha, p, n-p, lower.tail = 0)) c = c*sqrt((n+1)/n) phi = 2*pi*seq(from = 0, to = 1, length = 201) # 2*pi*(0,1/200, 2/200, ... 199/200, 1) z = rbind(cos(phi), sin(phi)) # 201 columns, each of norm 1 mu = xbar23 + (c/sqrt(n))*t(U)%*%z dev.new() plot(mu[1,], mu[2,], xlab = "Vb", ylab = "Sc", type = 'l') ## Superimpose the previous confidence ellipsoid on the means: c = sqrt(((n-1)*p/(n-p))*qf(alpha, p, n-p, lower.tail = 0)) phi = 2*pi*seq(from = 0, to = 1, length = 201) # 2*pi*(0,1/200, 2/200, ... 199/200, 1) mu = xbar23 + (c/sqrt(n))*t(U)%*%z #dev.new() lines(mu[1,], mu[2,], type = 'l') ######################################## ## Test mu = c(500, 50, 25) using Tsqd and the large sample test p = 3 mu0 = c(500, 50, 25) Tsqd = n*t(xbar-mu0)%*%solve(S,xbar-mu0) Fcalc = (n-p)*Tsqd/(p*(n-1)) pval = pf(Fcalc, p, n-p, lower.tail = 0) # prob(F > Fcalc) Lambda = (1+Tsqd/(n-1))^(-n/2) chisqd.calc = -2*log(Lambda) pval.approx = pchisq(chisqd.calc, p, lower.tail = 0) pval pval.approx