Skip to contents

Bayesian inference

Source: vignettes/mcmc.qmd

Overview

This vignette demonstrates how to conduct Bayesian inference for the CCSBT stock assessment model using the SparseNUTS package and the sparse inverse Hessian generated by RTMB (Monnahan et al. 2026) (https://github.com/noaa-afsc/SparseNUTS). We illustrate the setup, sampling procedure, and posterior diagnostics. If you have not done so already, you will need to install the packages below:

remotes::install_github("janoleko/RTMBdist")
remotes::install_github("andrjohns/StanEstimators")
remotes::install_github("noaa-afsc/SparseNUTS")

The model vignette finishes with an optimised sbt_fit. Save that object there, then read it here so the MCMC analysis continues from the same fit.

library(sbt)
library(tidyverse)
library(reshape2)
library(bayesplot)
library(SparseNUTS)

theme_set(theme_bw())

article_extdata <- function(file) {
  candidates <- c(
    file.path("inst", "extdata", file),
    file.path("..", "inst", "extdata", file),
    system.file("extdata", file, package = "sbt")
  )
  candidates <- candidates[nzchar(candidates) & file.exists(candidates)]
  if (!length(candidates)) {
    stop("Could not find extdata file: ", file, call. = FALSE)
  }
  candidates[[1]]
}

Markov chain Monte Carlo

Read the fit saved at the end of the model vignette. sbt_mcmc() uses its MLE for initialisation, runs SparseNUTS, and stores the posterior and diagnostics in the same object.

fit <- sbt_fit_read("base-model.sbt.rds")
fit <- sbt_mcmc(
  fit,
  num_samples = 500L,
  num_warmup = 1000L,
  chains = 4L,
  cores = 4L,
  control = list(adapt_delta = 0.999)
)
sbt_fit_save(fit, "base-model.sbt.rds", overwrite = TRUE)

The diagnostic sections use a bundled full posterior so that the sampler is not rerun whenever the article is rendered.

The demonstration MCMC was split into 4 chains, each with a warm-up of 1000 iterations followed by 500 retained iterations. This resulted in 2000 posterior samples for 1553 parameters. This MCMC took 49 minutes to complete.

If sampling is unnecessarily slow, a lower adapt_delta can reduce runtime, but it may increase divergent transitions. If divergences persist, increase warmup and/or adapt_delta and recheck all diagnostics. Consider thinning large posterior outputs only to reduce storage or plotting cost; thinning does not repair poor convergence.

Diagnostics

Sampler parameters are shown in Figure 1.

plot_sampler_params(fit = mcmc, plot = TRUE)
Figure 1: Sampler parameters for the bundled full posterior.

A pairs plot for the five slowest-mixing parameters is shown in Figure 2. Adjust pars in pairs_rtmb() to inspect specific parameter blocks.

pairs_rtmb(fit = mcmc, order = "slow", pars = 1:5)
Figure 2: Pairs plot for the five slowest-mixing parameters in the bundled full posterior. The MLE point and asymptotic SE ellipse are overlaid when available. Divergent transitions are identified as separate points if present.

Bayesian and frequentist uncertainty estimates are compared in Figure 3.

plot_uncertainties(fit = mcmc, log = TRUE, plot = TRUE)
Figure 3: Comparison of Bayesian and frequentist uncertainty estimates in the bundled full posterior.

Marginal distributions are compared in Figure 4.

plot_marginals(fit = mcmc, pars = 1:6)
Figure 4: Marginal distributions in the bundled full posterior compared with the approximation around the mode.

MCMC trace plots are shown in Figure 5.

post <- as.data.frame(mcmc)
pars <- mcmc$par_names[1:6]
mcmc_trace(x = post, pars = pars)
Figure 5: Trace plots for the bundled full posterior.

The NUTS energy diagnostic is shown in Figure 6, with additional diagnostics in Figure 7, Figure 8, Figure 9, Figure 10, and Figure 11. The energy plot compares energy__ with the change in energy__ (Betancourt 2017).

np <- extract_sampler_params(fit = mcmc) |>
  pivot_longer(-c(chain, iteration), names_to = "Parameter") |>
  select(Iteration = iteration, Parameter, Value = value, Chain = chain) |>
  mutate(Parameter = factor(Parameter),
         Iteration = as.integer(Iteration),
         Chain = as.integer(Chain)) |>
  as.data.frame()

mcmc_nuts_energy(x = np)
`stat_bin()` using `bins = 30`. Pick better value `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value `binwidth`.
Figure 6: NUTS energy diagnostic. 
lp <- extract_samples(fit = mcmc, inc_lp = TRUE, as.list = TRUE) |>
  melt(value.name = "Value") |>
  filter(Var2 == "lp__") |>
  mutate(Chain = as.integer(L1), Iteration = as.integer(Var1)) |>
  select(Chain, Iteration, Value)
mcmc_nuts_acceptance(x = np, lp = lp)
`stat_bin()` using `bins = 30`. Pick better value `binwidth`.
Figure 7: NUTS acceptance. 
mcmc_nuts_divergence(x = np, lp = lp)
Warning: Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Groups with fewer than two datapoints have been dropped.
ℹ Set `drop = FALSE` to consider such groups for position adjustment purposes.
Figure 8: NUTS divergence. 
mcmc_nuts_stepsize(x = np, lp = lp)
Figure 9: NUTS stepsize. 
mcmc_nuts_treedepth(x = np, lp = lp)
Figure 10: NUTS treedepth.
Figure 11: NUTS treedepth.

Leave-One-Out Information Criterion

Leave-one-out information criterion (LOO-IC) estimates out-of-sample predictive performance from a Bayesian model (Vehtari et al. 2017; Yao et al. 2018; Magnusson et al. 2019). Exact leave-one-out cross-validation would refit the model after omitting each data point. get_loo() instead uses Pareto-smoothed importance sampling (PSIS) through the loo package to approximate those leave-one-out distributions from the existing posterior draws (Vehtari et al. 2024).

Lower LOO-IC indicates better estimated predictive performance only when comparing models fitted to the same observations and likelihood partition. The pointwise Pareto-() diagnostics must also be checked: influential or unreliable points can invalidate a simple ranking. LOO-IC is therefore one model-comparison diagnostic, not a universal test of model adequacy or a substitute for biological and sampler diagnostics.

Composition rows with zero input effective N remain in the model data and prediction/output structures, but get_loo() excludes them from the pointwise LOO matrix. plot_loo() also omits non-finite Pareto shape estimates rather than displaying them at an artificial infinity axis break.

looic <- get_loo(fit)
looic

Computed from 2000 by 5464 log-likelihood matrix.

         Estimate    SE
elpd_loo  -2074.6 117.1
p_loo       282.6  24.2
looic      4149.1 234.1
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 1.5]).

Pareto k diagnostic values:
                         Count Pct.    Min. ESS
(-Inf, 0.7]   (good)     5357  98.0%   75
   (0.7, 1]   (bad)        52   1.0%   <NA>
   (1, Inf)   (very bad)   55   1.0%   <NA>
See help('pareto-k-diagnostic') for details.
plot_loo(x = looic, data = data, exclude = "pop")
Figure 12: Pareto-smoothed importance-sampling (PSIS) diagnostic plot. The plot displays the estimated Pareto shape parameter (k̂) for each finite pointwise diagnostic across data sets. AF and LF points are grouped by fishery so influential points can be traced to fisheries such as LL1. Composition rows with zero input effective N remain in model I/O but are excluded from the LOO calculation, and non-finite Pareto shape estimates are omitted from the plot. Dashed lines at 0.5, 0.7, and 1.0 provide fixed visual references; the loo package uses the current sample-size-dependent reliability threshold described by Vehtari et al. (2024).

References

Betancourt, Michael. 2017. A Conceptual Introduction to Hamiltonian Monte Carlo. https://arxiv.org/abs/1701.02434.
Magnusson, Måns, Michael Andersen, Johan Jonasson, and Aki Vehtari. 2019. “Bayesian Leave-One-Out Cross-Validation for Large Data.” In Proceedings of the 36th International Conference on Machine Learning, edited by Kamalika Chaudhuri and Ruslan Salakhutdinov, vol. 97. Proceedings of Machine Learning Research. PMLR. https://proceedings.mlr.press/v97/magnusson19a.html.
Monnahan, Cole C., Kasper Kristensen, James T. Thorson, and Bob Carpenter. 2026. Leveraging Sparsity to Improve No-u-Turn Sampling Efficiency for Hierarchical Bayesian Models. https://arxiv.org/abs/2603.02437.
Vehtari, Aki, Andrew Gelman, and Jonah Gabry. 2017. “Practical Bayesian Model Evaluation Using Leave-One-Out Cross-Validation and WAIC.” Statistics and Computing 27 (5): 1413–32. https://doi.org/10.1007/s11222-016-9696-4.
Vehtari, Aki, Daniel Simpson, Andrew Gelman, Yuling Yao, and Jonah Gabry. 2024. “Pareto Smoothed Importance Sampling.” Journal of Machine Learning Research 25 (72): 1–58. https://jmlr.org/papers/v25/19-556.html.
Yao, Yuling, Aki Vehtari, Daniel Simpson, and Andrew Gelman. 2018. “Using Stacking to Average Bayesian Predictive Distributions (with Discussion).” Bayesian Analysis 13 (3): 917–1007. https://doi.org/10.1214/17-BA1091.