--- title: "Unbalanced trials and environment-specific variances" 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{Unbalanced trials and environment-specific variances} %\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 Real trials lose plots, and some genotypes are not grown in every environment. The error variance also differs from one site to the next. This vignette shows how `bayes_ammi()` handles both: missing plots and empty cells by data augmentation, and a separate error variance for every environment with `.het = TRUE`. It ends with a comparison of the two variance models on `Maiz`. # Missing plots and empty cells ## The idea With $B$ the largest number of plots in a cell, a cell $(i, j)$ with $n_{ij} < B$ plots has $B - n_{ij}$ missing plots. A cell with $n_{ij} = 0$ is an empty cell. At every iteration the sampler draws the missing plots from the model, $$y^{\mathrm{mis}}_{ijl} \mid \vtheta \sim N(\eta_{ij}, \sigma^2_j), \qquad \eta_{ij} = \mu + \tau_i + \delta_j + \sum_{k=1}^{K}\lambda_k\alpha_{ik}\gamma_{jk},$$ and then draws the parameters from the conditional distributions of the balanced design (vignette 2) given the completed data. This is the data augmentation of Tanner and Wong (1987). The draws of the parameters come from their posterior given the observed plots only, so nothing is invented: a missing plot adds no information, and an empty cell is predicted from the row and column structure of the model. The missing plots are assumed to be missing at random (Little and Rubin 2019). The sampler needs only the sums, sums of squares and counts of the observed plots in every cell, so the completed data are never stored. ## An example From `Maiz` we remove 60 plots at random and every plot of three genotype-environment combinations: ```{r unbal} set.seed(seed = 21) mz <- Maiz[-sample(x = nrow(Maiz), size = 60), ] gone <- (mz$entry == 2 & mz$site == 5) | (mz$entry == 7 & mz$site == 11) | (mz$entry == 9 & mz$site == 18) mz <- mz[!gone, ] table(plots = table(factor(mz$entry):factor(mz$site))) ``` ```{r unbal-fit} fu <- bayes_ammi( .data = mz , .y = y , .gen = entry , .env = site , .rep = rep , .nIter = 3000 , .k = 2 , .burnin = 500 , .thin = 1 , .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 = 3 , .align = TRUE , .model = "AMMI" , .het = FALSE , .waic = TRUE ) fu ``` The print shows the number of empty cells. `bam_predict()` gives the posterior of every cell mean $\eta_{ij}$, including the empty cells: ```{r unbal-pred} pu <- bam_predict( fit = fu , .level = 0.95 ) pu$cells[n == 0] Maiz_cells <- data.table::as.data.table(x = Maiz)[, list(full = mean(y)), by = list( Gen = factor(entry) , Env = factor(site) )] merge(x = pu$cells[n == 0], y = Maiz_cells, by = c("Gen", "Env"))[, list( Gen , Env , mean , lower , upper , full )] ``` ```{r cover, include = FALSE} mm <- merge(x = pu$cells[n == 0], y = Maiz_cells, by = c("Gen", "Env")) inside <- sum(mm$lower <= mm$full & mm$full <= mm$upper) miss1 <- mm[!(lower <= full & full <= upper)] sdr <- mean(pu$cells[n == 0, sd]) / mean(pu$cells[n == 4, sd]) ``` The column `full` holds the cell means of the complete data, which the fit never saw. `r inside` of the 3 HPD intervals contain them.`r if (nrow(miss1)) paste0(" The miss, genotype ", miss1$Gen[1], " in environment ", miss1$Env[1], ", is a cell whose interaction is not carried by the two bilinear terms: the model predicts an empty cell only from the pattern of the rest of the trial.") else ""` The posterior standard deviations of the empty cells are on average `r round(sdr, 1)` times those of the cells with all four plots. ```{r unbal-plot, fig.height = 5} plot(x = pu, which = "cells") ``` # Environment-specific error variances ## The model With `.het = TRUE` the plots of environment $j$ have their own variance, $e_{ijl} \sim N(0, \sigma^2_j)$, and every $\sigma^2_j$ has the scaled inverse chi-squared prior of the common variance. Write $w_j = 1/\sigma^2_j$. The conditional distributions change as follows. * $\sigma^2_j$: scaled inverse chi-squared with $\nu + BG$ degrees of freedom and the residual sum of squares of environment $j$. * $\tau_i$: the data precision becomes $B\sum_j w_j$, the same for every genotype. * $\mu$ and $\delta_j$: no longer independent, so $(\mu, \mathbf{d})$ is drawn jointly from its multivariate normal conditional. * $\lambda_k$: mean and variance weighted by $w_j$. * $\va_k$: still von Mises-Fisher, because the quadratic term $\sum_j w_j\gamma_{jk}^2$ does not depend on $\va_k$. * $\vg_k$: the quadratic term $\sum_j w_j\gamma_{jk}^2$ now varies on the sphere. The conditional is a Bingham-von Mises-Fisher distribution (Hoff 2009), drawn one coordinate at a time in the eigenbasis of the quadratic form by slice sampling. Perez-Elizalde et al. (2012) give the general linear-bilinear model with these conditionals. ## The fit ```{r het} fh <- bayes_ammi( .data = Maiz , .y = y , .gen = entry , .env = site , .rep = rep , .nIter = 3000 , .k = 2 , .burnin = 500 , .thin = 1 , .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 = 4 , .align = TRUE , .model = "AMMI" , .het = TRUE , .waic = TRUE ) fh$summary[block == "sigma2", list(parameter, mean, lower, upper, rhat)] ``` ```{r het-range, include = FALSE} v <- fh$summary[block == "sigma2", mean] rr <- max(v) / min(v) ``` The posterior means range from `r format(round(min(v)), big.mark = ",")` to `r format(round(max(v)), big.mark = ",")`, a ratio of `r round(rr)`. One common variance cannot describe both extremes. ```{r het-plot, fig.height = 4.6} ggplot2::ggplot( data = fh$summary[block == "sigma2"] , mapping = ggplot2::aes(x = reorder(label, mean), y = mean, ymin = lower, ymax = upper) ) + ggplot2::geom_pointrange(colour = "#2c7fb8") + ggplot2::scale_y_continuous(breaks = scales::pretty_breaks(n = 6), labels = scales::label_comma()) + ggplot2::labs(x = "Environment", y = "Error variance (posterior mean, 95% HPD)") + ggplot2::theme_bw() ``` # Which variance model? ```{r het-common} fc <- bayes_ammi( .data = Maiz , .y = y , .gen = entry , .env = site , .rep = rep , .nIter = 3000 , .k = 2 , .burnin = 500 , .thin = 1 , .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 = 4 , .align = TRUE , .model = "AMMI" , .het = FALSE , .waic = TRUE ) cp <- bam_compare( common = fc , environment_var = fh ) cp ``` ```{r dw, include = FALSE} dw <- cp$table$dWAIC[2] se <- cp$table$dWAIC_se[2] l2 <- cp$table$lambda_last_lower ``` The environment-specific variances lower WAIC by `r round(dw)`, `r round(dw / se, 1)` times the standard error of the difference. They also sharpen the interaction: the lower HPD limit of $\lambda_2$ rises from `r round(l2[2])` to `r round(l2[1])`, because the plots of the noisy environments no longer blur the second term. The posterior predictive check of the common-variance model tells the same story through the ratio of the largest to the smallest environment variance: ```{r het-ppc} pp <- bam_ppc( fit = fc , .draws = 500 , .seed = 1 ) pp ``` ```{r ppcp, include = FALSE} pc <- pp$checks pv <- pc[statistic == "largest / smallest environment variance", p_value] ``` The observed ratio is `r round(pc[statistic == "largest / smallest environment variance", observed])`, against `r round(pc[statistic == "largest / smallest environment variance", mean_rep], 1)` on average in the replicated trials of the common-variance model (posterior predictive p-value `r format(pv, nsmall = 3)`). The other statistics are reproduced well. Vignette 9 returns to these tools. # 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 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 Hoff, P. D. (2009). Simulation of the matrix Bingham-von Mises-Fisher distribution, with applications to multivariate and relational data. *Journal of Computational and Graphical Statistics*, 18(2), 438--456. Tanner, M. A. and Wong, W. H. (1987). The calculation of posterior distributions by data augmentation. *Journal of the American Statistical Association*, 82(398), 528--540. Little, R. J. A. and Rubin, D. B. (2019). *Statistical Analysis with Missing Data*, 3rd edn. Wiley, Hoboken.