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