#R commands for split-plot design, data in Table 3.29
rm(list = ls())  # clear the memory
resist = read.table("http://www.stat.ualberta.ca/~wiens/stat568/datasets/WaterResistance_table_3_29.dat", header = F, quote= "")

A = as.factor(resist[,2])
B = as.factor(resist[,3])
R = as.factor(resist[,4])
y = resist[,5]


fit <- aov(y ~ (R/A) + (A + B)^2)
summary(fit)

MSA <- 782.04 ;   dfA <- 1  
MSRA <- 199.19;   dfRA <- 2
F.A <- MSA/MSRA
p.A <- 1-pf(F.A, dfA, dfRA)

cat("F.A =", F.A, "with p =",p.A,"\n")


MSB <- 88.97 ;   dfB <- 3  
MSE <- 12.71;   dfE <- 12
F.B <- MSB/MSE
p.B <- 1-pf(F.B, dfB, dfE)

cat("F.B =", F.B, "with p =",p.B,"\n")


# Here is a way to get everything from one command:
fit2 <- aov(y ~ (A + B)^2 + Error(R/A))
summary(fit2)

# The aov command has fitted a model y~(A + B)^2, 
# taken the residuals from this model, and fitted 
# R and R*A to these residuals to get the correct 
# value of SSE.