# Compute main effects and interaction effects with the data in Table 6-10
# Open this file in R, click on 'Save', then on Edit -> Run all.

A <- as.factor(rep(c(-1,1), times=8))
B <- as.factor(rep(c(-1,1), each = 2, times=4))
C <- as.factor(rep(c(-1,1), each = 4, times=2))
D <- as.factor(rep(c(-1,1), each = 8))
y <- c(45, 71, 48, 65, 68, 60, 80, 65, 43, 100, 45, 104, 75, 86, 70, 96)

data <- data.frame(A, B, C, D, y)
data

g <- lm(y ~ (A+B+C+D)^4)  
# Easier than, but equivalent to, 
# y ~ A + B + C + D + A*B + A*C + A*D + B*C + B*D + C*D + A*B*C + A*B*D + B*C*D + A*B*C*D
anova(g)

g$effects
effects <- abs(g$effects[-1])


qq <- qqnorm(effects, type="n")  # "n" means no plotting
text(qq$x, qq$y, labels = names(effects))

h <- lm(y ~ (A+C+D)^3)
anova(h)

newdata <- data[,-2]
newdata

windows()
par(mfrow=c(2,2))
plot.design(newdata)
interaction.plot(A,C,y)
interaction.plot(A,D,y)
interaction.plot(C,D,y)

windows()
par(mfrow=c(3,2))
resids <- h$residuals
fits <- h$fitted.values
plot(c(A), resids)
abline(h=0)
plot(c(B), resids)
abline(h=0)
plot(c(C), resids)
abline(h=0)
plot(c(D), resids)
abline(h=0)
plot(fits, resids)
abline(h=0)
qqnorm(resids)
qqline(resids)