Posterior Predictive Checks and the Bayesian Workflow

A Poisson regression for daily support tickets per customer converges cleanly: \(\hat{R} = 1.00\), no divergences, tight intervals. Its fitted mean is about 2 tickets. A Poisson distribution with mean 2 puts \(e^{-2} \approx 13.5\%\) of its mass at zero, and 40% of real customer-days have no tickets. The sampler worked perfectly on a model that cannot represent the data, and nothing in the convergence diagnostics could have said so. Only comparing data the model generates with the data that exist exposes it.

The check, precisely

The posterior predictive distribution generates replicated datasets from the fitted model:

In practice, for each posterior draw \(\theta^s\), simulate one dataset \(y^{\text{rep},s}\) of the same size and structure as \(y\). Then choose a test quantity \(T\), which may depend on parameters as well as data, and compare \(T(y, \theta^s)\) with \(T(y^{\text{rep},s}, \theta^s)\) across draws. The posterior predictive p-value is

estimated as the fraction of draws where the replicate is at least as extreme (Gelman, Meng and Stern, 1996, Posterior Predictive Assessment of Model Fitness via Realized Discrepancies, Statistica Sinica 6). In the ticket example, \(T\) = proportion of zeros gives a replicated range around 13% against an observed 40%, and \(p_B \approx 0\).

Choosing what to test

The test quantity is the entire design of the check. Statistics that the model fits directly are useless: a normal model with a free mean reproduces the sample mean every time, so \(T = \bar{y}\) always passes. Useful statistics target features the model does not explicitly fit: proportion of zeros, maximum, skewness, the variance-to-mean ratio for counts, lag-1 autocorrelation of residuals, the spread of group means in a hierarchical model. Each test should correspond to a way the model could be wrong that would change a decision.

Graphics usually beat p-values. Overlaying the density of the observed data on densities of dozens of replicates, or plotting the distribution of \(T(y^{\text{rep}})\) with the observed value marked, shows how the model fails, not just that it does (Gabry, Simpson, Vehtari, Betancourt and Gelman, 2019, Visualization in Bayesian workflow, JRSS-A 182(2)).

The rest of the workflow

Posterior predictive checks are one stage of an iterative loop (Gelman et al., 2020, Bayesian Workflow, arXiv:2011.01808).

Prior predictive checks simulate data from the prior alone, before touching \(y\). If a logistic regression's priors imply that most predicted probabilities sit at 0 or 1, the prior is making a strong claim nobody intended.

Fake-data simulation fixes parameters, simulates data, fits, and checks recovery. If the model cannot recover known parameters from data it generated itself, fitting real data is pointless.

Simulation-based calibration makes that systematic. Draw \(\tilde\theta\) from the prior, simulate \(\tilde y\), draw \(L\) posterior samples given \(\tilde y\), and record the rank of \(\tilde\theta\) among them. If the computation is correct, the rank is uniform on \(\{0, \dots, L\}\), and characteristic departures (a U shape for under-dispersed posteriors, a skew for bias) point at the failure (Talts, Betancourt, Simpson, Vehtari and Gelman, 2018, Validating Bayesian Inference Algorithms with Simulation-Based Calibration, arXiv:1804.06788).

Computational diagnostics such as \(\hat R\), effective sample size and divergences (covered in the MCMC concept on this track) establish that the posterior was computed, not that the model is adequate.

Model expansion and comparison follows failure: a zero-inflated or negative binomial model for the tickets, then PSIS-LOO to compare predictive performance.

When it breaks

Posterior predictive p-values are not calibrated. Because \(y\) is used both to fit and to check, replicates are pulled toward the observed data, and under a correct model \(p_B\) concentrates near 0.5 rather than being uniform. The check is conservative: an extreme value is strong evidence of misfit, but a moderate value is weak evidence of fit (Bayarri and Berger, 2000, P Values for Composite Null Models, JASA 95(452)). Bayarri and Berger proposed calibrated alternatives; Gelman and colleagues argue that the check is an exploratory tool for finding misfit, not a hypothesis test, so calibration matters less. That disagreement has not been settled.

Flexible models pass everything. A model with enough parameters reproduces any statistic of the training data, including its noise. Passing posterior predictive checks says nothing about out-of-sample prediction; that is a cross-validation question.

Iterating on the data is a garden of forking paths. Each round of check, expand, refit makes the final model a function of the data in ways its posterior does not account for, so the reported uncertainty is too narrow. The workflow's advocates accept the cost and recommend documenting the sequence of models; critics see it as the Bayesian version of specification search.

Simulation-based calibration is expensive. A useful SBC run needs hundreds of simulated datasets, each a full fit. For a model that takes ten minutes to sample, 500 replications is over three days of compute, which is why it is typically run on a simplified model or a subset of parameters.

Hierarchical replicates need a choice. Replicating new observations within existing groups and replicating entirely new groups test different things, and a model can pass one while failing the other.