Chapter 12: Visualizing posteriors with bayesplot

1. Purpose

bayesplot provides ggplot2-based plotting for Markov chains and Bayesian workflow diagnostics. glmbayes does not depend on bayesplot (install it from CRAN if you want to run the code below). A former S3 pp_check() method for glmb lives in legacy_code/pp_check.glmb.R.

After fitting models in earlier chapters (Chapters 07–08 overview; Chapter 09 binomial illustration), visual summaries help assess mixing, posterior shape, and how well replicated data match the observed responses.

2. Minimal binomial fit (Menarche)

We reuse the Menarche grouped binomial illustration with default Prior_Setup() machinery.

data(menarche, package = "MASS")

Age2 <- menarche$Age - 13
Menarche_Model_Data <- data.frame(
  Menarche  = menarche$Menarche,
  Total     = menarche$Total,
  Age2      = Age2
)

ps <- Prior_Setup(
  cbind(Menarche, Total - Menarche) ~ Age2,
  family = binomial(link = "logit"),
  data   = Menarche_Model_Data
)

fit_logit <- glmb(
  cbind(Menarche, Total - Menarche) ~ Age2,
  family  = binomial(link = "logit"),
  pfamily = dNormal(mu = ps$mu, Sigma = ps$Sigma),
  data    = Menarche_Model_Data,
  n             = 800,
  use_parallel  = FALSE
)
coef_draws <- as.matrix(fit_logit$coefficients)

3. Coefficient trace and density overlays

Because glmb() returns i.i.d. draws, you can think of the coefficient matrix as one posterior sample stream. With bayesplot installed separately from CRAN, mcmc_trace() and mcmc_dens() give univariate views (functions that need multiple MCMC chains, such as mcmc_dens_overlay(), are not appropriate here).

4. Posterior predictive checks

For binomial models glmbayes works on proportions so y and yrep line up. The chunk below builds y_obs and y_rep for optional bayesplot::ppc_* plots (or the pp_check.glmb wrapper in legacy_code/).

## Observed success proportions (aligned with simulate.glmb for binomial)
y_obs <- Menarche_Model_Data$Menarche / Menarche_Model_Data$Total

pred_resp <- predict(fit_logit, type = "response")
nd <- max(1L, min(50L, nrow(pred_resp)))
pred_sub <- pred_resp[seq_len(nd), , drop = FALSE]

y_rep <- stats::simulate(
  fit_logit,
  pred = pred_sub,
  prior.weights = fit_logit$prior.weights
)

stopifnot(nrow(y_rep) == nd, ncol(y_rep) == length(y_obs))

Try other ppc_* helpers via the fun argument when bayesplot is installed (see ?bayesplot-helpers).

See also

• Chapter 08 — full glmb() walk-through and interpretation of printed summaries.
  
• ?simulate.glmb, ?predict.glmb — mechanics behind replicated responses; optional pp_check wrapper in legacy_code/pp_check.glmb.R.