A Posterior Predictive Check Cannot See What You Adjusted For

Matthew Reda· August 2, 2026· Bayesian workflow

The posterior predictive suite that ships in this framework opens with the mean. It draws replicated KPI series from the fitted posterior, compares their means to the observed mean, and returns a Bayesian p-value. Across 400 simulated weekly panels and five model specifications, one of which gets its search contribution wrong by a factor of three, that check rejected zero times. Not rarely. Zero. The variance check did the same. Both results follow from arithmetic. A linear model with an intercept reproduces the sample mean of its own training data, and a model with a free residual scale reproduces the sample variance. A check aimed at a quantity the fit reproduces has already passed before you run it. The principle generalizes past the mean. Residuals are orthogonal to every column of the design, so a check built from the correlation between residuals and a variable you adjusted for is dead on arrival too. What survives is a thin band of test quantities that can see the shape of the fitted mean over time. They see the shape. They never see the split.

What a Model Absorbs, It Cannot Fail

The default check list lives in PPCConfig (validation/config.py:26-38) and reads ("mean", "variance", "autocorrelation", "skewness", "extremes"), run at a significance level of 0.05 over 500 replicated datasets. Each entry resolves to a class in validation/posterior_predictive.py:37-263, and every class does the same thing: compute a statistic on the observed series and on each replicated series, then count how often the replicated value lands at or above the observed one. The comment above the first p-value in the file states the convention plainly.

Fit a linear model with an intercept, and the fitted values average to the observed average by construction. The replicated series differ from that fitted mean only by simulated noise, centered on zero with a standard error of \( \sigma/\sqrt{n} \). So the observed mean sits in the dead center of the replicated distribution, every time, on any dataset, under any specification carrying a constant. The p-value concentrates at 0.5. Its spread is far tighter than the check's own rejection thresholds can reach.

The same holds for the variance, since a model with a free residual scale sets that scale to whatever variation is left over. Two of the five checks in the default suite are structurally incapable of returning a failure, and they are the two an analyst reads first.

Absorption is the general rule, and it extends well past the intercept. Write the fitted design as \( X \) and its hat matrix as \( H = X(X'X)^{-1}X' \). The least-squares residual is \( (I - H)y \), which is orthogonal to every column of \( X \) by definition. Any test quantity that reduces to a linear function of the residuals inside the span of \( X \) evaluates to approximately zero on the observed data, no matter how badly the model is misspecified, while the replicated version fluctuates around zero at random. The observed value sits at the center of the replicated distribution and the check passes. That covers the correlation between residuals and any control you included. It comes close on the average residual in weeks a channel went dark, where the fitted columns explain about half of the dark-week indicator and leave the check with almost nothing to read.

The Quantity Under Test

Gelman, Meng and Stern (1996) drew the distinction that makes the rest of this tractable. A test statistic \( T(y) \) is a function of the data alone. A discrepancy \( T(y, \theta) \) may use the parameters as well, which lets it measure something like a standardized residual that has no meaning until you fix a value of \( \theta \). The Bayesian version compares the realized discrepancy to the replicated one draw by draw, holding \( \theta \) fixed within each comparison.

\( p_B \) is not a frequentist p-value. The data enter twice, once to form the posterior and again to form the observed side of the comparison, so its null distribution is not uniform and its concentration around 0.5 varies by test quantity. That concentration is the property you have to measure before you can read a single value. A check whose p-value only ever occupies the interval from 0.45 to 0.55 will never fire the rule in the file. Modrák and colleagues (2023) make the equivalent point for simulation-based calibration: the choice of test quantity determines sensitivity, and a suite of them can be blind to a fault that one well-chosen quantity exposes.

So I measured it. Below is the null calibration of ten test quantities, computed on 400 panels drawn from a planted world and fitted with the correct specification. The first five mirror the checks that ship. Realized discrepancies built from standardized residuals make up the other five, and they include the three that the earlier post on holdout error named as the ones worth building.

Five of the ten cannot reject a correctly specified model, which is the good news. Four are the absorbable quantities the orthogonality argument predicts: the mean, the variance, the residual correlation with search adstock, which is a column of the design, and the demand index, in the fitted design here and absorbed for the same reason. The fifth is the mean residual in dark weeks, which the design comes close to spanning without containing. Five survivors run between 3.2% and 7.0% against a nominal 10%, so every one of them is conservative.

Ten Checks, Four Faults

Calibration tells you the false-alarm rate. Severity is the other half, and Gelman and colleagues are blunt about how little general guidance exists for it. The workflow paper behind their 2026 book says outright that “there is no general way to choose which checks one should perform on a model,” and points at Mayo’s severe-test criterion for the shape of an answer: what you want are “checks that are likely to fail if the model would give misleading answers to the questions we care most about” (§6.1). Likely to fail is a probability. Estimate it against a given fault and you have a number to put beside a check instead of a hope.

The world is a three-year weekly panel. A brand channel runs in flights and is dark in about 62 of the 156 weeks. A search channel has its budget set against a category demand index, which also raises sales directly, at a correlation of about 0.95 between spend and demand. Media contributes 9.1% of the KPI in truth. Everything is linear in fixed transforms, so the estimating model is a Bayesian linear regression under a reference prior with an exact posterior, no sampler and no prior doing hidden work. Four faults, each one a competent analyst could ship: drop the demand index, assume no carryover on brand when the truth carries at 0.5, add a placebo channel correlated with brand at 0.8 and no true effect, or replace the two-harmonic seasonal basis with eight harmonics and a cubic trend.

Median over 400 panels. Two of the four faults are expensive. Dropping the demand index costs 2.7 points of in-sample R². It triples the search contribution and doubles working media, from 9.1% of the KPI to 18.6%. Assuming away carryover halves the brand estimate. The other two barely move the answer in this world. That matters at the next table, where an invisible fault that costs nothing is not a scandal.

Read the rows first. Four rows sit at 0.0% across every column, including the column where the model triples a channel, and the residual against the demand index joins them everywhere the demand index is in the design. Nothing about the size of the fault rescues them, because the fit reproduces the quantity being tested and the fault has no way to enter. Widening the rejection thresholds would not move them either. The best-performing check in the whole set is the lag-1 autocorrelation of the KPI, at 33.5% and 37.2% against the two expensive faults, which means it misses each of them about two panels in three. Its residual counterpart, the statistic I would have picked on theory, does worse on both.

Then read the bottom row. Ten checks at a nominal 10% each flag a correctly specified model on 21.8% of panels, well under the 65.1% that ten independent uniform p-values would produce, because the conservatism runs deep and the checks are dependent. Against the confounded model the family rate is 52.5%. A team running the full suite, in a world where their media number is triple the truth, sees nothing unusual on about half of all panels, and an autocorrelation flag on four panels in ten. That flag says the temporal structure is off. It says nothing about the search coefficient absorbing a demand index.

That is the honest summary of what these checks buy. They have real power against faults that misstate the shape of the fitted mean over time, which is worth having, and the previous post on holdout error ranked them fourth of five checks by the failure space each one covers. They have no power at all against a fault that lives entirely inside the column space.

Every Control You Add Is a Check You Give Up

The one row with any signal against confounding is the residual correlated against the demand index, and it manages 4.2%. Figure 14D of the workflow paper illustrates the same instinct with a predictive check grouped by a covariate left out of the model, and the reason this one barely works here is worth the algebra.

Turn the demand-chasing dial on the search budget and both quantities move, in ways that have nothing to do with each other. The parameter \( \kappa \) below is how strongly weekly search spend tracks the demand index. Raising it makes the confounder easier for the fitted columns to explain, which is what \( R^2 \) measures here, and the check goes blind somewhere between 0.80 and 0.93.

Median over 250 panels per row. The first row is the honest case: no confounding, an error of 6.3% that is mostly the noise floor of a weakly identified coefficient, and a check that fires on every single panel because the omitted variable is real and independent. It is working perfectly and reporting a fault that costs nothing. Three rows down the estimate is wrong by a factor of three and the check fires 5.2% of the time, which is its false-alarm rate. It has stopped being a test.

There is a practical rule buried in that table, and it is uncomfortable. Adjusting for a variable and testing against it are competing uses of the same column. Put the demand index in the design and the residual is orthogonal to it, so no residual check can ask whether it belonged there. Leave it out while the other columns explain little of it, and the check fires on every panel. The fault it flags costs 6%. Leave it out while a channel explains 93% of it, and the estimate moves by a factor of three with the check silent. Same variable in all three, and only the modeling decision changed.

A check suite could hand that trade back to the analyst. The information a posterior predictive check needs has to come from outside the fitted columns, so name the outside explicitly.

# illustrative
# Declare, before the fit, covariates you will NOT adjust for, and check
# the residuals against them. The test is bought with an adjustment, so
# the declaration belongs in the pre-registered spec.
check = HeldOutCovariateCheck(
    columns=["category_demand_index", "competitor_price"],
    statistic="residual_correlation",
)
result = check.compute(y_obs=y, y_rep=y_rep, residuals=r, frame=panel_extras)
# result.p_value only reads next to result.power_against: the rejection
# rate this statistic reaches on a planted fault of the size you care about.

That second comment is where the field stops reading. A p-value from a check whose power against the fault of interest is 5% carries no information, and that power is computable whenever you will plant the fault in a simulation, which the framework's synthetic worlds already do for other purposes.

Leaving a Week Out Does Not Help

One repair for a p-value that uses the data twice is to stop using it twice. The leave-one-out predictive interval transform, LOO-PIT, evaluates each observation against a posterior predictive distribution fitted without that observation, and under a well-specified model the resulting values are uniform on the unit interval by construction. Vehtari, Gelman and Gabry (2017) give the importance-sampling machinery that makes this affordable for a single fit, along with the Pareto-k diagnostic that says when the approximation itself is straining.

For the Gaussian linear model above, LOO-PIT is available in closed form from the externally studentized residuals, so I computed it exactly on the same 400 panels and tested each panel's PIT values against uniformity with a Kolmogorov-Smirnov statistic. Calibrating the critical value on the correct specification puts it at KS = 0.0706 for a 5% false-alarm rate. Against the four faults the rejection rates are 4.8% for the omitted demand index, 6.5% for the missing carryover, 4.8% for the placebo channel and 5.5% for the flexible baseline. The double use of the data is gone and the blindness is untouched, because LOO-PIT is still a statement about how well the model predicts the KPI. It is a sharper instrument aimed at the same quantity.

One caveat on the method. Leave-one-out assumes exchangeable observations, and weekly MMM data with adstocked regressors is not exchangeable, so dropping a single week leaves most of its information sitting in its neighbours. Bürkner, Gabry and Vehtari (2020) give the leave-future-out construction that respects the time ordering, and it is the right tool on real panels. Ordinary LOO is exact in my planted world, whose errors are independent by construction. The blindness above does not come from the cross-validation scheme.

What This Does Not Establish

Absorption is exact in a linear-Gaussian model fitted under a flat prior, which is what I simulated. A real MMM has informative priors, and a coefficient held down by a tight prior cannot absorb as freely, so a strongly regularized model leaves more of a fault in its residuals than the tables above suggest. Saturation curves have a similar effect, since the residual is then only approximately orthogonal to the transformed columns. Both push toward more detectability, and I have not measured how far.

The rejection rates are power against four faults I chose, on one data-generating process, at one signal-to-noise ratio. A different catalogue reorders the table, and the two faults that came out invisible here were also nearly harmless here, so this says nothing about whether placebo channels and over-flexible baselines are generally safe. The holdout-error post planted a world where an over-flexible baseline moved the estimand by 20%.

Nothing here says posterior predictive checking is worthless. It caught the carryover misspecification more often than anything else in the suite, and a model that reproduces the KPI's mean and variance while failing on autocorrelation is telling you something true about its temporal structure. The claim is narrower, and it is about coverage of the failure space: this family cannot audit the decomposition, and a suite whose output is a row of green ticks invites the reader to believe it did.

Where This Goes Next

Two repairs are buildable now. One is the held-out covariate check sketched above, with the held-out set pinned in the pre-registered spec so that the choice is made before anyone has seen a residual. The other is a power column beside every check in the validation report. Plant a fault of a declared size in a synthetic world matched to the panel, count rejections, print the rate next to the p-value. That turns each check from a claim into a measurement.

A third repair stays out of reach. A discrepancy computed on the decomposition rather than on the KPI would need a replicated quantity to compare the split against, and the only sources of one are a planted truth or an experiment. Both sit outside the posterior predictive machinery.

The prior predictive and calibration post closes by handing the question of whether the model resembles reality to posterior predictive checks on real data. No simulation-based check can do that job, so the handoff is right. The receiving check can audit the temporal shape, and nothing about the split.

References

• Gelman, A., Meng, X.-L., & Stern, H. (1996). Posterior Predictive Assessment of Model Fitness via Realized Discrepancies. Statistica Sinica, 6(4), 733–807.
• Gelman, A., Vehtari, A., Simpson, D., Margossian, C. C., Carpenter, B., Yao, Y., Kennedy, L., Gabry, J., Bürkner, P.-C., & Modrák, M. (2020). Bayesian Workflow. arXiv:2011.01808. Section 6.1 on posterior predictive checking and severe tests.
• Gelman, A., Vehtari, A., McElreath, R., and colleagues (2026). Bayesian Workflow. CRC Press, 538 pp. ISBN 978-0-367-49014-0.
• Mayo, D. G. (2018). Statistical Inference as Severe Testing: How to Get Beyond the Statistics Wars. Cambridge University Press.
• Gabry, J., Simpson, D., Vehtari, A., Betancourt, M., & Gelman, A. (2019). Visualization in Bayesian Workflow. Journal of the Royal Statistical Society, Series A, 182(2), 389–402.
• Vehtari, A., Gelman, A., & Gabry, J. (2017). Practical Bayesian Model Evaluation Using Leave-One-Out Cross-Validation and WAIC. Statistics and Computing, 27(5), 1413–1432.
• Bürkner, P.-C., Gabry, J., & Vehtari, A. (2020). Approximate Leave-Future-Out Cross-Validation for Bayesian Time Series Models. Journal of Statistical Computation and Simulation, 90(14), 2499–2523.
• Modrák, M., Moon, A. H., Kim, S., Bürkner, P.-C., Huurre, N., Faltejsková, K., Gelman, A., & Vehtari, A. (2023). Simulation-Based Calibration Checking for Bayesian Computation: The Choice of Test Quantities Shapes Sensitivity. Bayesian Analysis.