## Not run: 

  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, yield ~ col*row,
          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 from emmeans package and asreml::predict

  # --- 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)

  # --- asreml4  ---
  libs(asreml)
  m5a <- asreml(yield ~ nitrof + gen,
                random = ~ block + block:gen, data=dat2)
  libs(lucid)
  vc(m5l)
  vc(m5a)

  emmeans::emmeans(m5l, "gen")
  predict(m5a, data=dat2, classify="gen")$pvals

  # ----------

  if(0){

    # Demonstrate use of regress package, compare to lme

    libs(regress)
    m6 <- regress(yield ~ nitrof + gen, ~block + I(block:gen), identity=TRUE,
                  verbose=1, data=dat)
    summary(m6)
    ## Variance Coefficients:
    ##                Estimate Std. Error
    ##   block         214.468    168.794
    ##   I(block:gen)  109.700     67.741
    ##   In            162.558     32.189

    # ordinal causes clash with VarCorr
    if(is.element("package:ordinal", search())) detach(package:ordinal)

    m7 <- lme(yield ~ nitrof + gen, random = ~ 1|block/gen, data=dat)
    lme4::VarCorr(m7)
    ##             Variance     StdDev
    ## block =     pdLogChol(1)
    ## (Intercept) 214.4716     14.64485
    ## gen =       pdLogChol(1)
    ## (Intercept) 109.6929     10.47344
    ## Residual    162.5591     12.74987
  }


## End(Not run)