Perez-Elizalde et al. (2012, Section 3.2) analysed the second year of the maize trial with the first year as prior and all 11 bilinear terms. This vignette repeats their analysis and compares it with Table B.1 (the 38 linear effects), Table 2 (\(\sigma\), the singular values and the significant scores) and Figure 4 (the biplot).
2 The fit
The paper used two chains of 20,000 iterations. Here two chains of 5000 after a burn-in of 1000 suffice for the comparison.
All 38 posterior means agree with the published ones to within 0.035 posterior standard deviations, and the standard deviations agree too.
3.1 Why the posterior means are averages of the two years
With \(n_0 = B = 3\) the earlier trial weighs as much as the new one, so the posterior mean of every linear effect is the average of its two least-squares estimates:
both <-bys_ls_means(.data =rbind( cultivo2008[c("site", "rep", "entry", "y")] , cultivo2009[c("site", "rep", "entry", "y")] ) , .y = y , .gen = entry , .env = site )c(average_mu = both$grand_mean, published =3848)
average_mu published
3848.275 3848.000
4 Table 2: the residual standard deviation and the singular values
sig <-sqrt(unlist(fm$tau1))lam <- s[block =="lambda"]t2 <- data.table::data.table(parameter =c("sigma", sprintf("lambda%d", 1:11)) , published = pub$mean[match(x =c("sigma", sprintf("lambda[%d]", 1:11)), table = pub$parameter)] , here =round(c(mean(sig), lam$mean)) , sd_pub = pub$sd[match(x =c("sigma", sprintf("lambda[%d]", 1:11)), table = pub$parameter)] , sd_here =round(c(sd(sig), lam$sd)) )t2[, ratio :=round(x = here / published, digits =2)]t2
The residual standard deviation agrees: 1305.2 against 1311.00, 0.15 posterior standard deviations apart.
The singular values do not agree: here they are 1.33 to 2.07 times the published ones. By eq. 2.15 the posterior mean of \(\lambda_k\) is close to \((n_0\lambda_{0k} + B\,\va_k^\top\mY\vg_k)/(n_0 +
B)\). Here \(\lambda_{01} =\) 5306 is the first singular value of the first year. With equal weights, a published \(\lambda_1 = 2572\) below \(\lambda_{01}/2 =\) 2653 would need \(\va_1^\top\mY\vg_1\) to be negative on average, that is, the second-year data pulling the first term against itself. The leading published values are close to \(\lambda_{0k}/2\), as if the second year added nothing to them.
Two explanations were checked with the sampler of this package (dev/lambda-checks.R in the source) and neither gives the published values:
Score columns that do not move (no pair moves, as in the column-by-column sampler of the paper), started at the first year’s fit or at random orthonormal scores: \(\lambda_1\) is then 3613 or 3313, not 2572, because the environment scores still follow the data.
A prior mean of zero for the singular values with the same weight: \(\lambda_1\) is then 1203.
The linear effects and \(\sigma\), which come from the same draws, are reproduced; the sampler gives the same draws as an independent R version of eq. 2.15 and recovers the true values of simulated trials (validation vignette). The published singular values therefore disagree with the method as the paper states it, and the difference lies in the article, through a step it does not report or an error in its Table 2, rather than in the package.
5 Table 2: the significant scores
The paper lists the scores whose bivariate 0.95 HPD region excludes the origin: genotype 12 on the second term and environments 5, 9, 23 and 25 on the first.
On the first term the environment scores agree in sign and size with the published values (0.288, -0.277, 0.278 and 0.459 for environments 5, 9, 23 and 25; here 0.324, -0.257, 0.197, 0.505). Terms beyond the first are weakly identified when 11 terms of similar size are fitted: their order and orientation can differ between analyses.
Perez-Elizalde, S., Jarquin, D. and Crossa, J. (2012). A general Bayesian estimation method of linear-bilinear models applied to plant breeding trials with genotype x environment interaction. Journal of Agricultural, Biological, and Environmental Statistics, 17, 15–37. doi:10.1007/s13253-011-0063-9
Crossa, J., Perez-Elizalde, S., Jarquin, D., Cotes, J. M., Viele, K., Liu, G. and Cornelius, P. L. (2011). Bayesian estimation of the additive main effects and multiplicative interaction model. Crop Science, 51(4), 1458–1469. doi:10.2135/cropsci2010.06.0343
