if (FALSE) {
  
  library(agridat)
  data(perry.springwheat)
  dat <- perry.springwheat
  libs(lattice)
  xyplot(yield~yor|env, dat, type=c('p','r'), xlab="year of release",
         main="perry.springwheat")
  
  # Show a table of sites*year
  # libs(latticeExtra)
  # useOuterStrips(xyplot(yield~yor|site*factor(year), dat,
  #                       type=c('p','r')))
  # Perry reports a rate of gain of 5.8 kg/ha/year.  No model is given.
  # We fit a model with separate intercept/slope for each env
  m1 <- lm(yield ~ env + yor + env:yor, data=dat)
  # Average slope across environments
  mean(c(coef(m1)[21], coef(m1)[21]+coef(m1)[22:40]))
  ## [1] 5.496781
  # ----------
  # Now a mixed-effects model.  Fixed overall int/slope.  Random env int/slope.
  # First, re-scale response so we don't have huge variances
  dat$y <- dat$yield / 100
  libs(lme4)
  # Use || for uncorrelated int/slope.  Bad model.  See below.
  # m2 <- lmer(y ~ 1 + yor + (1+yor||env), data=dat)
  ## Warning messages:
  ## 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  ##   Model failed to converge with max|grad| = 0.55842 (tol = 0.002, component 1)
  ## 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  ##   Model is nearly unidentifiable: very large eigenvalue
  ##  - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
  ##  - Rescale variables?
  # Looks like lme4 is having trouble with variance of intercepts
  # There is nothing special about 1800 years, so change the
  # intercept -- 'correct' yor by subtracting 1800 and try again.
  dat$yorc <- dat$yor - 1800
  m3 <- lmer(y ~ 1 + yorc + (1+yorc||env), data=dat)
  
  # Now lme4 succeeds.  Rate of gain is 100*0.0549 = 5.49
  fixef(m3)
  ## (Intercept)        yorc
  ##  5.87492444  0.05494464
  if(require("asreml", quietly=TRUE)){
    libs(asreml,lucid)
    m3a <- asreml(y ~ 1 + yorc, data=dat, random = ~ env + env:yorc)
    lucid::vc(m3)
    ##      grp        var1 var2     vcov   sdcor
    ##      env (Intercept)  11.61    3.407
    ##    env.1        yorc   0.00063 0.02511
    ## Residual           3.551   1.884
    
    lucid::vc(m3a)
    ##           effect component std.error z.ratio      con
    ##      env!env.var  11.61     4.385        2.6 Positive
    ## env:yorc!env.var   0.00063  0.000236     2.7 Positive
    ##       R!variance   3.551    0.2231      16   Positive
  }
  
}
Run the code above in your browser using DataLab