# Data from "Analysis of Messy Data, Vol. 1" by Milliken and Johnson 2nd edition page 477 ;
 # see the diagram for this lecture to show the layout of the experiment ;


data <- matrix(scan(),ncol=4,byrow=TRUE)
 1 1 1 55  1 1 2 71  1 3 1 69  1 3 2 78  1 2 1 62  1 2 2 77
 2 3 2 81  2 3 1 77  2 1 2 77  2 1 1 63  2 2 2 79  2 2 1 66
 3 2 1 70  3 2 2 78  3 3 1 79  3 3 2 80  3 1 1 63  3 1 2 77
 4 3 1 76  4 3 2 79  4 2 1 66  4 2 2 76  4 1 1 65  4 1 2 75
  

data

SplitBlk <- data.frame(data)
colnames(SplitBlk) <- c("plot","nitro","irrig","y")

SplitBlk$plot   <- as.factor(SplitBlk$plot)
SplitBlk$nitro       <- as.factor(SplitBlk$nitro)
SplitBlk$irrig       <- as.factor(SplitBlk$irrig)

rm(data)


 # using the afex package 

library(afex)

SplitBlk.afex1  <- mixed(y ~ nitro +irrig +nitro:irrig  +
                      (1|plot) +(1|nitro:plot) +(1|irrig:plot), data = SplitBlk)

summary(SplitBlk.afex1)

nice(SplitBlk.afex1)

lsmip(SplitBlk.afex1, nitro ~ irrig)

lsmip(SplitBlk.afex1, irrig ~ nitro)