Your MMM Is Reporting the Wrong Cross-Validation

Matthew Reda· August 2, 2026· Bayesian workflow

Leave-one-out on a weekly series holds out one week and trains on every week that came after it. That is a prediction nobody in marketing ever makes. On the 156-week simulation below, deleting the two weeks either side of the held-out one reverses which of two models wins, and the forward fold agrees with the reversal.

The Fold You Ran and the Fold You Meant

An MMM validation run in this framework calls az.loo and stores four scalars: elpd_loo, its standard error, p_loo, and a count of observations whose Pareto k̂ exceeds 0.7, plus the pointwise Pareto-k and pointwise elpd arrays (validation/validator.py:367-385). Every one of those is computed correctly. The question is what they are computed about.

Leave-one-out estimates the expected log predictive density for a new observation given every other observation in the dataset, and on a weekly panel indexed by time that phrase quietly includes the weeks after the one you dropped. Hold out week 78 of a three-year series and the model still sees weeks 79 through 156. It knows where the baseline went. A media planner asking what happens next quarter has no such thing, and neither does the backtest that grades a plan.

How much this matters depends on how local the model's time basis is. A parametric baseline of intercept, linear trend and two annual harmonics is estimated from all 156 weeks at once, so removing one of them barely moves it. Contrast a local-level or random-walk baseline, which is what most MMM trend components are: at week 78 it is pinned almost entirely by weeks 77 and 79, and dropping week 78 from the likelihood leaves those two neighbours to rebuild it out of information a forecaster would never hold. The flexible model is being scored on a task its flexibility cannot fail.

I ran that comparison exactly, with no sampler and no importance sampling anywhere. The world is 156 weekly observations: three channels put through geometric adstock and a Hill curve at fixed transform parameters, an AR(1) residual with autocorrelation 0.6, and fifteen promotional weeks whose noise standard deviation is four times the rest. Both candidates are linear in their parameters with Gaussian priors and a fixed observation standard deviation of 0.391, so every posterior and every refit is a closed form. The local-level candidate carries 159 parameters (three media coefficients plus a weekly random walk), against nine for the parametric one. Both are scored on the same 52 weeks, numbers 52 through 103, so the only thing that changes across the rows below is which weeks were available for training.

Δelpd is local-level minus parametric, so a positive number favours the flexible baseline. Read the first row and the last. Row-wise LOO puts the local-level model ahead by 8.8 nats with a standard error of 5.1, the kind of gap that gets written into a deck as "the state-space specification fits materially better." An honest one-step-ahead forward fold puts them 0.5 nats apart with a standard error of 7.9, which is nothing at all.

The second row is the one I would put on a slide. A leave-3-out block deletes the held-out week and its two immediate neighbours, so the training set shrinks from 155 weeks to 153. Two observations out of 155. The verdict goes from +8.8 to −1.7 and changes sign. Nothing else about the model, the data or the scored weeks moved.

One simulated world proves nothing, so I repeated it across 20 seeds. Row-wise LOO favoured the local-level model in 12 of the 20 worlds and leave-future-out favoured it in 4; the two schemes disagreed about the winner in 10 of 20. Between LOO and leave-3-out the direction of the shift was the same every time: LOO scored the local-level model higher relative to the parametric one in 20 of 20 worlds, by a median of 2.2 nats and never by less than 0.4. Broken out per model, the amount that row-wise LOO adds over leave-3-out has median 3.8 nats for the local-level candidate and 1.6 for the parametric one, and the flexible model got the larger boost in all 20.

A shift that keeps its sign across 20 worlds is a bias, and it points at the flexible model.

Time Is One Direction. Geography Is the Other.

Gelman and co-authors are explicit about this in the cross-validation section of Bayesian Workflow. Writing about hierarchical structure, they say that "cross validation for multilevel (hierarchical) models requires more thought. Leave-one-out is still possible, but it does not always match our inferential goals," and point at leave-one-unit-out and leave-one-group-out as the alternatives whose folds match the question. They add that importance sampling for those folds "tends to be much harder than in the leave-one-out case," which is a polite way of saying you will be doing refits.

A geo MMM has both structures at once. Along time, the fold that matches the question is leave-future-out: train to a cutoff, predict past it, roll the cutoff forward (Bürkner, Gabry and Vehtari, 2020). Along geography, the fold that matches the question is leave-one-group-out: hold out a DMA entirely and ask whether the model would have got that market right from the others. Those answer different business questions. One is whether a plan built today survives contact with next quarter. The other is whether the partial pooling across markets is real, which is exactly the claim a per-geo media coefficient makes when vary_media_by_geo is on.

Row-wise LOO on a geo panel answers neither. Its fold drops one geo-week cell and keeps the same week in every other market and the same market in 155 other weeks, so both hierarchies leak into the prediction at once. The elpd it returns is a real quantity about a real predictive task. That task is interpolating a hole in a panel you already have.

What the Scalar Throws Away

The same section of Bayesian Workflow lists three uses for cross-validation, and the headline number is the least of them. One is calibration of the cross-validated predictive distribution, read off the uniformity of the LOO-PIT values. The other two are identifications: which observations the model finds hardest to predict, and which ones are most influential, the second of those read off the importance weights themselves. A report that quotes elpd_loo and stops has paid for all three and kept the thinnest.

Take calibration. In the simulation above, the noise is heteroskedastic by construction and both candidates assume it is constant, so their predictive distributions are the wrong width in a specific pattern. The exact LOO-PIT values over all 156 weeks make that visible. For the local-level model the Kolmogorov-Smirnov statistic against uniform is 0.126 with p = 0.013; its central 50% predictive band holds 71% of the weeks rather than 50%, and its outer 10% band holds 6% rather than 10%. Both readings say the intervals are too wide for a typical week. Then look at which weeks did land outside: all 9 of them are promotional weeks. The intervals are simultaneously too wide for the ordinary weeks and too narrow for the ten percent of weeks that matter to a promotional plan. A single elpd number cannot express that, and neither can a MAPE.

This framework already computes the diagnostic. validation/calibration.py:48 (loo_pit_check) returns the PIT values, a KS uniformity test and a calibrated flag, and reporting/interactive/facts.py:424 feeds it real PSIS log-weights plus simultaneous ECDF-difference bands for the interactive results report. What it computes them on is the row-wise fold. LOO-PIT built on a fold that trains on the future inherits the same complaint as the elpd built on it: the calibration it certifies is the calibration of an interpolation.

The 0.7 Threshold Is a Function of Your Draw Count

Everyone quotes k̂ ≤ 0.7. That number is a large-sample rule. The JMLR treatment of Pareto smoothed importance sampling (Vehtari, Simpson, Gelman, Yao and Gabry, 2024) makes the threshold depend on how many draws you have:

What makes a k̂ above threshold fatal has nothing to do with a known amount of bias. The number of draws needed for a reliable Monte Carlo estimate grows infeasibly fast past that point, which is why the posterior package now ships pareto_min_ss() and pareto_khat_threshold() alongside pareto_khat(). Read honestly, the diagnostic reports a required sample size rather than a pass or a fail.

In the simulation, the local-level model at S = 500 draws has three weeks above the 0.629 threshold for that budget and only one above 0.7. At S = 2000 it has two above 0.697, with a maximum k̂ of 0.90. The parametric model never exceeds 0.52. So the flexible baseline is the one straining the approximation, and at a modest draw budget the folklore threshold clears two of its worst observations.

Five places in this codebase hard-code 0.7: validation/validator.py:378, the influential-observation flag at validation/validator.py:1254, the pareto_k_threshold default at validation/config.py:95, validation/spec_curve.py:441, and reporting/interactive/facts.py:466. For a default NUTS fit of 4 chains by 1000 draws that is correct, and for a backtest refit at 4 chains by 500 (BacktestConfig.draws) it is correct to within 0.003. Below 2000 it is not. Trim the draws for a reconnaissance fit, for a reduced-budget spec-curve sweep, or on a 40-minute-per-fit production model to get it through the night, and the threshold the code checks against has quietly become too permissive. Quoting 0.7 without stating S is quoting a rule you have not checked applies.

Where the Shortcut Breaks, and What Does Not Repair It

PSIS-LOO worked fine in my simulation. At S = 2000 it estimated the local-level model's elpd within 0.53 nats of the exact leave-one-out answer, and the parametric model within 0.10. Two linear models with Gaussian posteriors are an easy case, and I want that on the record before the next paragraph, because the failure mode is structural rather than universal.

Bindoff (2026) studies where it fails and finds the failures are predictable in advance. The Gelman-Pardoe pooling factor for group \(j\),

together with each observation's pull on its own group's random effect, predicts which folds will come back with k̂ > 0.7 at an AUC of 0.96 in Gaussian linear mixed models and 0.81 in replicated logistic GLMMs. In the Gaussian case the per-observation Fisher information is constant, so the map reduces to group size and you can draw it before fitting anything. Groups observed once or twice are where the leave-one-out posterior stops resembling the full posterior, and no amount of smoothing on the weights bridges a gap that large.

On the brms::epilepsy model with an observation-level overdispersion effect, PSIS-LOO fails on 97 of 236 folds. Moment matching, the standard repair, leaves 37 of them uncorrected and scores an elpd RMSE of 0.658 against exact refits, against raw PSIS-LOO's 0.638. The repair bought nothing. Marginalising the random-effect block and importance-sampling only the base parameters scored 0.041 and reproduced an 82-minute exact refit at no optimisation cost. The block here is Poisson, so that marginalisation is a one-dimensional quadrature rather than the closed-form downdate a Gaussian block would get.

Then the part that reaches a decision. Comparing that model against a negative-binomial alternative, default PSIS-LOO reports Δelpd = +39.3 ± 8.0, a z of 4.9 and a verdict anyone would call settled. Running reloo at the recommended k̂ > 0.7 cutoff brings it to +11.0 ± 3.2. The exact fold-wise answer is +3.1 ± 3.1, a z of 1.0 and no verdict at all. reloo lands optimistic because PSIS is already biased on the folds between 0.55 and 0.7 that it declines to refit.

Read the caveats before you act on it. That paper is a July 2026 single-author preprint, its package is R-only, and the author is explicit that the marginalisation itself is prior art (Vehtari 2016, Merkle 2019, Bürkner 2021) rather than a new idea. What the preprint contributes is the diagnostic map and the head-to-head evidence, and both are worth knowing whatever happens to the package.

A geo MMM with per-geo channel coefficients and a modest number of weeks per market is the singleton-heavy regime this describes. "LOO prefers the hierarchical specification" is a sentence MMM reports produce routinely. The sharp part is the asymmetry: the bias does not cancel when you difference two models that differ in latent structure, and pooled against hierarchical, or with against without a latent mediator, is precisely the comparison an MMM makes.

How Much of the Gap Is One Week

Back to the +8.8 nat row-wise gap in the table, and where it came from. Its largest single week contributes +2.74, which is 31% of the total, and the three largest contribute 81% between them. Delete only the largest and the gap falls to +6.0 ± 4.3, moving z from +1.73 to +1.39, while thirty of the 52 scored weeks were already pointing the opposite way from the total they sum to.

A handful of weeks are choosing the model.

Two recent results say to take that seriously. Sivula, Magnusson, Alonzo Matamoros and Vehtari (2025) study when the normal approximation behind loo_compare's standard error is well calibrated, and offer a probability-that-one-model-is-better framing in place of a bare z. Bindoff, whose own decision flip is above, describes the same standard error as "itself uncertain and mildly anti-conservative" and reads his z-values as descriptive only. Small n with near-identical candidates is the awkward case. An MMM comparison is 104 or 156 observations between two specifications that differ by one channel's transform, which is that case exactly.

The other result attacks the score itself. Choi, Spencer and Dey (2026) point out that ELPD is a log score, so an observation the model assigns near-zero predictive density to contributes an unbounded amount. Their fix is to use the same β-divergence Bregman score \(S_\varphi\) in all three places at once: for the generalized-Bayes posterior update, for forming the posterior predictive, and for the leave-one-out evaluation. That gives a generalized ELPD, with β = 1 recovering the ordinary one and β > 1 bounding what any single observation can contribute. They prove the matching is load-bearing: a procedure that evaluates with a Bregman score on top of an ordinary Bayesian update is consistent for a target that is neither the model closest under the evaluation divergence nor the one closest under the updating divergence. In their forensic footwear application the models without a contact-surface component win under the log score, driven almost entirely by one observation the contact models assign essentially zero probability, and for every β > 1 the ranking reverses.

What It Would Take Here

Refitting is the expensive half, and it is already built. run_backtest refits the full specification at every rolling origin and records one row per origin and horizon: y_true, the posterior-mean y_pred, interval bounds and coverage flags at 50/80/95, plus a naive and seasonal-naive baseline. It grades those with MAPE and MASE, which is a point-forecast question. The posterior predictive draws that produced y_pred also give the log predictive density of y_true, and summing that column over the records is an exact leave-future-out elpd with no importance sampling anywhere, because the refits already happened.

# illustrative
from mmm_framework.validation.backtest import BacktestConfig, run_backtest

result = run_backtest(model, BacktestConfig(min_train_size=104, horizon=1, step=1))
# result.records already carries y_true and y_pred per (origin, horizon).
# One more column, log p(y_true | data before the origin), turns the same
# refits into elpd_lfo, and its per-record values into a forward LOO-PIT.

The geo direction needs more: a cell-subsetting counterpart to _slice_panel_prefix so a fold can hold out a market rather than a suffix, and a decision about whether the held-out market's own intercept is drawn from the population or held fixed, which is the difference between asking about a new market and asking about a market you stopped measuring. That choice changes the estimand, so it belongs in the pre-registered spec rather than in the diagnostic.

Cost is the old objection to all of this. Cooper, Vehtari, Forbes, Simpson and Kennedy (2024) weaken it: running every cross-validation fold as a parallel chain reaches roughly the wall-clock of a single full-data fit on accelerator hardware, with partitioning schemes for non-exchangeable data and support for a range of scoring rules. The approximation inside leave-future-out exists because exact refits were unaffordable. On a GPU with 52 folds running as 52 chains, that premise is weaker than it was in 2020.

Fixing the fold makes a predictive score honest about prediction. Ranking causal estimands is beyond any predictive score, and Stop Validating Your MMM With Holdout Error is the argument for why. The simulation here is a small illustration of that too: the model row-wise LOO preferred reads total media effect 16.3% above the planted truth, and the one it rejected reads 6.9% below. Which specification to bill a channel against is still decided by recovery, refutation, coverage and a held-out experiment.

References

• Gelman, A., Vehtari, A., & McElreath, R., with Simpson, D., Margossian, C. C., Yao, Y., Kennedy, L., Gabry, J., Bürkner, P.-C., Modrák, M., & Leos Barajas, V. (2026). Bayesian Workflow. CRC Press. ISBN 978-0-367-49014-0.
• 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. Quotations here are from §6.2 of this preprint.
• 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.
• Vehtari, A., Simpson, D., Gelman, A., Yao, Y., & Gabry, J. (2024). Pareto Smoothed Importance Sampling. Journal of Machine Learning Research, 25(72), 1–58.
• 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.
• Bindoff, A. D. (2026). Partial pooling predicts cross-validation reliability: a closed-form triage and Rao–Blackwellised cure for hierarchical LOO. arXiv:2607.18836.
• Sivula, T., Magnusson, M., Alonzo Matamoros, A., & Vehtari, A. (2025). Uncertainty in Bayesian leave-one-out cross-validation based model comparison. Bayesian Analysis. arXiv:2008.10296.
• Choi, J., Spencer, N. A., & Dey, D. K. (2026). Robust Bayesian Predictive Model Selection using Bregman Divergence. arXiv:2606.10409.
• Cooper, A., Vehtari, A., Forbes, C., Simpson, D., & Kennedy, L. (2024). Bayesian cross-validation by parallel Markov chain Monte Carlo. Statistics and Computing, 34. arXiv:2310.07002.