# MJ 24.1 Split Plot
spl <- read.table("split.dat",header=T)
spl$plot <- 4*(spl$block-1) + spl$fert
spl$fert <- factor(spl$fert)
spl$block <- factor(spl$block)
spl$variety <- factor(spl$variety)
spl$plot <- factor(spl$plot)

# split plot with block and block:fert as fixed effects
spl.fit <- aov(yield~block*fert+ variety+ variety:fert,spl)
# split plot using block:fert as random effect
# see Venables & Ripley (1994, sec. 6.7) or Chambers & Hastie (1992, sec. 5.2.1)
spl.bfit <- aov(yield~fert+Error(block+block:fert)+variety+variety:fert,spl)
print(summary(spl.bfit))

# split plot using LME (Splus on ALPHA computers only)
# see Lindstrom & Bates (1988) JASA 83:1014-1022
library(nlme)
spl.lme <- lme(yield~fert+variety+variety:fert,random=~block+block:fert,
	cluster=~plot, data=spl, est.method="RML", re.structure="identity")
print(summary(spl.lme))
print(fitted.values(spl.lme))