if (FALSE) {
library(agridat)
data(yates.oats)
dat <- yates.oats
## # Means match Rothamsted report p. 144
## libs(dplyr)
## dat
## summarize(grain=mean(grain)*80/112,
## straw=mean(straw)*80/112)
libs(desplot)
# Experiment design & yield heatmap
desplot(dat, block ~ col*row, col.regions=c("black","yellow"),
out1=block, num=nitro, col=gen,
cex=1, aspect=511/176, # true aspect
main="yates.oats")
# Roughly linear gradient across the field. The right-half of each
# block has lower yield. The blocking is inadequate!
libs("lattice")
xyplot(yield ~ col|factor(nitro), dat,
type = c('p', 'r'), xlab='col', as.table = TRUE,
main="yates.oats")
libs(lme4)
# Typical split-plot analysis. Non-significant gen differences
m3 <- lmer(yield ~ factor(nitro) * gen + (1|block/gen), data=dat)
# Residuals still show structure
xyplot(resid(m3) ~ dat$col, xlab='col', type=c('p','smooth'),
main="yates.oats")
# Add a linear trend for column
m4 <- lmer(yield ~ col + factor(nitro) * gen + (1|block/gen), data=dat)
# xyplot(resid(m4) ~ dat$col, type=c('p','smooth'), xlab='col')
## Compare fits
AIC(m3,m4)
## df AIC
## m3 9 581.2372
## m4 10 557.9424 # Substantially better
# ----------
# Marginal predictions
# --- nlme ---
libs(nlme)
libs(emmeans)
# create unbalance
dat2 <- yates.oats[-c(1,2,3,5,8,13,21,34,55),]
m5l <- lme(yield ~ factor(nitro) + gen, random = ~1 | block/gen,
data = dat2)
# asreml r 4 has a bug with asreml( factor(nitro))
dat2$nitrof <- factor(dat2$nitro)
# --- asreml ---
if(require("asreml", quietly=TRUE)){
libs(asreml,lucid)
m5a <- asreml(yield ~ nitrof + gen,
random = ~ block + block:gen, data=dat2)
lucid::vc(m5l)
lucid::vc(m5a)
emmeans::emmeans(m5l, "gen")
predict(m5a, data=dat2, classify="gen")$pvals
}
}
Run the code above in your browser using DataLab