--- title: "Reproducing Perez-Elizalde et al. (2012)" author: "Muhammad Yaseen" date: last-modified format: html: self-contained: true number-sections: true toc: true toc-depth: 3 highlight-style: pygments code-tools: true html-math-method: mathjax vignette: > %\VignetteIndexEntry{Reproducing Perez-Elizalde et al. (2012)} %\VignetteEngine{quarto::html} %\VignetteEncoding{UTF-8} execute: echo: true warning: false message: false --- {{< include _baystability-macros.txt >}} ```{r setup, include = FALSE} # HTML figures are drawn as PNG and kept as WebP: lossless (the same pixels in about half the space of PNG), or # lossy at quality 90 when that is less than half the size, as for traces and densities of thousands of draws. The # PDF figures stay vector graphics. Any failure of magick (not installed, no WebP, or a broken installation whose # system libraries are missing) keeps the PNG webp <- function(path) { if (knitr::is_latex_output() || !grepl(pattern = "[.]png$", x = path)) return(path) tryCatch( { if (!requireNamespace("magick", quietly = TRUE) || !isTRUE(magick::magick_config()$webp)) return(path) img <- magick::image_read(path) out <- sub(pattern = "[.]png$", replacement = ".webp", x = path) tmp <- tempfile(fileext = ".webp") magick::image_write( image = img , path = out , format = "webp" , quality = 100 , defines = c(`webp:lossless` = "true") ) magick::image_write(image = img, path = tmp, format = "webp", quality = 90) if (file.size(tmp) < file.size(out) / 2) file.copy(from = tmp, to = out, overwrite = TRUE) unlink(c(path, tmp)) out } , error = function(e) path ) } knitr::opts_chunk$set(collapse = FALSE, comment = "", fig.width = 7, fig.height = 6, dev = if (knitr::is_latex_output()) "cairo_pdf" else "ragg_png", dpi = 96, fig.process = webp) if (requireNamespace("baystability", quietly = TRUE)) library(baystability) else pkgload::load_all(export_all = FALSE, helpers = FALSE, quiet = TRUE) ``` # Abstract 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). # 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. ```{r fit} fm <- bayes_ammi( .data = cultivo2009 , .y = y , .gen = entry , .env = site , .rep = rep , .nIter = 5000 , .prior_data = cultivo2008 , .n0 = NULL , .model = "AMMI" , .k = NULL , .burnin = 1000 , .thin = 1 , .chains = 2 , .priors = bys_priors( mu_mean = 0 , mu_var = 1e+8 , tau_mean = 0 , tau_var = 1e+8 , delta_mean = 0 , delta_var = 1e+8 , lambda_mean = 0 , lambda_var = 1e+8 , nu = 1 , s2 = NULL ) , .level = 0.95 , .seed = 2012 , .align = TRUE , .waic = TRUE ) ``` ```{r conv} s <- fm$summary c(max_rhat = max(s$rhat), min_ess_linear = min(s[block %in% c("mu", "tau", "delta"), ess])) ``` # Table B.1: the linear effects ```{r b1} pub <- utils::read.csv(file = system.file("extdata", "perez2012.csv", package = "baystability")) pos <- c("mu", sprintf("tau[%d]", 1:12), sprintf("delta[%d]", 1:25)) lin <- s[block %in% c("mu", "tau", "delta")] b1 <- data.table::data.table( parameter = c("mu", sprintf("alpha%d", 1:12), sprintf("beta%d", 1:25)) , published = pub$mean[match(x = pos, table = pub$parameter)] , here = round(x = lin$mean, digits = 1) , sd_pub = pub$sd[match(x = pos, table = pub$parameter)] , sd_here = round(x = lin$sd, digits = 2) ) b1[, diff_in_sd := round(x = (here - published) / sd_pub, digits = 3)] b1 ``` ```{r b1f, include = FALSE} mz <- max(abs(b1$diff_in_sd)) ``` All 38 posterior means agree with the published ones to within `r mz` posterior standard deviations, and the standard deviations agree too. ## 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: ```{r avg} 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) ``` # Table 2: the residual standard deviation and the singular values ```{r t2} 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 ``` ```{r t2f, include = FALSE} zs <- (mean(sig) - 1311) / 39.13 l01 <- fm$prior$fit$lambda[1] ``` The residual standard deviation agrees: `r round(mean(sig), 1)` against 1311.00, `r round(abs(zs), 2)` posterior standard deviations apart. The singular values do not agree: here they are `r min(t2$ratio[-1])` to `r max(t2$ratio[-1])` 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} =$ `r round(l01)` is the first singular value of the first year. With equal weights, a published $\lambda_1 = 2572$ below $\lambda_{01}/2 =$ `r round(l01 / 2)` 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. # 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. ```{r scores} sc <- s[block %in% c("alpha", "gamma") & axis %in% 1:2] sc[(block == "alpha" & label == levels(factor(cultivo2009$entry))[12]) | (block == "gamma" & label %in% c("5", "9", "23", "25")), list( parameter , axis , mean , sd , lower , upper )] ``` ```{r pubsc, include = FALSE} g <- s[block == "gamma" & axis == 1 & label %in% c("5", "9", "23", "25"), mean] ``` 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 `r paste(round(g, 3), collapse = ", ")`). Terms beyond the first are weakly identified when 11 terms of similar size are fitted: their order and orientation can differ between analyses. # Figure 4: the biplot ```{r fig4, fig.width = 9} bp <- bayes_ammi_biplot( fit = fm , burnin = 0 , thin = 1 , pb = c(0.10, 0.05) , plot_stable = FALSE , plot_unstable = TRUE , ncolors = 5 , .which = "both" , .grid = 100 ) bp$biplot ``` # References 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