Split-plot experiment of oats
yates.oats.RdThe yield of oats from a split-plot field trial conducted at Rothamsted in 1931.
Varieties were applied to the main plots.
Manurial (nitrogen) treatments were applied to the sub-plots.
Each plot is 1/80 acre = 28.4 links * 44 links.
Field width: 4 plots * 44 links = 176 links.
Field length: 18 rows * 28.4 links = 511 links
The 'block' numbers in this data are as given in the Rothamsted Report. The 'grain' and 'straw' values are the actual pounds per sub-plot as shown in the Rothamsted Report. Each sub-plot is 1/80 acre, and a 'hundredweight (cwt)' is 112 pounds, so converting from sub-plot weight to hundredweight/acre needs a conversion factor of 80/112.
The 'yield' values are the values as they appeared in the paper by Yates, who used 1/4-pounds as the units (i.e. he multiplied the original weight by 4) for simpler calculations.
Format
rowrow
colcolumn
yieldyield in 1/4 pounds per sub-plot, each 1/80 acre
nitronitrogen treatment in hundredweight per acre
gengenotype, 3 levels
blockblock, 6 levels
graingrain weight in pounds per sub-plot
strawstraw weight in pounds per sub-plot
Source
Report for 1931. Rothamsted Experiment Station. Page 143. https://www.era.rothamsted.ac.uk/eradoc/article/ResReport1931-141-159
References
Yates, Frank (1935) Complex experiments, Journal of the Royal Statistical Society Suppl 2, 181-247. Figure 2. https://doi.org/10.2307/2983638
Examples
# \dontrun{
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")
#> Warning: pseudoinverse used at 2
#> Warning: neighborhood radius 1
#> Warning: reciprocal condition number 0
#> Warning: There are other near singularities as well. 1
#> Warning: pseudoinverse used at 2
#> Warning: neighborhood radius 1
#> Warning: reciprocal condition number 0
#> Warning: There are other near singularities as well. 1
#> Warning: pseudoinverse used at 2
#> Warning: neighborhood radius 1
#> Warning: reciprocal condition number 0
#> Warning: There are other near singularities as well. 1
#> Warning: pseudoinverse used at 2
#> Warning: neighborhood radius 1
#> Warning: reciprocal condition number 0
#> Warning: There are other near singularities as well. 1
#> Warning: pseudoinverse used at 2
#> Warning: neighborhood radius 1
#> Warning: reciprocal condition number 0
#> Warning: There are other near singularities as well. 1
# 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 15 559.0285
#> m4 16 533.3343
## 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)
#> Model fitted using the gamma parameterization.
#> ASReml 4.1.0 Fri Dec 17 15:17:52 2021
#> LogLik Sigma2 DF wall cpu
#> 1 -202.393 297.150 57 15:17:52 0.0
#> 2 -198.059 229.020 57 15:17:52 0.0
#> 3 -194.582 178.247 57 15:17:52 0.0
#> 4 -193.036 152.265 57 15:17:52 0.0
#> 5 -192.591 140.059 57 15:17:52 0.0
#> 6 -192.551 137.448 57 15:17:52 0.0
#> 7 -192.550 137.136 57 15:17:52 0.0
libs(lucid)
vc(m5l)
#> effect variance stddev
#> block = NA NA
#> (Intercept) 186 13.64
#> gen = NA NA
#> (Intercept).1 179.6 13.4
#> Residual 137.1 11.71
vc(m5a)
#> effect component std.error z.ratio bound %ch
#> block 186 167.4 1.1 P 0.1
#> block:gen 179.6 99.63 1.8 P 0
#> units!R 137.1 29.93 4.6 P 0
emmeans::emmeans(m5l, "gen")
#> gen emmean SE df lower.CL upper.CL
#> GoldenRain 103.4 8.37 5 81.9 125
#> Marvellous 108.9 8.18 5 87.8 130
#> Victory 97.3 8.18 5 76.3 118
#>
#> Results are averaged over the levels of: nitro
#> Degrees-of-freedom method: containment
#> Confidence level used: 0.95
predict(m5a, data=dat2, classify="gen")$pvals
#> Model fitted using the gamma parameterization.
#> ASReml 4.1.0 Fri Dec 17 15:17:52 2021
#> LogLik Sigma2 DF wall cpu
#> 1 -192.550 137.129 57 15:17:52 0.0
#> 2 -192.550 137.129 57 15:17:52 0.0
#> 3 -192.550 137.129 57 15:17:52 0.0
#>
#> Notes:
#> - The predictions are obtained by averaging across the hypertable
#> calculated from model terms constructed solely from factors in
#> the averaging and classify sets.
#> - Use 'average' to move ignored factors into the averaging set.
#> - The simple averaging set: nitrof
#> - The ignored set: block
#>
#>
#> gen predicted.value std.error status
#> 1 GoldenRain 103.42039 8.368025 Estimable
#> 2 Marvellous 108.86384 8.183540 Estimable
#> 3 Victory 97.30395 8.183515 Estimable
# ----------
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
}
# }