--- title: "Reproducing Crossa et al. (2011)" 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 Crossa et al. (2011)} %\VignetteEngine{quarto::html} %\VignetteEncoding{UTF-8} execute: echo: true warning: false message: false --- {{< include _bayesammi-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 = 4.6, dev = if (knitr::is_latex_output()) "cairo_pdf" else "ragg_png", dpi = 96, fig.process = webp) if (requireNamespace("bayesammi", quietly = TRUE)) library(bayesammi) else pkgload::load_all(export_all = FALSE, helpers = FALSE, quiet = TRUE) ``` # Abstract Crossa et al. (2011) fitted the Bayesian AMMI model with two bilinear terms to the maize trial `Maiz`. They ran two chains of 35000 iterations, discarded a burn-in of 5000 and thinned by 2. This vignette repeats that run and compares the posterior summaries with their Tables 1 and B1. The overall mean, the error variance, the main effects, the singular values and the scores of the first term agree within Monte Carlo error. The second scores differ, for a reason the paper does not discuss: the sign of a weak term. # The run of the paper ```{r fit} fit <- bayes_ammi( .data = Maiz , .y = y , .gen = entry , .env = site , .rep = rep , .nIter = 30000 , .k = 2 , .burnin = 5000 , .thin = 2 , .chains = 2 , .priors = bam_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 = 2011 , .align = TRUE , .model = "AMMI" , .het = FALSE , .waic = TRUE ) fit$settings$seconds ``` # Table 1 ```{r table1} pub <- utils::read.csv(file = system.file("extdata", "crossa2011.csv", package = "bayesammi")) s <- fit$summary t1 <- merge(x = pub[pub$table == "1", ], y = s[, list(parameter, mean, sd, lower, upper)], by = "parameter", suffixes = c(".paper", ".here"), sort = FALSE) t1$dev_sd <- round(x = (t1$mean.here - t1$mean.paper) / t1$sd.paper, digits = 3) num <- c("mean.paper", "mean.here", "sd.paper", "sd.here", "lower.paper", "lower.here", "upper.paper", "upper.here") t1[num] <- lapply(X = t1[num], FUN = signif, digits = 5) t1[, c("parameter", "mean.paper", "mean.here", "sd.paper", "sd.here", "dev_sd")] t1[, c("parameter", "lower.paper", "lower.here", "upper.paper", "upper.here")] ``` The column `dev_sd` is the difference of the posterior means in units of the posterior standard deviation; its largest absolute value is `r max(abs(t1$dev_sd))`. Both columns are Monte Carlo estimates, so they agree up to Monte Carlo error. ## The error variance The paper reports a prior scale of $S_0 = 10^8$ with $\nu = 1$, but its posterior of $\sigma^2$ (mean 763600) is what a vague prior gives. Taken literally, with $S^2_0 = 10^8$, the prior would add $10^8/721$ to the posterior mean: ```{r s2} fit_s2 <- bayes_ammi( .data = Maiz , .y = y , .gen = entry , .env = site , .rep = rep , .nIter = 3000 , .k = 2 , .burnin = 500 , .thin = 1 , .chains = 1 , .priors = bam_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 = 1e8 ) , .level = 0.95 , .seed = 1 , .align = TRUE , .model = "AMMI" , .het = FALSE , .waic = TRUE ) c(default = s[parameter == "sigma2", mean], s2_1e8 = fit_s2$summary[parameter == "sigma2", mean]) ``` The package therefore defaults to the least-squares residual mean square as the scale, which reproduces Table 1. # Table B1: the first term ```{r b1} b1 <- merge( x = pub[pub$table == "B1" & grepl(pattern = ",1\\]$", x = pub$parameter), ] , y = s[, list(parameter, mean, sd, lower, upper)] , by = "parameter" , suffixes = c(".paper", ".here") , sort = FALSE ) b1$dev_sd <- round(x = (b1$mean.here - b1$mean.paper) / b1$sd.paper, digits = 3) b1[c("mean.here", "sd.here")] <- lapply(X = b1[c("mean.here", "sd.here")], FUN = round, digits = 4) b1[, c("parameter", "mean.paper", "mean.here", "sd.paper", "sd.here", "dev_sd")] ``` The largest absolute deviation is `r max(abs(b1$dev_sd))` posterior standard deviations. # Table B1: the second term ```{r b2} b2 <- merge( x = pub[pub$table == "B1" & grepl(pattern = ",2\\]$", x = pub$parameter), ] , y = s[, list(parameter, mean, sd)] , by = "parameter" , suffixes = c(".paper", ".here") , sort = FALSE ) b2[c("mean.here", "sd.here")] <- lapply(X = b2[c("mean.here", "sd.here")], FUN = round, digits = 4) b2[, c("parameter", "mean.paper", "mean.here", "sd.paper", "sd.here")] ``` Here the means differ, mostly in sign, and the published means are smaller in absolute value. The second term is weak: the HPD interval of $\lambda_2$ reaches down to near zero, both in the paper (2.09) and here (`r round(s[parameter == "lambda[2]", lower], 2)`). A chain without sign alignment then moves between the two mirror images $(\va_2, \vg_2)$ and $(-\va_2, -\vg_2)$. Its posterior means are averages over both images and shrink towards zero: ```{r noalign} fit_free <- bayes_ammi( .data = Maiz , .y = y , .gen = entry , .env = site , .rep = rep , .nIter = 30000 , .k = 2 , .burnin = 5000 , .thin = 2 , .chains = 1 , .priors = bam_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 = 2011 , .align = FALSE , .model = "AMMI" , .het = FALSE , .waic = TRUE ) d <- fit_free$draws[[1]][, "alpha[8,2]"] c(share_positive = mean(d > 0), mean = mean(d), paper = pub$mean[pub$parameter == "alpha[8,2]"]) ``` ```{r share, include = FALSE} am <- s[parameter == "alpha[8,2]", mean] pp <- pub$mean[pub$parameter == "alpha[8,2]"] ``` The published value, `r pp`, is smaller in absolute value than the aligned posterior mean, `r round(am, 3)`, and has the opposite sign. A chain that spends a share $q$ of its run in the image opposite to the least-squares orientation has a mean of about $(1 - 2q)$ times the aligned one. The published value corresponds to $q \approx$ `r round((1 - pp / am) / 2, 2)`: the published chain spent most of its run in the other image. The aligned summaries in this package describe one orientation of the term, and the scores of genotypes then keep their meaning relative to each other. # Figure 1: traces and histograms of the singular values ```{r fig1, fig.height = 5} plot(x = fit, which = "trace", pars = c("lambda[1]", "lambda[2]")) plot(x = fit, which = "density", pars = c("lambda[1]", "lambda[2]")) ``` # Figures 2 to 4: biplots ```{r fig2, fig.height = 6} plot(x = fit, which = "biplot") ``` The coordinates are $\lambda_k^{1/2}\alpha_{ik}$ and $\lambda_k^{1/2}\gamma_{jk}$, as the paper defines them. The axes of its figures run from -1 to 1, so they show a rescaled version; the relative positions are the same. ```{r fig34, fig.width = 9, fig.height = 5.5} bp <- biplots( model = fit , burnin = 0 , thin = 0.25 , pb = c(0.05, 0.01) , plot_stable = TRUE , plot_unstable = TRUE , ncolors = 5 , .which = "both" , .grid = 100 ) bp$biplot bp$biplot_data[axis1_significant == TRUE, list(byplot, name, x, axis1_lower, axis1_upper, group)] ``` On the first term, genotypes 4, 5, 6 and 9 form a positive group and genotypes 1 and 8 a negative group, as in the paper's Figure 3. Sites S8 (negative) and S11, S13 and S18 (positive) stand out, as in its Figure 4. With the intervals of Table B1, genotypes 2 and 3 and some further sites also exclude zero at 95%; the paper's figures show a subset of them. # References 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