dat <- kempton.slatehall
desplot(yield ~ col * row, dat, num=gen, out1=rep,
main="kempton.slatehall")
dat <- transform(dat, xf=factor(col), yf=factor(row))
# Incomplete block model of Gilmour et al 1995
require(lme4)
m1 <- lmer(yield ~ gen + (1|rep) + (1|rep:yf) + (1|rep:xf), data=dat)
print(VarCorr(m1), comp=c("Variance","Std.Dev."))
## Groups Name Variance Std.Dev.
## rep:xf (Intercept) 14811.6 121.703
## rep:yf (Intercept) 15595.0 124.880
## rep (Intercept) 4262.4 65.287
## Residual 8061.8 89.788
# Incomplete block model of Gilmour et al 1995
require("asreml")
m2 <- asreml(yield ~ gen, random = ~ rep/(xf+yf), data=dat)
summary(m2)$varcomp
# Table 4
predict(m2, classify="gen")$predictions$pvalsRun the code above in your browser using DataLab