# NOT RUN {
data(kempton.rowcol)
dat <- kempton.rowcol
dat <- transform(dat, rowf=factor(row), colf=factor(col))
if(require(desplot)){
desplot(yield~col*row|rep, dat,
num=gen, out1=rep, # unknown aspect
main="kempton.rowcol")
}
# }
# NOT RUN {
# Model with rep, row, col as random. Kempton, page 62.
# Use "-1" so that the vcov matrix doesn't include intercept
require(lme4)
m1 <- lmer(yield ~ -1 + gen + rep + (1|rep:rowf) + (1|rep:colf), data=dat)
# Variance components match Kempton.
print(m1, corr=FALSE)
# Standard error of difference for genotypes. Kempton page 62, bottom.
covs <- as.matrix(vcov(m1)[1:35, 1:35])
vars <- diag(covs)
vdiff <- outer(vars, vars, "+") - 2 * covs
sed <- sqrt(vdiff[upper.tri(vdiff)])
min(sed) # Minimum SED
mean(sed) # Average SED
max(sed) # Maximum SED
# }
# NOT RUN {
# }
Run the code above in your browser using DataCamp Workspace