#R example; acetylene data
# Data from Montgomery & Peck Example 8.1
# Response variabe (conv) is % conversion of n-heptane to acetylene
# Explanatory variables temp (reactor temperature), mole (ratio of H2 to n-heptane), cont (contact time)
# Enter the data:
conv =c(49,50.2,50.5,48.5,47.5,44.5,28,31.5,34.5,35,38,38.5,15,17,20.5,29.5)
temp = c(rep(1300,6),rep(1200,6),rep(1100,4))
mole = c(7.5,9,11,13.5,17,23,5.3,7.5,11,13.5,17,23,5.3,7.5,11,17)
cont = c(120,120,115,130,135,120,400,380,320,260,340,410,840,980,920,860)/10000

# Put the data into a "frame"; look at all pairs of plots
acet = data.frame(conv, temp, mole, cont)
pairs(acet) # Note that the predictors cont and temp are highly correlated

# Fit a full second order model
x = cbind(temp, mole, cont, temp*mole, temp*cont, mole*cont, temp^2, mole^2, cont^2)
dimnames(x) = list(NULL, c("T", "M", "C", "T*M", "T*C", "M*C", "T2", "M2", "C2"))

fit1 = lsfit(x, conv)
ls.print(fit1)
d = ls.diag(fit1); print(d)
n = nrow(x)
p = ncol(x)+1
SSE.full = d$std.dev^2*(n-p)

# Reduced model, without columns 3,5,6,9 of X
fit2 = lsfit(x[ ,-c(3,5,6,9)], conv)
ls.print(fit2)
d2 = ls.diag(fit2); print(d2)
p2 = p-4
SSE.red = d2$std.dev^2*(n-p2)

F = ((SSE.red - SSE.full)/(p-p2))/(SSE.full/(n-p))
p.to.drop.C = 1-pf(F, p-p2, n-p)
cat("p-value of test to drop all terms involving 'C' is", p.to.drop.C, "\n")


fits = cbind(1,x)%*%fit1$coef
plot(fits,d$std.res, ylab="std.res")

xtx = crossprod(cbind(1,x))
det(xtx)
solve(xtx)

dd = d$cov.unscaled # This is xtx-inverse
dd%*%xtx
solve(dd)
diag(dd)

#> v = eigen(xtx)$values; v
# [1] 3.543840e+13 6.843691e+08 1.133772e+05 1.156081e+04 1.371599e+03 1.794139e+00
# [7] 1.046738e-02 4.513771e-04 3.203306e-05 2.608609e-07
#> max(v)/min(v)
#[1] 1.358517e+20