split_plot.R

library(lme4)
oats <- read.csv("oats.csv")
attach(oats)
nitrof <- factor(nitro)
              
xyplot(yield~nitrof|variety)
xyplot(yield~nitrof|variety,groups=replicate)
xyplot(yield~nitrof|variety,groups=replicate,aspect="xy",type="o")
xyplot(yield~variety|nitrof,groups=replicate,aspect="xy")
xyplot(yield~variety|nitrof,groups=replicate,aspect="xy",type="a")

res2.m <- lmer(yield~variety*nitrof+(1|replicate:variety))
summary(res2.m)
anova(res2.m)

1-pf(.6114,2,18-(3-1))  # Number of whole plots - (a-1)
1-pf(37.69,3,72-3*(4-1)-18) # Nt -a(b-1)-whole_plot_EU
1-pf(.30,6,72-3*(4-1)-18)

summary(res2.m)

n0 <- 80.00
n2 <- n0 + 18.5
n4 <- n0 + 34.6667
n6 <- n0 + 44.8333

lsd <- 1/sqrt(2)*qtukey(.95,4,72-3*(4-1)-18)*7.6830
lsd
n0-n2
n0-n4
n2-n4 # etc.

#
# Example 2
#

# There are types of grasses (cult) that make up the whole plot
# There three kinds of inoculation (inoc) treatments within each plot
# The response variable is dry weight of the grasses
# 8 fields were planted

data.2 <- read.csv("sp_grass.csv")
attach(data.2)

xyplot(drywt~cult|inoc,groups=rep,aspect="xy",type="o")
xyplot(drywt~inoc|cult,groups=rep,type="o")

head(data.2)  # Must make rep a factor

repf <- factor(rep)

res2.m <- lmer(drywt~factor(cult)*inoc+(1|cult:repf))

summary(res2.m)
anova(res2.m)

# Find p-values
# Find lsd for significant effects