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,])