# This file is posted as  grass.ssc

#  Enter the data into a data frame and 
#  define factors

grass <- read.table("grass.dat", 
   col.names=c("Block","Culti","Innoc","Yield"))
grass$Block <- as.factor(grass$Block)
grass$Culti <- as.factor(grass$Culti)
grass$Innoc <- as.factor(grass$Innoc)

# Multistratum models may be fitted using aov() if
# the design is balanced and all factors are random. 
# Models are specified by a model formula
#
#  response ~ mean.formula + Error(strata.formular)

summary(grass.aov <- aov(Yield ~ Culti*Innoc + 
          Error(Block/Culti),data=grass))


# Use the lme function to fit the split-plot
# model with random block effects.  Here
# each cultivar is used only once in each
# block, so cultivars within blocks can be
# used to denote the whole plots within 
# blocks.  The code  Block/Culti  in the
# random statement designates random
# block effects and random effects for
# cultivars within blocks (whole plots),  

options(contrasts=c(factor="contr.treatment",
                    ordered="contr.poly"))
grass.lme <- lme( Yield ~ Culti*Innoc , data=grass,
    random = ~ 1 | Block/Culti )

# Print F-tests for fixed effects
anova(grass.lme)

# Print estimates of variance components
VarCorr(grass.lme)

# Print estimates of fixed effects
fixef(grass.lme)

# Print predictions of random effects
ranef(grass.lme)

# Print approximate confidence intervals
intervals(grass.lme)

# Plot residuals versus fitted values
par(fin=c(6,6),cex=1.5,mkh=1.2,mex=1.5)
plot(grass.lme, resid(.,type="p") ~ fitted(.),
      id=0.01, adj=-0.3)

# Plot observations versus fitted values
plot(grass.lme, Yield ~ fitted(.),
      id=0.01, adj=-0.3)

# Normal probability plot of within whole
# plot residuals
qqnorm( grass.lme, ~resid(.))

# Normal probability plot of random
# effects
qqnorm( grass.lme, ~ranef(.,level=1), id=.01)
qqnorm( grass.lme, ~ranef(.,level=2), id=.01)