#   A split-split plot design with a randomized block design for the whole plots.
 #   This is Example 16.7 from Oehlerts text 
 #   "A First Course in Design and Analysis of Experiments" 
 #   Here nitrogen (n) is the whole plot factor with tables as blocks in a RCB design,
 #   weed (w) treatments are the split plot factor, 
 #   and clipping times (c) are the split-split factor


data <- matrix(scan(),ncol=7,byrow=TRUE)
    1    1    1    1    1    1   87.2
    1    1    1    2    1    1   88.8
    1    1    2    1    1    2   70.4
    1    1    2    2    1    2   75.7
    1    1    3    1    1    3   75.9
    1    1    3    2    1    3   80.6
    1    2    1    1    2    4   80.5
    1    2    1    2    2    4   83.8
    1    2    2    1    2    5   59.2
    1    2    2    2    2    5   61.5
    1    2    3    1    2    6   59.5
    1    2    3    2    2    6   62.5
    1    3    1    1    3    7   76.8
    1    3    1    2    3    7   80.8
    1    3    2    1    3    8   47.8
    1    3    2    2    3    8   49.5
    1    3    3    1    3    9   48.4
    1    3    3    2    3    9   52.9
    1    4    1    1    4   10   77.7
    1    4    1    2    4   10   81.5
    1    4    2    1    4   11   35.7
    1    4    2    2    4   11   37.3
    1    4    3    1    4   12   38.3
    1    4    3    2    4   12   42.4
    2    1    1    1    5   13   78.2
    2    1    1    2    5   13   80.5
    2    1    2    1    5   14   65.1
    2    1    2    2    5   14   68.3
    2    1    3    1    5   15   65.3
    2    1    3    2    5   15   66.6
    2    2    1    1    6   16   79.8
    2    2    1    2    6   16   85.2
    2    2    2    1    6   17   57.6
    2    2    2    2    6   17   61.4
    2    2    3    1    6   18   58.5
    2    2    3    2    6   18   61.6
    2    3    1    1    7   19   82.4
    2    3    1    2    7   19   83.1
    2    3    2    1    7   20   50.5
    2    3    2    2    7   20   54.0
    2    3    3    1    7   21   51.6
    2    3    3    2    7   21   54.7
    2    4    1    1    8   22   75.5
    2    4    1    2    8   22   78.7
    2    4    2    1    8   23   39.0
    2    4    2    2    8   23   43.9
    2    4    3    1    8   24   41.9
    2    4    3    2    8   24   45.1

data

RCBSplitsplit <- data.frame(data)
colnames(RCBSplitsplit) <- c("table","n","w","c","tray","wetland","biomass")

RCBSplitsplit$table   <- as.factor(RCBSplitsplit$table)
RCBSplitsplit$n       <- as.factor(RCBSplitsplit$n)
RCBSplitsplit$w       <- as.factor(RCBSplitsplit$w)
RCBSplitsplit$c       <- as.factor(RCBSplitsplit$c)

rm(data)


 # gives most of the same information as Proc Mixed for this Split-split plot design

library(lme4)
library(lmerTest)

RCBSpsp.lme1  <- lmer(biomass ~ n +w +c +n:w +n:c +w:c + n:w:c +
                      (1|table) +(1|n:table) +(1|n:w:table), data = RCBSplitsplit)

summary(RCBSpsp.lme1)

anova(RCBSpsp.lme1)

lsmeans(RCBSpsp.lme1)

difflsmeans(RCBSpsp.lme1)