--- title: "Priors and the number of bilinear terms" 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{Priors and the number of bilinear terms} %\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.4, 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 The defaults of `bam_priors()` are vague: the data dominate the posterior. This vignette shows how to set informative priors from earlier trials and how the posterior of the last singular value guides the choice of the number of bilinear terms $K$. # Vague priors ```{r vague} f0 <- 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 = NULL ) , .level = 0.95 , .seed = 6 , .align = TRUE , .model = "AMMI" , .het = FALSE , .waic = TRUE ) f0$summary[block %in% c("lambda", "sigma2"), list(parameter, mean, sd, lower, upper)] ``` # Informative priors Suppose earlier trials in the region give an error variance of about $8 \times 10^5$, with the information of 50 plots, and a first singular value near 5000: ```{r informative} f1 <- 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 = 5000 , lambda_var = 250^2 , nu = 50 , s2 = 8e5 ) , .level = 0.95 , .seed = 6 , .align = TRUE , .model = "AMMI" , .het = FALSE , .waic = TRUE ) pm <- rbind(vague = f0$summary[block %in% c("lambda", "sigma2"), mean], informative = f1$summary[block %in% c("lambda", "sigma2"), mean]) colnames(pm) <- f0$summary[block %in% c("lambda", "sigma2"), parameter] round(pm) ``` The prior on $\lambda_k$ applies to every term. An informative prior suited to the first term therefore also pulls the second towards 5000 (here from `r round(pm[1, 3])` to `r round(pm[2, 3])`), a reason to keep `lambda_var` large unless the prior information is about all terms. # How many terms? The paper fitted $K = 2$ terms. A term that the data do not support has a singular value whose posterior reaches down to zero: ```{r k} fits <- lapply( X = 1:3 , FUN = function(K) { bayes_ammi( .data = Maiz , .y = y , .gen = entry , .env = site , .rep = rep , .nIter = 3000 , .k = K , .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 = NULL ) , .level = 0.95 , .seed = 6 , .align = TRUE , .model = "AMMI" , .het = FALSE , .waic = TRUE ) } ) data.table::rbindlist(l = lapply(X = fits, FUN = function(f) { f$summary[block == "lambda", list(K = f$settings$k, parameter, mean, sd, lower, upper)] })) ``` ```{r klow, include = FALSE} lw <- vapply(X = fits, FUN = function(f) { l <- f$summary[block == "lambda"] l$lower[nrow(l)] }, FUN.VALUE = numeric(length = 1)) ``` With $K = 2$ the HPD interval of $\lambda_2$ reaches down to `r signif(lw[2], 2)`, and with $K = 3$ that of $\lambda_3$ reaches down to `r signif(lw[3], 2)`. In both cases the last term is barely supported, whereas the interval of $\lambda_1$ stays far from zero. With a third term in the model, $\lambda_2$ moves away from zero because the order $\lambda_2 \ge \lambda_3$ gives the third term the weakest pattern. The least-squares singular values (5923, 3070, 2552) do not show any of this, because they are biased upwards. Cornelius et al. (1996) found two terms significant for these data, and Liu (2001, unpublished) suggested two or three. # The posterior of the error variance with more terms ```{r s2} vapply( X = fits , FUN = function(f) f$summary[parameter == "sigma2", mean] , FUN.VALUE = numeric(length = 1) ) ``` Each further term takes $G + E - 1 - 2k$ degrees of freedom from the residual, so the error variance decreases only slowly once the terms are weak. # Information criteria WAIC compares the three fits by their estimated out-of-sample fit (vignette 9): ```{r waic} names(fits) <- paste0("K", 1:3) cp <- do.call(what = bam_compare, args = fits) cp ``` ```{r waicw, include = FALSE} best <- cp$table$model[1] ``` WAIC prefers `r best`. `dWAIC` is the difference from the best model and `dWAIC_se` its standard error, computed from the plot-by-plot differences between the two models; a difference of less than about two of its standard errors does not separate the models. # 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 Cornelius, P. L., Crossa, J. and Seyedsadr, M. S. (1996). Statistical tests and estimators of multiplicative models for genotype-by-environment interaction. In M. S. Kang and H. G. Gauch (eds), *Genotype-by-Environment Interaction*. CRC Press, Boca Raton, 199--234. Cornelius, P. L. and Crossa, J. (1999). Prediction assessment of shrinkage estimators of multiplicative models for multi-environment cultivar trials. *Crop Science*, 39(4), 998--1009. Hobert, J. P. and Casella, G. (1996). The effect of improper priors on Gibbs sampling in hierarchical linear mixed models. *Journal of the American Statistical Association*, 91(436), 1461--1473. Watanabe, S. (2010). Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. *Journal of Machine Learning Research*, 11, 3571--3594. Gelman, A., Hwang, J. and Vehtari, A. (2014). Understanding predictive information criteria for Bayesian models. *Statistics and Computing*, 24(6), 997--1016.