Posterior Predictive Checks For Regression Models

Why Check a Fitted Regression?

A regression model can produce plausible coefficient summaries and still fail to describe important features of the observed outcomes. Its mean structure may miss curvature, its errors may have heavier tails than assumed, or its variance may differ across groups. Posterior predictive checks ask a direct question: if this fitted model were a useful description of the data, would data generated from it look like the data we observed?

The check is a model diagnostic, not a test that proves a model true or false. A mismatch can identify a feature the model struggles to reproduce, but it does not by itself reveal the unique reason for that mismatch. The choice of diagnostic matters: a check can only be informative about features it actually measures.

Model and Posterior Predictive Distribution

For a basic regression check, suppose the observed outcomes $y_1,\ldots,y_n$ are modeled conditional on a fixed design matrix $X$. In a Gaussian linear regression, the model is

y_i = \mathbf{x}_i^\top\boldsymbol{\beta}+\varepsilon_i,\qquad \varepsilon_i\mid\sigma^2 \overset{\mathrm{ind}}{\sim} N(0,\sigma^2),\qquad i=1,\ldots,n.

Here $\mathbf{x}_i$ is the regressor vector for observation $i$, $\boldsymbol{\beta}$ is the coefficient vector, and the errors are conditionally independent, have mean zero, and share the same variance. The mean is linear in the included regressors, and the Gaussian likelihood specifies the conditional outcome distribution. These are modeling assumptions to examine, not facts guaranteed by fitting the regression.

Let $\theta=(\boldsymbol{\beta},\sigma^2)$. Given a prior and the observed data, the posterior is $p(\theta\mid y,X)$. A replicated outcome vector $y^{\mathrm{rep}}$ is generated by first drawing $\theta$ from that posterior, then drawing outcomes from the likelihood at the observed regressor values.

Definition: The posterior predictive distribution conditional on the observed design is $$ p(y^{\mathrm{rep}}\mid y,X) = \int p(y^{\mathrm{rep}}\mid\theta,X)\,p(\theta\mid y,X)\,d\theta. $$ It averages the likelihood for new outcomes over posterior uncertainty about the regression parameters. In a check conditional on $X$, the regressors are held fixed; they are not regenerated.

For the Gaussian regression, one replicate uses

y_i^{\mathrm{rep},(m)} = \mathbf{x}_i^\top\boldsymbol{\beta}^{(m)} +\varepsilon_i^{\mathrm{rep},(m)},\qquad \varepsilon_i^{\mathrm{rep},(m)}\sim N(0,\sigma^{2(m)}),

where $\theta^{(m)}=(\boldsymbol{\beta}^{(m)},\sigma^{2(m)})$ is a posterior draw. Repeating this for many draws approximates the posterior predictive distribution. The posterior draws should come from chains with adequate convergence and mixing; otherwise apparent diagnostic patterns may reflect Monte Carlo error rather than the posterior.

Discrepancy Measures and Their Interpretation

A discrepancy measure $D(y,\theta)$ summarizes a feature of the data that matters for the model assessment. It may depend on parameters so that, for each posterior draw, the observed data and its replicate are evaluated on a comparable scale. Examples include the absolute slope of residuals against a regressor, the largest absolute standardized residual, or a ratio of group-specific residual variances.

Definition: A posterior predictive probability for a discrepancy is $$ p_{\mathrm{ppc}} = P\!\left(D(y^{\mathrm{rep}},\theta)\geq D(y,\theta)\mid y,X\right). $$ With $M$ posterior draws and one replicate per draw, estimate it by $$ \widehat p_{\mathrm{ppc}} = \frac{1}{M}\sum_{m=1}^{M} I\!\left\{D(y^{\mathrm{rep},(m)},\theta^{(m)}) \geq D(y,\theta^{(m)})\right\}. $$

This probability is a measure of how often replicated data show at least as large a discrepancy as the observed data, under the fitted posterior predictive distribution. It is not the probability that the model is true, and it is not a classical p-value with a universal null distribution or a guaranteed uniform calibration. Since the observed data were used to form the posterior and also enter the comparison, posterior predictive probabilities are often concentrated away from extreme values. Interpret them alongside plots and substantive knowledge, not using a rigid cutoff as a decision rule.

Checking Residual Patterns

For each posterior draw, calculate observed residuals $r_i^{(m)}=y_i-\mathbf{x}_i^\top\boldsymbol{\beta}^{(m)}$ and replicated residuals $r_i^{\mathrm{rep},(m)}=y_i^{\mathrm{rep},(m)}-\mathbf{x}_i^\top\boldsymbol{\beta}^{(m)}$. Plot each set against fitted values, relevant regressors, or observation order. Under the stated linear mean structure, a systematic trend in observed residuals that is rarely reproduced can indicate that the mean function is inadequate. A changing residual spread can instead suggest that constant variance is inadequate.

Worked Example: A Residual Trend Check

Suppose a hypothetical model relates an outcome to a regressor $x$, and the observed values are $x=(0,1,2,3,4)$ and $y=(2,4,7,9,12)$. For illustration, consider one posterior draw with intercept $2$, slope $2$, and error standard deviation $1.5$. The fitted values are $(2,4,6,8,10)$, giving observed residuals $(0,0,1,1,2)$. The data are being checked conditional on these five observed regressor values, under the linear, independent, constant-variance Gaussian model.

A simple discrepancy is the absolute least-squares slope of residuals against $x$. Here $\bar{x}=2$, the residual mean is $0.8$, and

\widehat b_r = \frac{\sum_i(x_i-\bar{x})(r_i-\bar{r})}{\sum_i(x_i-\bar{x})^2} = \frac{5}{10} = 0.5.

Thus the observed discrepancy for this draw is $|\widehat b_r|=0.5$ outcome units per unit of $x$. Replicates are generated at the five listed $x$ values from the draw's fitted means and Gaussian error distribution. Across 100 paired posterior-draw comparisons, suppose the replicated absolute residual slope was at least the corresponding observed slope in 18 comparisons. The estimated posterior predictive probability is $18/100=0.18$. This is not by itself decisive evidence for or against the model; the residual plot and the pattern across draws should also be examined.

The residual check should target a plausible failure. A smooth curve in residuals against a regressor points toward a mean-structure issue; a widening band may point toward nonconstant variance. A plot against observation order can reveal structure the model omits, although independence is a substantive assumption and should not be inferred from a plot alone.

Checking Tails and Group Variation

The Gaussian likelihood implies more than a conditional mean: it also implies a particular tail behavior and, in this model, common error variance. Checks based on extreme residuals or group-specific spread can reveal a poor match in these features. For each posterior draw, standardized residuals can be formed as $z_i^{(m)}=r_i^{(m)}/\sigma^{(m)}$. For a tail check, a discrepancy might be $D(y,\theta)=\max_i|z_i|$, or the number of standardized residuals exceeding a specified magnitude. The same calculation is applied to the replicate.

Worked Example: Checking an Extreme Residual

Consider a hypothetical intercept-only Gaussian regression for eight outcomes. In one posterior draw, let the fitted mean be $10$ outcome units for each observation and let $\sigma=1.2$ units. Suppose the observed outcomes are $(10.2,9.6,10.1,10.5,9.7,10.7,9.4,13.0)$. The model assumes conditionally independent Gaussian errors with this shared standard deviation; the example illustrates a tail check under that assumption.

The residuals are $(0.2,-0.4,0.1,0.5,-0.3,0.7,-0.6,3.0)$, so the largest absolute standardized residual is

D(y,\theta)=\max_i\left|\frac{y_i-10}{1.2}\right| =\frac{3.0}{1.2}=2.5.

Generate replicated sets of eight outcomes from the draw's fitted model and calculate the same maximum for each set. Suppose that across 100 paired posterior comparisons, 6 replicated maxima were at least as large as the observed maximum for their corresponding draws. Then $\widehat p_{\mathrm{ppc}}=6/100=0.06$. In this hypothetical check, the observed extreme is relatively unusual among the replicates, warranting closer inspection of the tail fit; it does not establish that the Gaussian model is false.

For group variation, first decide whether the model includes group-specific means or other group structure. A useful check of a common-variance assumption compares residual spread across groups after accounting for the modeled means. A discrepancy might be the largest group residual variance divided by the smallest. If a group has too few observations for a stable variance estimate, the statistic itself may be noisy; show the group-level data and interpret the check cautiously.

Worked Example: Comparing Residual Variances Across Groups

Suppose a hypothetical regression includes group-specific means but assumes one common Gaussian error variance. After subtracting the fitted mean for each group, the observed residuals are $(-1,1)$ in group A and $(-3,3)$ in group B, in outcome units. Each group has two observations. The model's common-variance and conditional-independence assumptions are the features being checked; the group means themselves have already been allowed to differ.

Using the sample variance within each group, group A has $s_A^2=\{(-1)^2+1^2\}/(2-1)=2$, while group B has $s_B^2=\{(-3)^2+3^2\}/(2-1)=18$. The observed variance-ratio discrepancy is therefore $18/2=9$. Suppose that among 100 paired posterior predictive comparisons, 4 replicated ratios were at least as large as the observed ratio. The estimated posterior predictive probability is $4/100=0.04$. This suggests the observed difference in group spread is rarely reproduced by this fitted common-variance model, but the very small groups also make the sample variances unstable.

Practical Workflow and Common Pitfall

State the feature of interest.
Choose a residual pattern, tail feature, group comparison, or other discrepancy tied to a specific model assumption.

Generate conditional replicates.
For each posterior draw, simulate outcomes using the fitted likelihood and hold observed regressors fixed unless the check explicitly models their distribution.

Compare on a like-for-like basis.
Calculate the same discrepancy for the observed data and the paired replicate, using the same draw-specific parameters where the statistic depends on parameters.

Inspect plots and uncertainty.
Look for systematic differences across many replicated data sets, and consider Monte Carlo uncertainty and the stability of the chosen statistic.

A common mistake is to say, “The posterior predictive probability is $0.04$, so there is a 4% probability the regression model is correct.” That reverses the meaning of the calculation. The value estimates how often the specified replicated discrepancy reaches or exceeds the observed discrepancy, conditional on the fitted posterior and the chosen statistic. It does not assign a probability to model truth. Another mistake is to check only a single overall fit summary: a model may reproduce the outcome mean while missing curvature, tail behavior, or group variation. A small set of targeted checks is more informative than an unexplained collection of numbers.

Posterior predictive checks matter because regression conclusions depend on the likelihood and mean structure as well as the posterior coefficient summaries. If a check reveals misfit, report which feature is poorly reproduced and consider whether the model needs a different mean function, variance structure, dependence structure, or outcome distribution. The check is diagnostic evidence, not a substitute for a well-justified model or for uncertainty-aware interpretation.

Check Your Understanding

Answer each question by referring to the model, the discrepancy, and what the comparison can establish.

In a posterior predictive check conditional on $X$, which quantities are held fixed, and which are simulated?
Why is a residual slope against a regressor a relevant discrepancy for checking a linear conditional mean?
A posterior predictive probability is estimated as $0.08$. State what this number estimates, and one thing it does not mean.
If observed group residual variances differ substantially but replicated variances usually do not, which model feature may deserve examination?
Why should the observed and replicated discrepancies be calculated using paired posterior draws when the discrepancy depends on parameters?