--- title: "Validation and corrections" 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{Validation and corrections} %\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.2, 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 Version 1.0.0 rewrites the Gibbs sampler of version 0.3.0. This vignette lists what changed and why, and checks the new sampler five ways: on its building blocks, against a plain R version of itself, on its information criteria, on simulated data with known parameters, and on the published analysis (vignette 3). # Corrections to version 0.3.0 | Function | Version 0.3.0 | Now | |---|---|---| | `bayes_ammi` | singular values computed once and never sampled (`lambdas1` constant) | drawn every iteration from their truncated normal conditional (eq. 7) | | `bayes_ammi` | `tau1` held the constant precision of the least-squares fit; $\sigma^2$ not returned | `tau1` holds the draws of $\sigma^2$ | | `bayes_ammi` | residuals for $\sigma^2$ against the environment of data row $j$ (with `Maiz` always site 1) | residuals of every cell | | `bayes_ammi` | within-cell variances summed in place of sums of squares | within-cell sum of squares | | `bayes_ammi` | covariance of $\vt$, $\vd$ = square root of the variance; $BE$ for $\vd$ | variances; $BG$ for $\vd$ | | `bayes_ammi` | scores drawn on the full Stiefel manifold, then orthogonalised to $\vone$ | von Mises-Fisher draws on the constrained spheres (eq. 8-9) | | `bayes_ammi` | all $\min(G, E) - 1$ terms; improper prior for $\sigma^2$; no burn-in, chains or seed | `.k` (default 2), proper priors, `.burnin`, `.thin`, `.chains`, `.seed` | | `biplots` | contours of the singular vectors, fixed bandwidth 0.2; "stable" by bounding box | HPD regions of $\lambda_k^{1/2}\alpha_{ik}$, normal reference bandwidths; region contains the origin | | `ge_model`, `ge_ammi` | `lme4::lmer()` for the error variance | closed form (`ge_model()` keeps `lme4`, now suggested) | | `Maiz` | help page of another data set | 720 rows: site, entry, folds, rep, y | | `bayes_ammi` | balanced trials only; one error variance; AMMI only; R loops | missing plots and cells, `.het`, `.model = "GGE"`, DIC and WAIC; compiled sampler | | new | none | `bam_predict`, `bam_winners`, `bam_superiority`, `bam_ranks`, `bam_stability`, `bam_compare`, `bam_ppc`, `bam_app` | | new | none | `bayes_stability`, `bam_stab_priors`, `bam_indices`, `bam_ammi_indices`, `bam_cite` | # The tables agree with version 0.3.0 `ge_mean()`, `ge_var()`, `g_eff()`, `e_eff()`, `ge_eff()`, `matrix_k()` and `orthnorm()` return the same numbers as before (test `test-helpers.R`). For example: ```{r tables} ge_eff( .data = Maiz , .y = y , .gen = entry , .env = site )$ge_svd$d[1:3] ``` # The building blocks The von Mises-Fisher sampler reproduces the mean resultant length $A_p(\kappa) = I_{p/2}(\kappa) / I_{p/2-1}(\kappa)$: ```{r vmf} set.seed(seed = 7) p <- 18 vapply( X = c(1, 10, 100) , FUN = function(kap) { X <- t(replicate(n = 3000, expr = bayesammi:::.bam_rvmf(v = c(kap, rep(0, p - 1))))) c( kappa = kap , simulated = mean(X[, 1]) , theory = besselI(x = kap, nu = p / 2, expon.scaled = TRUE) / besselI( x = kap , nu = p / 2 - 1 , expon.scaled = TRUE ) ) } , FUN.VALUE = numeric(length = 3) ) ``` # The compiled sampler and the R reference The sampler runs in C++. The package keeps a plain R version of the balanced, common-variance sampler, written directly from the formulas of vignette 2. Both take their random numbers from R in the same order, so from the same seed they must return the same draws: ```{r same} ns <- asNamespace(ns = "bayesammi") cd <- ns$.bam_prep( "check" , Maiz , quote(y) , quote(entry) , quote(site) , quote(rep) , environment() , allow_missing = TRUE ) ls <- ns$.bam_ls(Y = cd$Y, K = 2) pr <- 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 ) pr$s2 <- (cd$ssw + 4 * ls$rss) / 600 st <- ns$.bam_start(ls, pr$s2, 1) set.seed(seed = 5) t_r <- system.time(r <- ns$.bam_chain(cd$Y, 4, cd$ssw, 2, pr, st, list(alpha = ls$alpha), 500, 1:500)) set.seed(seed = 5) t_c <- system.time(cc <- ns$.bam_run_chain(cd, 2, FALSE, FALSE, pr, st, ls$alpha, 500, rep(1L, 500), TRUE, FALSE)) c(scalars = max(abs(cc$scalar - r$scalar)), alpha = max(abs(cc$alpha - r$alpha)), gamma = max(abs(cc$gamma - r$gamma))) c(R = t_r[["elapsed"]], compiled = t_c[["elapsed"]]) ``` The largest difference is at the level of rounding error, and the compiled sampler is `r round(t_r[["elapsed"]] / max(t_c[["elapsed"]], 0.001))` times faster here. # The stability sampler The sampler of `bayes_stability()` has a plain R version too, and the two must give the same draws for all three models: ```{r same-stab} vapply( X = 0:2 , FUN = function(code) { ls <- ns$.bam_stab_ls(cd, code) pr <- ns$.bam_stab_fill_priors( bam_stab_priors( mu_var = NULL , g_var = NULL , h_var = NULL , b_var = 100 , nu = 1 , s2 = NULL , nu_d = 1 , s2_d = NULL ) , cd , ls ) st <- ns$.bam_stab_start(ls, 1L, code) set.seed(seed = 9) r <- ns$.bam_stab_chain(cd$S, cd$Q, cd$N * 1, code, pr, st, 200, rep(1L, 200)) set.seed(seed = 9) cc <- ns$.bam_stab_cpp( cd$S , cd$Q , cd$N * 1 , code , pr , st , 200L , rep(1L, 200) , cd$data$y , cd$data$cell - 1L , FALSE ) max(abs(cc$mu - r$mu), abs(cc$g - r$g), abs(cc$b - r$b), abs(cc$d - r$d)) } , FUN.VALUE = numeric(length = 1) ) ``` The test `test-stability.R` also recovers known slopes and deviation variances from simulated trials with missing plots and an empty cell. # The information criteria WAIC computed during the run agrees with WAIC computed afterwards by **loo** from the pointwise log-likelihood of the kept draws (test `test-extensions.R`, to $10^{-8}$). DIC uses the plug-in deviance at the posterior means. # Simulated data Data are simulated from the model with known parameters, and the 95% HPD intervals should cover them: ```{r sim} set.seed(seed = 8) G <- 8 E <- 12 B <- 3 A <- qr.Q(qr(cbind(1, matrix(data = rnorm(G * 2), nrow = G))))[, 2:3] Gm <- qr.Q(qr(cbind(1, matrix(data = rnorm(E * 2), nrow = E))))[, 2:3] tau <- as.vector(scale(x = rnorm(n = G, mean = 0, sd = 6), scale = FALSE)) del <- as.vector(scale(x = rnorm(n = E, mean = 0, sd = 10), scale = FALSE)) M <- 100 + outer(X = tau, Y = del, FUN = "+") + A %*% (c(80, 40) * t(Gm)) sim <- expand.grid(rep = 1:B, gen = 1:G, env = 1:E) sim$y <- M[cbind(sim$gen, sim$env)] + rnorm(n = nrow(sim), mean = 0, sd = 5) fs <- bayes_ammi( .data = sim , .y = y , .gen = gen , .env = env , .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 = 9 , .align = TRUE , .model = "AMMI" , .het = FALSE , .waic = TRUE ) s <- fs$summary c(tau_covered = mean(s[block == "tau", lower] <= tau & tau <= s[block == "tau", upper]), delta_covered = mean(s[block == "delta", lower] <= del & del <= s[block == "delta", upper])) s[block %in% c("sigma2", "lambda"), list(parameter, truth = c(25, 80, 40), mean, lower, upper)] ``` # The published analysis Vignette 3 compares all 33 parameters of Table 1 and the 58 first scores of Table B1 of Crossa et al. (2011). Every posterior mean is within a quarter of a posterior standard deviation of the published one, and the test `test-crossa2011.R` repeats this comparison. ```{r ref} bam_reference( source = "Crossa" , evaluate = TRUE , level = 1 )[, list(quantity, value, computed, agrees)] ``` # 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 Wood, A. T. A. (1994). Simulation of the von Mises Fisher distribution. *Communications in Statistics - Simulation and Computation*, 23(1), 157--164. Mardia, K. V. and Jupp, P. E. (2000). *Directional Statistics*. Wiley, Chichester.