Hundreds of Models, One Winner: Multiple Comparisons and the Winner's Curse
Fit enough model variants and one of them will look spectacular. The uncomfortable part is that its measured performance — the tight fit, the significant channel effects, the flattering ROI — is inflated by the very act of selecting it. This post explains why that happens, why the classical p-value corrections don't really fix it, and why the durable answer is to shrink the estimates before you compare them.
The workflow is familiar to anyone who has built a media mix model under deadline: four adstock forms, three saturation curves, a handful of control sets, two trend specifications — a grid of hundreds of runs. Somebody scores them, sorts the table, and the top row becomes "the model." The same pattern recurs at every scale of measurement: the best creative out of forty test cells, the strongest lift out of twelve channels, the highest measured ROI out of two hundred geos.
In every case the selection operates on noisy estimates, and the maximum of a set of noisy estimates is biased upward as an estimate of the winner's true value. Auction theorists call this the winner's curse: the bidder who wins is, on average, the one who most overestimated the object's worth. Statistics has an exact analogue, worked out in detail by Gelman and collaborators; its pieces — the garden of forking paths, Type S and Type M errors, and partial pooling — assemble into a practical answer for model builders.
The Winner's Curse
Suppose each of \( m \) candidate models (or channels, or test cells) has some true quality \( \tau_j \), and the analysis produces an estimate \( \hat{\tau}_j = \tau_j + \varepsilon_j \) with noise \( \varepsilon_j \). Ranking by \( \hat{\tau}_j \) and taking the top row does not select the model with the largest \( \tau_j \); it selects the model with the largest sum of true quality and luck. Conditional on being the winner, the noise term is almost surely positive — and the more candidates you screened, the larger the winner's luck component is likely to be.
This is not a subtle bias. When the true qualities are close together — exactly the situation when comparing hundreds of minor variations on one model — the ranking is driven almost entirely by noise, and the winner's measured edge is close to pure artifact. Its performance will not replicate, its "significant" coefficients carry no evidential weight, and decisions calibrated to its point estimates are calibrated to luck.
Pick the winner, watch it regress
Every candidate spec estimates the same true effect (the dotted line), differing only by luck. Ranking picks the largest estimate — the winner (★) — which is almost always sitting above the truth. Re-run the winning spec on fresh data and its replication value (◇) falls back toward the line. Drag the number of specs you screened and hit Resample: more candidates, higher winner, bigger the drop.
The winner's selected estimate is biased upward — it is the largest sum of true quality and luck. Its true value is just τ, so on replication it regresses by the whole gap between the star and the line. The exaggeration is manufactured by the selection, not the estimator.
Naming the problem "multiple comparisons" makes it sound like a technicality about p-values. It is better understood as a statement about estimation: selection converts honest estimates into dishonest ones. Any fix that leaves the estimates untouched and only adjusts the significance threshold is treating the symptom.
The expected selection bias — how far above the truth the winner sits, in units of the estimation noise σ — climbs relentlessly with the number of candidates screened. It is the expected maximum of m mean-zero errors: near 2σ by a hundred specs and still rising past three hundred. Screen more, and the winner you crown is exaggerated by more.
The Garden of Forking Paths
How much multiplicity does a model search actually contain? The raw arithmetic is bad enough: with \( m \) independent tests each run at level \( \alpha \),
$$ P(\text{at least one false positive}) \;=\; 1 - (1 - \alpha)^m , $$which at \( \alpha = 0.05 \) reaches 99.4% for a hundred models and is effectively certain for a thousand. But Gelman and Loken (2013) argue the count of tests you ran understates the problem. They distinguish three levels of data-contingent analysis:
- A pre-specified test, \( T(y; \phi) \) with the procedure \( \phi \) fixed before seeing data. Here the classical calibration holds.
- Researcher degrees of freedom, \( T(y; \phi(y)) \): one analysis is performed, but its choices — transformation, controls, outlier handling — were made after seeing the data. Different data would have produced a different test, so the multiplicity is real even though only one model was ever fit.
- Explicit search, \( T(y; \phi^{\text{best}}(y)) \): many analyses are run and the best is reported.
Selecting specifications by fit or significance is structurally case 3, and no bad intent is required — it is what a diligent analyst does when the workflow is "iterate until it looks right." Worse, the effective model space is the combinatorial one: every variable, transformation, and interaction that could have been tried. The hundreds of models you fit sample a space of millions you implicitly searched. A p-value computed as if the winning specification had been pre-specified answers a question nobody asked.
Type S and Type M Errors
The classical error framework — Type 1 and Type 2 — is a poor lens for this problem, because it is organized around a null hypothesis that is rarely believed. In marketing measurement, no channel's true effect is exactly zero; the practical questions are whether we have the direction right and whether the magnitude is anywhere near the truth. Gelman and Tuerlinckx (2000) reframed the errors accordingly, and Gelman and Carlin (2014) built a design-analysis workflow around them:
Definition: Type S and Type M errors
A Type S (sign) error occurs when the analysis claims an effect is in one direction — \( \hat{\tau} \gt 0 \), or \( \hat{\tau}_j \gt \hat{\tau}_k \) for a comparison — when the true effect runs the other way.
A Type M (magnitude) error occurs when the estimated size is far from the true size. Its natural measure is the exaggeration ratio
$$ \text{exaggeration ratio} \;=\; \frac{|\hat{\tau}|}{|\tau|}, $$the factor by which the reported effect overstates (or understates) the real one.
The connection to low power is the key intuition. A result is declared significant when \( |\hat{\tau}| \gt 1.96\,\mathrm{se} \). If the true effect is small relative to the standard error, the only estimates that clear that bar are the ones that landed far from the truth — the significance filter selects for overestimates. Compare two estimators when the true effect is near zero: with a standard deviation of 1, the significant estimates are modest; with a standard deviation of 3, the estimates that clear the bar are far larger in absolute value, giving the impression of a big effect. Large significant estimates from noisy analyses are a byproduct of the noise, not evidence of large effects. And with a wide sampling distribution straddling zero, the sign of a significant estimate approaches a coin flip — Type S errors climb alongside Type M.
⚠️ The winner is the biggest Type M error
Combine the significance filter with a search over hundreds of specifications and the biases compound: the specification crowned "best" is disproportionately the one whose noise was most flattering. The most significant result from a heavily searched space is not the strongest finding — in expectation, it is the most exaggerated one.
This matters operationally because subgroup cuts are where MMM insight is supposed to come from. Splitting a national read into channel × region × audience cells shrinks the sample behind each estimate and inflates each standard error — raising Type S and Type M rates at exactly the granularity where reallocation decisions get made.
Treating the Symptom: Classical Corrections
The textbook response to multiplicity is to adjust the significance threshold:
| Method | Rule at \( m = 100 \) | What it controls |
|---|---|---|
| Bonferroni | test each at \( p \le 0.0005 \) | Familywise error rate (any false positive) |
| Holm | step-down variant, slightly more powerful | Familywise error rate |
| Benjamini–Hochberg | largest \( k \) with \( p_{(k)} \le (k/m)\,\alpha \) | False discovery rate (share of false positives) |
These procedures have legitimate homes — Bonferroni and Holm for a few high-stakes pre-specified tests, Benjamini–Hochberg (1995) for thousands of parallel tests with a follow-up screen, as in genomics — but for model selection they fail on three counts.
First, they are threshold adjustments only. The point estimates — the quantities the Type M framing says are the actual problem — pass through untouched. A Bonferroni-corrected analysis of the winning model reports the same exaggerated ROI, just with a wider interval around it. Second, they assume pre-specified tests — precisely what a model search violates; there is no honest \( m \) to plug in when the implicit space is combinatorial. Third, a stricter threshold cuts power further, and low power is the engine of Type S and Type M errors — a classical correction can make surviving estimates less trustworthy in sign and magnitude even as false positives get rarer. In the Infant Health and Development Program subgroup analysis discussed by Gelman, Hill and Yajima (2012), Bonferroni widening pushed four classical subgroup intervals to cover zero, while the multilevel intervals for the same subgroups sat stably away from it.
Correcting the Estimates: Partial Pooling
The Bayesian move, argued at length by Gelman, Hill and Yajima (2012), is to change the estimates rather than the thresholds. Model the many quantities being compared — channel effects, geo ROIs, subgroup lifts — as draws from a common population:
$$ \bar{y}_j \mid \theta_j \sim \mathrm{N}\!\left(\theta_j,\, \sigma_{\bar{y}}^2\right), \qquad \theta_j \sim \mathrm{N}\!\left(\mu,\, \sigma_\theta^2\right), \qquad j = 1, \dots, J. $$The posterior mean for each \( \theta_j \) is a precision-weighted compromise between that group's own data \( \bar{y}_j \) and the population mean \( \mu \) — partial pooling. Estimates measured noisily are pulled hard toward the center; estimates measured precisely barely move. The consequence for comparisons is the central algebraic fact of this literature: the posterior z-score of any pairwise comparison is the classical z-score times a shrinkage factor,
$$ z_{\text{Bayes}} \;=\; \underbrace{\frac{\bar{y}_j - \bar{y}_k}{\sqrt{2}\,\sigma_{\bar{y}}}}_{\text{classical } z} \;\cdot\; \underbrace{\frac{1}{\sqrt{1 + \sigma_{\bar{y}}^2 / \sigma_\theta^2}}}_{\text{shrinkage factor}\; \lt\; 1}. $$The factor is always below one, and its size is governed by the variance ratio \( \sigma_\theta^2 / \sigma_{\bar{y}}^2 \). When the groups are genuinely similar (\( \sigma_\theta^2 \) small), shrinkage is strong and comparisons are aggressively damped — precisely the regime where the winner's curse bites hardest. When the groups genuinely differ (\( \sigma_\theta^2 \) large), the factor approaches one and inference approaches the classical answer — little correction was needed.
💡 The correction is data-adaptive — in both directions
Bonferroni applies the same blunt penalty whether the effects are indistinguishable or wildly different. Partial pooling estimates \( \sigma_\theta \) from the data and corrects only as much as the data warrant. In Gelman, Hill and Yajima's analysis of state-level NAEP test scores, where true state differences are large, the multilevel model declared more confident comparisons than an FDR-corrected classical analysis — it refuses to over-correct when real differences are well supported.
The simulation evidence is stark. Using the eight-schools structure — eight effects drawn from \( \mathrm{N}(0, 5^2) \), standard errors from Rubin (1981), all 28 pairwise comparisons, 1,000 replications — the classical analysis finds at least one significant comparison in 47% of replications, and among the comparisons it declares significant only 63% have the correct sign: a 37% Type S rate on published claims. The hierarchical analysis declares far fewer comparisons (at least one in only 5% of replications), and the claims it does make have the right sign 89% of the time. Rerun with genuinely large effects, \( \mathrm{N}(0, 10^2) \), and the model adapts: shrinkage weakens, more comparisons survive, and sign accuracy rises to 96%. Fewer claims, better claims — and no explicit correction step anywhere in the procedure.
| Property | Classical corrections (Bonferroni / FDR) | Multilevel model (partial pooling) |
|---|---|---|
| What is adjusted | Threshold or interval width | The point estimates themselves, plus intervals |
| Point estimates | Unchanged (still exaggerated) | Shrunk toward the population mean |
| Interval width | Always wider | Can be narrower than classical |
| Adapts to the data | No — fixed penalty in \( m \) | Yes — via the variance ratio |
| Type S / Type M errors | Not addressed; can worsen via power loss | Directly reduced by shrinkage |
Why does this count as a multiple-comparisons correction? Multiplicity is only a problem because the estimates are uncertain, and classical inference ignores what the estimates say about each other. A hierarchical model treats the \( J \) estimates as evidence about a shared population, so each borrows strength from the rest — and the noisiest estimates, precisely those most prone to sign and magnitude errors, are pulled hardest toward the center. Selection applied to shrunken estimates is selection applied to already-honest quantities. Greenland and Robins (1991) made the same argument in epidemiology, and the James–Stein result (Efron & Morris, 1975) is its frequentist shadow: with three or more parallel estimates, pooling dominates outright.
Deep diveDeriving the z-score shrinkage factor
In the normal-normal model, standard conjugate algebra gives each group's posterior as normal with
$$ \mathbb{E}\left[\theta_j \mid \bar{y}\right] = \frac{\mu/\sigma_\theta^2 + \bar{y}_j/\sigma_{\bar{y}}^2}{1/\sigma_\theta^2 + 1/\sigma_{\bar{y}}^2}, \qquad \operatorname{sd}\left[\theta_j \mid \bar{y}\right] = \left(\frac{1}{\sigma_\theta^2} + \frac{1}{\sigma_{\bar{y}}^2}\right)^{-1/2}. $$For the comparison \( \theta_j - \theta_k \), the shared prior mean \( \mu \) cancels and (treating the group posteriors as conditionally independent given the hyperparameters):
$$ \mathbb{E}\left[\theta_j - \theta_k \mid \bar{y}\right] = \frac{\sigma_\theta^2}{\sigma_\theta^2 + \sigma_{\bar{y}}^2}\,(\bar{y}_j - \bar{y}_k), \qquad \operatorname{sd}\left[\theta_j - \theta_k \mid \bar{y}\right] = \frac{\sqrt{2}\,\sigma_{\bar{y}}\,\sigma_\theta}{\sqrt{\sigma_\theta^2 + \sigma_{\bar{y}}^2}}. $$Dividing mean by standard deviation and factoring out the classical z-score \( (\bar{y}_j - \bar{y}_k)/(\sqrt{2}\,\sigma_{\bar{y}}) \):
$$ z_{\text{Bayes}} = \frac{\bar{y}_j - \bar{y}_k}{\sqrt{2}\,\sigma_{\bar{y}}} \cdot \frac{\sigma_\theta}{\sqrt{\sigma_\theta^2 + \sigma_{\bar{y}}^2}} = z_{\text{classical}} \cdot \frac{1}{\sqrt{1 + \sigma_{\bar{y}}^2/\sigma_\theta^2}}. $$The limits bracket the behavior: as \( \sigma_\theta \to 0 \) the factor goes to zero (identical groups — declare nothing); as \( \sigma_\theta \to \infty \) it goes to one (unrelated groups — classical inference returns). In practice \( \sigma_\theta \) is estimated from the data, which is what makes the correction adaptive. Note also that the posterior standard deviation of \( \theta_j \) is strictly smaller than the classical standard error: the model earns narrower intervals and fewer significant comparisons at once, because the point estimates moved.
Choosing Among Model Specifications
Marketing mix modeling meets this problem in two distinct shapes, and it pays to treat them separately.
Many parallel estimates inside one model. Channel effects across geos, ROIs across products, lifts across segments — these are the eight schools in different clothes, and the prescription transfers directly: give them hierarchical structure. A per-geo channel effect modeled with partial pooling toward the channel's national mean is a per-geo estimate that has already paid its multiple-comparisons tax. Ranking shrunken geo ROIs and reallocating toward the top of the list is a defensible operation; ranking unpooled ones is a winner's-curse generator.
Many candidate specifications. The harder shape is the one this post opened with: hundreds of model variants, one winner. Here the honest answer is to stop framing it as a selection-by-significance problem at all.
- Replace the search with one encompassing model. Most of what the grid searches over — which controls, which transformations — can live inside one larger model with regularizing priors that encode skepticism about extreme effects: ridge-like normal priors, or sparsity-inducing priors such as the horseshoe for large candidate sets. Shrinkage priors play the same role for specification uncertainty that partial pooling plays for group estimates; the variable selection page covers the options in the MMM setting.
- If a smaller model is needed, project it out of the full one. Projection predictive selection (Piironen & Vehtari, 2017) fits the full regularized model once, then finds the smallest submodel that preserves its predictive distribution. The selection is guided by one posterior rather than by hundreds of independent significance tests, which is what removes the overfitting-by-search.
- Compare models on out-of-sample prediction, not in-sample fit. Cross-validated predictive accuracy (in practice PSIS-LOO) measures something selection cannot trivially inflate — though comparing hundreds of candidates reintroduces selection noise, an argument for a few structurally distinct models rather than a dense grid of near-duplicates.
- When the comparison is uncertain, average instead of selecting. Stacking predictive distributions (Yao et al., 2018) sidesteps the winner's curse entirely: no single winner is crowned, so no winner's luck is baked into the reported numbers.
- If you do explore a grid, report the grid. A multiverse analysis (Steegen et al., 2016) fits all plausible specifications and shows how the conclusion moves across them. If the ROI ranking is stable across the model space, the specification choice was never the risk; if it swings, that instability is the finding.
Underneath all five recommendations is a discipline point rather than a technique point: decide the analysis before the data decides it for you. Pre-specifying the model space, the priors, and the decision rule — and iterating transparently when the data forces a change — is the multiplicity control that no formula supplies. That discipline, and the diagnostics that support honest iteration, are the subject of the Bayesian workflow page; the forking-paths argument is ultimately an argument for workflow, not for a different correction factor.
Takeaways
- Selecting the best of many noisy estimates inflates the winner's measured performance — and the bias grows with the number of candidates and the noise.
- Multiplicity is broader than the tests you ran: data-contingent choices (the garden of forking paths) create multiplicity even in a single fitted model.
- Think in Type S (sign) and Type M (magnitude) errors, not Type 1: under low power, the significance filter selects for exaggerated estimates, and searches select for the most exaggerated of all.
- Bonferroni and FDR adjust thresholds but leave the exaggerated point estimates untouched — and by cutting power they can make surviving claims worse in sign and magnitude.
- Partial pooling shrinks estimates toward each other adaptively, so selection operates on already-honest quantities; for specification search, prefer one regularized encompassing model, projection predictive selection, stacking, or a reported multiverse over crowning a single significant winner.
References
- Gelman, A., & Carlin, J. (2014). Beyond Power Calculations: Assessing Type S (Sign) and Type M (Magnitude) Errors. Perspectives on Psychological Science, 9(6).
- Gelman, A., Hill, J., & Yajima, M. (2012). Why We (Usually) Don't Have to Worry About Multiple Comparisons. Journal of Research on Educational Effectiveness, 5(2).
- Gelman, A., & Loken, E. (2013). The Garden of Forking Paths: Why Multiple Comparisons Can Be a Problem, Even When There Is No "Fishing Expedition" or "p-Hacking" and the Research Hypothesis Was Posited Ahead of Time. Department of Statistics, Columbia University.
- Gelman, A., & Tuerlinckx, F. (2000). Type S Error Rates for Classical and Bayesian Single and Multiple Comparison Procedures. Computational Statistics, 15(3).
- Benjamini, Y., & Hochberg, Y. (1995). Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. Journal of the Royal Statistical Society, Series B, 57(1).
- Rubin, D. B. (1981). Estimation in Parallel Randomized Experiments. Journal of Educational Statistics, 6(4).
- Efron, B., & Morris, C. (1975). Data Analysis Using Stein's Estimator and Its Generalizations. Journal of the American Statistical Association, 70(350).
- Greenland, S., & Robins, J. M. (1991). Empirical-Bayes Adjustments for Multiple Comparisons Are Sometimes Useful. Epidemiology, 2(4).
- Piironen, J., & Vehtari, A. (2017). Comparison of Bayesian Predictive Methods for Model Selection. Statistics and Computing, 27(3).
- Yao, Y., Vehtari, A., Simpson, D., & Gelman, A. (2018). Using Stacking to Average Bayesian Predictive Distributions. Bayesian Analysis, 13(3).
- Steegen, S., Tuerlinckx, F., Gelman, A., & Vanpaemel, W. (2016). Increasing Transparency Through a Multiverse Analysis. Perspectives on Psychological Science, 11(5).
- 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.