# 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)
# 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)

# Put the data into a "frame"; look at all pairs of plots
acet = data.frame(conv, temp, mole)
pairs(acet) 

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

fit1 = lm(conv~x)
summary(fit1)
d = ls.diag(fit1); print(d)

win.graph()
par(mfrow=c(1,2))
plot(fit1$fitted.values,d$stud.res, 
ylab="stud.res")
qqnorm(d$stud.res, main="QQ Plot") #qq plot for studentized resid 

n = nrow(x)
p = ncol(x)+1
SSE.full = d$std.dev^2*(n-p)

# Reduced model, without T or its interaction with M
x2 = x[ ,-c(1,3,4)]
dimnames(x2) = list(NULL, c("M", "M2"))
fit2 = lm(conv~x2)
summary(fit2)
d2 = ls.diag(fit2); print(d2)
p2 = p-3
SSE.red = d2$std.dev^2*(n-p2)

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