#R commands for split-plot design, data in Table 14-14

y <- c(30, 35, 37, 36, 34, 41, 38, 42, 29, 26, 33, 36, 28, 32, 40, 41, 31,
  36, 42, 40, 31, 30, 32, 40, 31, 37, 41, 40, 35, 40, 39, 44, 32, 34, 39, 45)

R <- as.factor(rep(1:3, each=12))
A.pulp <- as.factor(rep(1:3, each = 4, times = 3))
B.temp <- as.factor(rep(c(200, 225, 250, 275), times=9))

g <- lm(y ~ (R/A.pulp) + (A.pulp + B.temp)^2)
anova(g)


MSA <- 64.19;   dfA <- 2  
MSRA <- 9.07;   dfRA <- 4
F.A <- MSA/MSRA
p.A <- 1-pf(F.A, dfA, dfRA)

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

data <- data.frame(R, A.pulp, B.temp, y)
par(mfrow=c(2,2))
plot.design(data)
interaction.plot(B.temp,A.pulp,y,lty=c(1,2,4))
interaction.plot(R,A.pulp,y,lty=c(1,2,4))
interaction.plot(B.temp,R,y,lty=c(1,2,4))



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

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

k<-lm(y~(A.pulp+B.temp)^2)#Gives MSA, MSB, MSAB
anova(k)
k2<-lm(k$resid~(R+A.pulp)^2)#Gives MSR, MSRA, SSE
#with SSA=0 (since A has already been accounted for).


anova(k2) #df(SSE) should be df(k$resid)-df(A)-df(RA)= 24-2-4=18.