Bayesian workflow in R: simulation, summaries, diagnostics, and reproducibility

NoteWhat you should be able to do after this chapter
  • use R as the orchestration layer for Bayesian workflows;
  • write reusable functions for conjugate updating and custom samplers;
  • compute posterior summaries from raw draws rather than only fitted-model printouts;
  • organize posterior predictive checks, diagnostics, and model comparison in R;
  • structure Bayesian analysis scripts for reproducibility and extension.

1 Why R matters in this course

R is not used here merely because it is common in statistics. It plays a specific technical role:

Even when the final sampler is Stan or JAX-based, R remains a highly effective environment for data preparation, posterior summarization, workflow diagnostics, and reporting.

2 Minimal workflow structure

A serious Bayesian R workflow usually has four layers:

  1. model definition: prior, likelihood, parameterization;
  2. computation: exact update, MCMC, Stan fit, or external engine;
  3. posterior object handling: draws, summaries, intervals, predictions;
  4. criticism: PPCs, ESS, \(\hat{R}\), LOO, sensitivity analysis.

The practical discipline is to keep those layers separable in code.

3 Conjugate functions as reusable building blocks

The fastest route to reproducible Bayesian code is to encapsulate exact model updates as functions.

For example, the Beta-Binomial update:

beta_binomial_update <- function(y, n, alpha, beta) {
  list(
    alpha_post = alpha + y,
    beta_post = beta + n - y,
    post_mean = (alpha + y) / (alpha + beta + n)
  )
}

This style matters because it makes exact calculations:

  • testable;
  • reusable in teaching materials;
  • comparable with approximate algorithms later.

See beta_binomial.R and conjugate_families.R.

4 Hand-coded MCMC in R

Writing one’s own sampler is educational because it forces clarity about:

  • target density;
  • proposal mechanism;
  • acceptance rule;
  • diagnostics.

The repository includes a random-walk Metropolis implementation in mh_sampler.R. A minimal structure is:

rw_metropolis <- function(init, n_iter, proposal_sd, log_target) {
  draws <- numeric(n_iter)
  draws[1] <- init
  accept <- 0L

  for (t in 2:n_iter) {
    proposal <- rnorm(1, draws[t - 1], proposal_sd)
    log_alpha <- log_target(proposal) - log_target(draws[t - 1])

    if (log(runif(1)) < log_alpha) {
      draws[t] <- proposal
      accept <- accept + 1L
    } else {
      draws[t] <- draws[t - 1]
    }
  }

  list(draws = draws, accept_rate = accept / (n_iter - 1))
}

This is obviously not industrial-strength software. Its value lies in transparency.

5 Working with posterior draws

Bayesian analysis should be draw-centric. Once posterior draws exist, many quantities follow by direct transformation.

5.1 Means and intervals

mean(draws)
quantile(draws, c(0.05, 0.5, 0.95))

5.2 Functionals

If a scientific quantity of interest is \(g(\theta)\), compute it on each draw:

g_draws <- exp(draws)
mean(g_draws)
quantile(g_draws, c(0.05, 0.95))

This draw-based approach avoids deriving new asymptotics for every transformed estimand.

6 Posterior predictive simulation in R

Posterior predictive checks become easy once draws are available:

y_rep <- replicate(length(mu_draws), {
  rnorm(n, mean = mu_draws[i], sd = sigma_draws[i])
})

Of course, in actual code one should index cleanly and vectorize where appropriate, but the principle is simple:

  • draw parameters from the posterior;
  • draw replicated data from the sampling model;
  • compare replicated and observed summaries.

See posterior_predictive_checks.R.

7 Diagnostics with posterior and bayesplot

When draws come from Stan or another modern MCMC backend, the posterior package standardizes storage and summaries.

Common workflow:

library(posterior)
library(bayesplot)

summarise_draws(draws_object)
mcmc_trace(draws_object, pars = c("beta[1]", "beta[2]"))
mcmc_acf(draws_object, pars = "beta[1]")

The advantage is consistency across projects and models, which matters more than the convenience.

8 Model comparison with loo

If pointwise log likelihood contributions are available, loo provides an efficient route to PSIS-LOO:

library(loo)
loo_result <- loo(log_lik_matrix)
print(loo_result)

This makes R the natural place to compare fitted models, aggregate predictive summaries, and present model criticism cleanly.

9 Organizing scripts cleanly

A strong repository layout usually separates:

  • labs/r/ for executable examples;
  • labs/stan/ for Stan programs;
  • chapters/ for course prose;
  • helper functions into compact reusable scripts when repetition begins.

Within an analysis, it is usually better to keep:

  • data generation or import explicit;
  • model functions pure when possible;
  • random seeds fixed in teaching examples;
  • plotting code separate from model code.

11 Existing labs in this repository

The repository now includes a broad R lab spine:

This is deliberate: R is the backbone of the course’s executable workflow.

12 Diagnostics and pitfalls

12.1 Treating fitted-model summaries as enough

Raw draws should be preserved and analyzed directly.

12.2 Mixing modeling, plotting, and data munging in one opaque script

This makes Bayesian work harder to debug and harder to teach.

12.3 Forgetting reproducible seeds

Stochastic examples should be reproducible by default in a course repository.

12.4 Hiding transformations

If a model quantity is standardized, transformed, or reparameterized, the code should make that explicit.

13 Exercises

  1. Write a reusable R function for a conjugate posterior update of your choice.
  2. Implement posterior predictive simulation from posterior draws and compare observed versus replicated summaries.
  3. Refactor an ad hoc analysis script into distinct sections for model definition, computation, summaries, and criticism.
  4. Use posterior draws to estimate a nonlinear functional and a credible interval for it.
  5. Explain why a draw-centric workflow is preferable to a summary-centric workflow.

14 Summary

R is the course’s operational glue:

  • exact examples are easy to encode;
  • custom algorithms are transparent;
  • posterior draws are easy to manipulate;
  • diagnostics and predictive checks are straightforward;
  • Quarto integration makes the workflow reproducible.

15 Further reading

The packages used here each have good documentation and it is worth reading it directly: Bürkner et al. (2024) for draws objects and the current convergence diagnostics, Gabry and Mahr (2024) for the diagnostic and posterior predictive graphics, Vehtari et al. (2024) for cross-validation and model comparison, and Gabry et al. (2024) for the Stan interface. R Core Team (2025) is the base language reference.

The methodological papers behind those packages are the ones to read for understanding rather than usage: Vehtari et al. (2021) for the diagnostics in posterior, Vehtari et al. (2017) for the estimators in loo, and Gabry et al. (2019) for the visualisations in bayesplot and the workflow they belong to. Gelman et al. (2020) ties the whole sequence together.

Back to top

References

Bürkner, Paul-Christian, Jonah Gabry, Aki Vehtari, et al. 2024. Posterior: Tools for Working with Posterior Distributions.
Gabry, Jonah, Rok Cešnovar, Andrew Johnson, Sean Bronder, Ben Goodrich, et al. 2024. Cmdstanr: Interface to ’CmdStan’.
Gabry, Jonah, and Tristan Mahr. 2024. Bayesplot: Plotting for Bayesian Models.
Gabry, Jonah, Daniel Simpson, Aki Vehtari, Michael Betancourt, and Andrew Gelman. 2019. “Visualization in Bayesian Workflow.” Journal of the Royal Statistical Society: Series A 182 (2): 389–402.
Gelman, Andrew, Aki Vehtari, Daniel Simpson, et al. 2020. “Bayesian Workflow.” arXiv Preprint arXiv:2011.01808.
R Core Team. 2025. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing.
Vehtari, Aki, Jonah Gabry, Yuling Yao, and Andrew Gelman. 2024. Loo: Efficient Leave-One-Out Cross-Validation and WAIC for Bayesian Models.
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.
Vehtari, Aki, Andrew Gelman, Daniel Simpson, Bob Carpenter, and Paul-Christian Bürkner. 2021. “Rank-Normalization, Folding, and Localization: An Improved \(\hat{R}\) for Assessing Convergence of MCMC.” Bayesian Analysis 16 (2): 667–718.