# Chapter 17 R Analysis of Split-Plot Experiments

library(lme4)
library(lmerTest)

y = fd$y b = factor(fd$block)
g = factor(fd$geno) f = factor(fd$fert/50 + 1)

o = lmer(y ~ g + f + g:f + (1|b) + (1|b:g))
summary(o)
anova(o)

ls_means(o)

betahat=fixef(o)
betahat

#Coefficients for geno 1 marginal mean
C1 = matrix(c(1, 0, 0, 1/4, 1/4, 1/4, 0, 0, 0, 0, 0, 0), nrow = 1)
#Coefficients for geno 1 - geno 2 marginal mean
C2 = matrix(c(0, -1, 0, 0, 0, 0, -1/4, 0, -1/4, 0, -1/4, 0), nrow = 1)
#Coefficients for geno 1 - geno 2 with no fertilizer
C3 = matrix(c(0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), nrow = 1)

C = rbind(C1, C2, C3)
contest(o, L = C, joint = F, confint = T)

a = anova(lm(y ˜ b + g + b:g + f + g:f))
a

MSbg = a[4,3]
MSe = a[6,3]
# an unbiased estimator of variance for the whole-plot random effects
(MSbg - MSe) / 4
MSg = a[2,3]
# The correct F statistic for testing for genotype main effects
MSg / MSbg

d = aggregate(y, by = list(b, g), FUN = mean)

anova(lm(wpaverage ˜ block + geno, data = d))