library(lme4)
library(lmerTest)
fd = read.delim("<https://dnett.github.io/S510/FieldSplitPlotData.txt>")
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))