rm(list = ls())  # clear the memory
###########################

data = read.table("http://www.stat.ualberta.ca/~wiens/stat575/datasets/T7-3.DAT")
Z = as.matrix(data[, 1:2])
Y = cbind(data[,3], c(301.8, 396.1, 328.2, 307.4, 362.4, 369.5, 229.1))

m = ncol(Y)
r = ncol(Z)
n = nrow(Y)

colnames(Z) = c("orders", "items")
colnames(Y) = c("cpu", "in/out")
pairs(cbind(Y, Z), main = "cpu and in/out vs. orders and items")

dev.new()
library(car)
# If you haven't yet installed this package: click on Packages on the top menu of the screen, 
# choose a mirror (E.G. in BC) and install the "cars" package - it contains the scatterplotMatrix 
# and linearHypothesis functions
#scatterplotMatrix(cbind(Y, Z), main = "cpu and in/out vs. orders and items")

fit = lm(Y~Z)
summary(fit)

est.B = coef(fit)  # B-hat
est.B

E = residuals(fit)

W = t(E)%*%E   # W[2,2] in text seems to be wrong
S = W/(n-r-1)
S

#help(linearHypothesis)
#linearHypothesis(model, hypothesis.matrix, rhs=NULL,
#    test=c("F", "Chisq"), vcov.=NULL, 
#	white.adjust=c(FALSE, TRUE, "hc3", "hc0", "hc1", "hc2", "hc4"), 
#	singular.ok=FALSE, ...)

# linearHypothesis(fit, hypothesis.matrix, P) # Tests ABP = D; with A = hypothesis.matrix 
# (must be there, can be the identity); P is assumed to be the identity if it's not there 
# (and must be called P if it is there); default D (the "rhs") is zero
# e.g. linearHypothesis(fit, hypothesis.matrix = c(0,1,0), P = cbind(c(1,0)), verbose = TRUE)
  # -- Test that "orders" can be dropped in first response (cpu); compare with summary(fit)

# Test overall significance (i.e. that both "orders" and "items" can be dropped 
A = rbind(c(0,1,0), c(0,0,1))
linearHypothesis(fit, hypothesis.matrix = A, verbose = TRUE) 

# 95% Confidence ellipsoid when z = (1, 130, 7.5)'
z0 = as.vector(c(1, 130, 7.5))
Z0 = cbind(1, Z)
quad = t(z0)%*%solve(t(Z0)%*%Z0, z0)

alpha = .05
csqd = (m*(n-r-1)/(n-r-m))*qf(alpha, m, n-r-m, lower.tail = 0)
t = sqrt(quad*csqd)

# Confidence ellipsoid
U = chol(S)
phi = 2*pi*seq(from = 0, to = 1, length = 200)
v = rbind(t*cos(phi), t*sin(phi))  # 200 columns, each of norm t
gammahat = t(est.B)%*%z0
gamma = gammahat%*%rep(1,200) - t(U)%*%v
plot(gamma[1,], gamma[2,], xlim = c(144, 159), ylim = c(339, 361), xlab = "cpu", ylab = "in/out", main = "Confidence and prediction ellipsoids")
points(gammahat[1], gammahat[2], pch='*')

# Prediction ellipsoid
t2 = sqrt((1+quad)*csqd)
v2 = rbind(t2*cos(phi), t2*sin(phi))
gamma2 = gammahat%*%rep(1,200) - t(U)%*%v2
lines(gamma2[1,], gamma2[2,])