06 Model Checking and Model Comparison

How to decide whether a Bayesian model is credible

1 Learning Outcomes

By the end of this module, participants should be able to:

Check basic MCMC diagnostics.
Use posterior predictive checks.
Compare models using LOO.
Separate computational problems from substantive model problems.

2 Two Kinds of Checking

Bayesian model checking is doing two jobs at once, and conflating them is the most common source of confusion when participants encounter their first warning message from Stan. The first job is computational: did the sampler explore the posterior properly, or did it fail to mix, get stuck, or signal a problem with the model parameterisation? The second job is substantive: even if the computation succeeded, does the fitted model describe the actual data well enough to support the research claim?

A model can fail either check independently. A computational failure produces unreliable posterior summaries no matter how good the model is conceptually; the summaries describe something other than the intended posterior. A substantive failure means the model is fitted correctly to a specification that does not match the data, so the posterior is reliable as a description of an inappropriate model. Both failures need to be diagnosed. Reporting one set of checks without the other leaves the reader unable to judge whether to trust the result.

3 Computational Diagnostics

The standard computational diagnostics for Hamiltonian Monte Carlo, the algorithm Stan uses, are R-hat, effective sample size, divergent transitions, and trace plots. Each one is checking a different aspect of whether the sampler did its job.

R-hat is the ratio of between-chain to within-chain variance for each parameter. If the chains are exploring the same distribution, that ratio should be very close to 1. An R-hat of 1.00 across all parameters is the standard target. An R-hat above 1.01 is a sign that the chains have not converged to the same place; one or more chains is exploring a different region of the posterior. The fix is usually to run longer chains, but persistent R-hat problems indicate that the model is poorly identified and the posterior has multiple modes the sampler cannot reconcile.

Effective sample size (ESS) is the number of effectively independent samples in the chain after accounting for autocorrelation. A nominal 4,000 posterior draws can have an ESS of only a few hundred if the chains move slowly through the posterior. Low ESS means posterior summaries, especially extreme quantiles like the 2.5% and 97.5% bounds of a credible interval, are noisy. The standard recommendation is at least 400 effective samples for any quantity the researcher will report.

Divergent transitions are a Stan-specific warning that indicates the sampler encountered a region of the posterior with very high curvature that it could not navigate. Even one or two divergences should not be ignored, especially in hierarchical models, where they often signal that the group-level standard deviation is poorly estimated. The usual fix is to increase the adapt_delta argument in brm() to a value like 0.99, or to reparameterise the model using a non-centred parameterisation.

Trace plots show the chain values across iterations, typically with four chains overlaid. A well-mixed model produces traces that look like a “hairy caterpillar”: the chains overlap, move around the same central region, and show no obvious trend or stickiness. Traces that diverge, drift, or get stuck in different regions are visible signs of the problems that R-hat and ESS quantify.

summary(bayes_model)
plot(bayes_model)

The summary() call returns R-hat and ESS for every parameter alongside the posterior summaries. The plot() call returns trace plots and density plots for each parameter, four chains overlaid. Reading these two outputs together is the routine first step after any brms fit.

4 Posterior Predictive Checks

Computational diagnostics tell the researcher whether the sampler explored the intended posterior. They cannot tell the researcher whether the intended posterior is a posterior worth interpreting. The substantive check for that is the posterior predictive check, which asks:

If the model were true, would simulated data look like the observed data?

The mechanic is straightforward. The fitted model is used to generate many simulated datasets, each one a plausible outcome under the model and the posterior. Those simulated datasets are then compared with the actual data on dimensions the researcher cares about: the shape of the distribution, the mean, the standard deviation, the proportion above a threshold, the variance by subgroup. Systematic disagreement between the simulations and the data is a sign the model is missing something important.

library(bayesplot)

pp_check(bayes_model, ndraws = 50)

The default pp_check() overlays the density of the observed outcome on the densities of 50 simulated outcomes. If the observed density sits naturally inside the cloud of simulated ones, the model is generating data that looks like the data the researcher actually has. If the observed density is in the tails of the simulated cloud, or has a feature (a second mode, a different skew, a hard upper bound) that the simulations do not reproduce, the model is misspecified in a way the analyst now has to investigate.

Targeted checks examine specific summary statistics:

pp_check(bayes_model, type = "dens_overlay", ndraws = 50)
pp_check(bayes_model, type = "stat", stat = "mean")
pp_check(bayes_model, type = "stat", stat = "sd")

The type = "stat" variant computes a summary statistic on each simulated dataset and on the observed data, then plots the distribution of simulated values with the observed value marked. A model that produces a posterior predictive distribution of means that does not contain the observed mean is a model that does not match the central tendency of the data. The same check on the standard deviation catches models that under- or over-estimate the spread, which is a common failure mode for models with the wrong likelihood family.

5 What Poor Fit Can Mean

A failed posterior predictive check is a diagnosis, not a verdict. It points to a mismatch between the model and the data without specifying which part of the model is wrong. The candidate explanations are usually one or more of the following: the wrong likelihood family for the outcome, a missing predictor whose effect is producing the unexplained structure, a non-linear relationship that the model is forcing into a linear shape, an unmodeled grouping that creates clusters the model does not see, outliers or heavy tails that the chosen likelihood cannot accommodate, priors that are constraining the parameters away from values the data supports, or measurement problems in the outcome itself.

Diagnosing which explanation applies requires looking at the specific way the model fails. If the simulated standard deviations cluster well below the observed standard deviation, the model is under-dispersed; a Negative Binomial in place of a Poisson is a common fix for count data with this symptom. If the simulated distributions miss a second mode in the observed data, the model is missing a categorical predictor that separates two groups. If individual residuals are systematically more extreme than the model expects, a likelihood with heavier tails (Student-t in place of Normal) may help.

The point is not to chase a perfect model. The point is to know whether the model is adequate for the substantive claim the research wants to make. A model that gets the central tendency right but the tails wrong may be fine for reporting a mean effect, but poor for predicting unusual cases. The check tells the researcher which kinds of claim the model supports.

6 Model Comparison with LOO

Posterior predictive checks evaluate a single model against the observed data. Model comparison asks a related but different question: among several candidate models, which one is expected to predict new data best? Leave-one-out cross-validation (LOO) approximates this by estimating how well the model would predict each observation if that observation had not been used to fit the model. Computationally, this is done efficiently using importance sampling on the posterior, without refitting the model for each held-out point.

library(loo)

loo_simple <- loo(bayes_model)
loo_hier <- loo(hier_model)

loo_compare(loo_simple, loo_hier)

The output of loo_compare() reports the difference in expected log predictive density between models, along with the standard error of that difference. The convention is to interpret a difference larger than 2 standard errors as substantively meaningful evidence that one model predicts better than the other. A smaller difference is consistent with the two models predicting comparably well, in which case the researcher can choose between them on grounds other than predictive performance (parsimony, interpretability, alignment with theory).

LOO answers a predictive question, not a universal-truth question. A model with better LOO performance is not automatically more interpretable or more theoretically defensible. A hierarchical model can predict only slightly better than a non-hierarchical one and still be the right choice because it correctly represents the dependence structure. The reverse can also hold: a richer model can predict better and still be the wrong choice if the additional parameters are not interpretable as the substantive quantities the research is meant to estimate.

7 Exercise

For a fitted model:

Check R-hat and effective sample size.
Inspect trace plots.
Run a posterior predictive density overlay.
Compare two models with LOO.
Write a short model-checking note.

8 Model-Checking Checklist

Check	Question	Result
R-hat	Are values close to 1.00?
ESS	Are effective sample sizes adequate?
Divergences	Are there divergent transitions?
Trace plots	Do chains mix well?
Posterior predictive check	Do simulated data resemble observed data?
Sensitivity	Do conclusions depend heavily on a questionable prior?
Substantive adequacy	Is the model good enough for the claim?

9 Takeaway

Bayesian analysis is not complete when the model finishes running. The output of a successful brm() call is the beginning of a checking process, not the end of an analysis. A credible workflow runs the computational diagnostics first to confirm the sampler explored the intended posterior, runs the posterior predictive checks to confirm the model describes the observed data adequately for the claim being made, and runs sensitivity checks on the prior choices that materially influence the conclusion. Only then are the posterior summaries safe to report. The reporting notebook covers what to actually write down.

--- title: "06 Model Checking and Model Comparison" subtitle: "How to decide whether a Bayesian model is credible" format: html: toc: true number-sections: true execute: warning: false message: false --- ```{=html} <nav class="module-route" aria-label="Notebook route"> <a href="01_bayesian_reasoning.html">01<br>Reasoning</a> <a href="02_priors.html">02<br>Priors</a> <a href="03_bayesian_regression.html">03<br>Regression</a> <a href="04_binary_and_ordinal_models.html">04<br>GLMs</a> <a href="05_hierarchical_models.html">05<br>Multilevel</a> <a class="active" href="06_model_checking.html">06<br>Checking</a> <a href="07_reporting.html">07<br>Reporting</a> </nav> ``` ## Learning Outcomes By the end of this module, participants should be able to: 1. Check basic MCMC diagnostics. 2. Use posterior predictive checks. 3. Compare models using LOO. 4. Separate computational problems from substantive model problems. ## Two Kinds of Checking Bayesian model checking is doing two jobs at once, and conflating them is the most common source of confusion when participants encounter their first warning message from Stan. The first job is computational: did the sampler explore the posterior properly, or did it fail to mix, get stuck, or signal a problem with the model parameterisation? The second job is substantive: even if the computation succeeded, does the fitted model describe the actual data well enough to support the research claim? A model can fail either check independently. A computational failure produces unreliable posterior summaries no matter how good the model is conceptually; the summaries describe something other than the intended posterior. A substantive failure means the model is fitted correctly to a specification that does not match the data, so the posterior is reliable as a description of an inappropriate model. Both failures need to be diagnosed. Reporting one set of checks without the other leaves the reader unable to judge whether to trust the result. ## Computational Diagnostics The standard computational diagnostics for Hamiltonian Monte Carlo, the algorithm Stan uses, are R-hat, effective sample size, divergent transitions, and trace plots. Each one is checking a different aspect of whether the sampler did its job. R-hat is the ratio of between-chain to within-chain variance for each parameter. If the chains are exploring the same distribution, that ratio should be very close to 1. An R-hat of 1.00 across all parameters is the standard target. An R-hat above 1.01 is a sign that the chains have not converged to the same place; one or more chains is exploring a different region of the posterior. The fix is usually to run longer chains, but persistent R-hat problems indicate that the model is poorly identified and the posterior has multiple modes the sampler cannot reconcile. Effective sample size (ESS) is the number of effectively independent samples in the chain after accounting for autocorrelation. A nominal 4,000 posterior draws can have an ESS of only a few hundred if the chains move slowly through the posterior. Low ESS means posterior summaries, especially extreme quantiles like the 2.5% and 97.5% bounds of a credible interval, are noisy. The standard recommendation is at least 400 effective samples for any quantity the researcher will report. Divergent transitions are a Stan-specific warning that indicates the sampler encountered a region of the posterior with very high curvature that it could not navigate. Even one or two divergences should not be ignored, especially in hierarchical models, where they often signal that the group-level standard deviation is poorly estimated. The usual fix is to increase the `adapt_delta` argument in `brm()` to a value like 0.99, or to reparameterise the model using a non-centred parameterisation. Trace plots show the chain values across iterations, typically with four chains overlaid. A well-mixed model produces traces that look like a "hairy caterpillar": the chains overlap, move around the same central region, and show no obvious trend or stickiness. Traces that diverge, drift, or get stuck in different regions are visible signs of the problems that R-hat and ESS quantify. ```{r} #| eval: false summary(bayes_model) plot(bayes_model) ``` The `summary()` call returns R-hat and ESS for every parameter alongside the posterior summaries. The `plot()` call returns trace plots and density plots for each parameter, four chains overlaid. Reading these two outputs together is the routine first step after any `brms` fit. ## Posterior Predictive Checks Computational diagnostics tell the researcher whether the sampler explored the intended posterior. They cannot tell the researcher whether the intended posterior is a posterior worth interpreting. The substantive check for that is the posterior predictive check, which asks: > If the model were true, would simulated data look like the observed data? The mechanic is straightforward. The fitted model is used to generate many simulated datasets, each one a plausible outcome under the model and the posterior. Those simulated datasets are then compared with the actual data on dimensions the researcher cares about: the shape of the distribution, the mean, the standard deviation, the proportion above a threshold, the variance by subgroup. Systematic disagreement between the simulations and the data is a sign the model is missing something important. ```{r} #| eval: false library(bayesplot) pp_check(bayes_model, ndraws = 50) ``` The default `pp_check()` overlays the density of the observed outcome on the densities of 50 simulated outcomes. If the observed density sits naturally inside the cloud of simulated ones, the model is generating data that looks like the data the researcher actually has. If the observed density is in the tails of the simulated cloud, or has a feature (a second mode, a different skew, a hard upper bound) that the simulations do not reproduce, the model is misspecified in a way the analyst now has to investigate. Targeted checks examine specific summary statistics: ```{r} #| eval: false pp_check(bayes_model, type = "dens_overlay", ndraws = 50) pp_check(bayes_model, type = "stat", stat = "mean") pp_check(bayes_model, type = "stat", stat = "sd") ``` The `type = "stat"` variant computes a summary statistic on each simulated dataset and on the observed data, then plots the distribution of simulated values with the observed value marked. A model that produces a posterior predictive distribution of means that does not contain the observed mean is a model that does not match the central tendency of the data. The same check on the standard deviation catches models that under- or over-estimate the spread, which is a common failure mode for models with the wrong likelihood family. ## What Poor Fit Can Mean A failed posterior predictive check is a diagnosis, not a verdict. It points to a mismatch between the model and the data without specifying which part of the model is wrong. The candidate explanations are usually one or more of the following: the wrong likelihood family for the outcome, a missing predictor whose effect is producing the unexplained structure, a non-linear relationship that the model is forcing into a linear shape, an unmodeled grouping that creates clusters the model does not see, outliers or heavy tails that the chosen likelihood cannot accommodate, priors that are constraining the parameters away from values the data supports, or measurement problems in the outcome itself. Diagnosing which explanation applies requires looking at the specific way the model fails. If the simulated standard deviations cluster well below the observed standard deviation, the model is under-dispersed; a Negative Binomial in place of a Poisson is a common fix for count data with this symptom. If the simulated distributions miss a second mode in the observed data, the model is missing a categorical predictor that separates two groups. If individual residuals are systematically more extreme than the model expects, a likelihood with heavier tails (Student-t in place of Normal) may help. The point is not to chase a perfect model. The point is to know whether the model is adequate for the substantive claim the research wants to make. A model that gets the central tendency right but the tails wrong may be fine for reporting a mean effect, but poor for predicting unusual cases. The check tells the researcher which kinds of claim the model supports. ## Model Comparison with LOO Posterior predictive checks evaluate a single model against the observed data. Model comparison asks a related but different question: among several candidate models, which one is expected to predict new data best? Leave-one-out cross-validation (LOO) approximates this by estimating how well the model would predict each observation if that observation had not been used to fit the model. Computationally, this is done efficiently using importance sampling on the posterior, without refitting the model for each held-out point. ```{r} #| eval: false library(loo) loo_simple <- loo(bayes_model) loo_hier <- loo(hier_model) loo_compare(loo_simple, loo_hier) ``` The output of `loo_compare()` reports the difference in expected log predictive density between models, along with the standard error of that difference. The convention is to interpret a difference larger than 2 standard errors as substantively meaningful evidence that one model predicts better than the other. A smaller difference is consistent with the two models predicting comparably well, in which case the researcher can choose between them on grounds other than predictive performance (parsimony, interpretability, alignment with theory). LOO answers a predictive question, not a universal-truth question. A model with better LOO performance is not automatically more interpretable or more theoretically defensible. A hierarchical model can predict only slightly better than a non-hierarchical one and still be the right choice because it correctly represents the dependence structure. The reverse can also hold: a richer model can predict better and still be the wrong choice if the additional parameters are not interpretable as the substantive quantities the research is meant to estimate. ## Exercise For a fitted model: 1. Check R-hat and effective sample size. 2. Inspect trace plots. 3. Run a posterior predictive density overlay. 4. Compare two models with LOO. 5. Write a short model-checking note. ## Model-Checking Checklist | Check | Question | Result | |---|---|---| | R-hat | Are values close to 1.00? | | | ESS | Are effective sample sizes adequate? | | | Divergences | Are there divergent transitions? | | | Trace plots | Do chains mix well? | | | Posterior predictive check | Do simulated data resemble observed data? | | | Sensitivity | Do conclusions depend heavily on a questionable prior? | | | Substantive adequacy | Is the model good enough for the claim? | | ## Takeaway Bayesian analysis is not complete when the model finishes running. The output of a successful `brm()` call is the beginning of a checking process, not the end of an analysis. A credible workflow runs the computational diagnostics first to confirm the sampler explored the intended posterior, runs the posterior predictive checks to confirm the model describes the observed data adequately for the claim being made, and runs sensitivity checks on the prior choices that materially influence the conclusion. Only then are the posterior summaries safe to report. The reporting notebook covers what to actually write down. ```{=html} <nav class="module-nav" aria-label="Module navigation"> <a href="05_hierarchical_models.html">← Previous: Multilevel Models</a> <a href="07_reporting.html">Next: Reporting →</a> </nav> ```