Fitting Many Models Is Not P-Hacking

Matthew Reda Bayesian workflow

This blog has argued for pre-specification and it has argued for the Bayesian workflow, and it has never put the two in the same room. The workflow says you will fit many models to one problem and that this is the point of it. Pre-specification says you commit to the analysis before you see the answer. What makes them compatible is narrower than it first looks: the estimand is fixed before any fit, the model space is declared before any fit, and no step that moves you between models is allowed to read the estimand. Gelman and his co-authors are unusually blunt about the exposure that still leaves. The one they name is subtler than picking a winner: fit a dozen defensible specifications, report the range they span, and treat that range as your uncertainty. On eight synthetic worlds where the causal answer is known by construction, the range spanned by sixteen defensible specifications contained the true ROAS in 24 of 32 channel-world pairs, and the width of the range carried no usable information about how far its centre sat from the truth.

None of that is an argument against iterating. Iteration needs a shape, and the shape has a name in the canon that this framework's own machinery already half implements without labelling it.

Two Rules That Look Opposite

The stated purpose of this framework, written at the top of its own developer guide, is “pre-specified analyses to reduce researcher degrees of freedom.” That sentence is doing real work. It explains the estimand registry, the registered_at timestamp on SpecSet, and the field-level lock list that records which parts of a spec a human touched. Read literally, it is also a promise to do the opposite of what the Bayesian workflow prescribes.

Section 8.1 of Bayesian Workflow opens by saying so out loud. “The key aspect of Bayesian workflow, which takes it beyond Bayesian data analysis, is that we are fitting many models while working on a single problem.” Many models, deliberately, as the method rather than as a confession. The same section then concedes what that costs: “we also have to be concerned with researcher degrees of freedom (Simmons et al., 2011), most directly from overfitting if a single best model is picked, or more subtly that if we are not careful, we can consider our inferences from a set of fitted models to bracket some total uncertainty, without recognizing that there are other models we could have fit.”

Read that second clause slowly. Most readers arrive expecting a different warning. The obvious sin is the one The Evidence a P-Value Actually Carries and Common Pitfalls in Statistical Modeling already cover from the frequentist side: you fit twenty things, you keep the one that flatters the channel your client cares about, and the reported uncertainty is a fiction. That failure is real and it is well policed. The subtler one is what happens after you have stopped cheating. You declare your spec set honestly, you fit every member of it, you refuse to pick a winner, you report the full spread, and what you have reported is how far apart the models you happened to think of turned out to be. Four sections later the book gives the tension its own headings, “Justification of iterative model building” (§12.2) and “Model selection and overfitting” (§12.3), which is a fair signal that the authors did not consider the matter settled by assertion.

What Has to Be Fixed

The reconciliation is a matter of ordering. Pre-specification does not require that you know which model you will end up believing. It requires that the things a motivated analyst could bend are decided before the data can bend them, and there are exactly two of those, plus one prohibition that ties them together.

The condition under which iteration is honest

Declared first: the estimand. “We will measure TV” does not qualify. What qualifies is the specific functional of the posterior, on a named scale, over a named window, with a named denominator. In this framework that is a registry entry rather than a paragraph: all_builtins() returns six of them (estimands/registry.py), of which three are on by default, and each carries its own numerator, denominator, window and causal_assumptions fields.

Declared first: the model space, which is a larger object than the model. SpecVariant (validation/spec_curve.py:88) is a serializable declaration over seven named axes: adstock form, saturation form, control set, pooling level, media-prior parameterization, trend family and seasonality, with an overrides escape hatch for anything else. SpecSet wraps a list of them with a rationale and a registered_at stamp.

Never allowed: a step that moves you between models by looking at the estimand. Predictive criteria may look at held-out \( y \). Convergence diagnostics may look at the sampler. Nothing in the loop may look at TV's ROI and use it to choose.

The prohibition does most of the work, and it rules out the most tempting shortcut. Stacking weights are chosen to maximize expected predictive utility, and this framework declines to apply them to a spec curve for exactly that reason. Two specifications can forecast the KPI equally well and still split the same fitted mean between media and baseline in ways that do not resemble each other, because the decomposition is the free parameter that prediction does not price. The module's own docstring puts it in one line: “Equal weights encode what pre-registration already asserted.” Stop Validating Your MMM With Holdout Error tells the story of that retraction in detail, and it stands.

Under those conditions, fitting a hundred models is a way of understanding one estimand, and it carries no more researcher freedom than fitting one. You may explore the whole space. You may not let the answer tell you where to stand in it. The ordering does not fix everything, though, because the space you declared has an edge.

The Bracket That Is Not a Bracket

A spec curve produces a range, and a range looks like an interval, and an interval looks like uncertainty. The three are different objects, and the gap between them is measurable if you have a world whose answer you already know.

So I built one. Eight worlds from synth/dgp.py, each a weekly panel of four media channels over three years, each sharing one clean core and breaking at most one modelling assumption: the clean positive control, unobserved confounding, reverse causality, multicollinearity, carryover misspecification, saturation misspecification, drifting effectiveness, and cross-channel synergy. Ground truth in these worlds is defined the way the model itself would report an effect, as the counterfactual of zeroing a channel's spend on the noiseless structural mean, so the true ROAS and the estimated ROAS are the same functional of the same quantity.

Against each world I ran sixteen specifications: four carryover forms (geometric with retention 0.2, 0.4 and 0.6, plus a Weibull) crossed with four saturation forms (exponential at three saturation rates, plus a Hill curve), every one of them carrying a linear trend and yearly Fourier seasonality. That is a generous declared set. The framework's own default_spec_variants returns four. Every member is a specification a competent analyst would sign off on before seeing results, and none of them is the true data-generating process for any world except by accident.

How the 128 fits were done

Each spec is fit by partialling the control block out of both the KPI and the transformed media block (Frisch-Waugh), then non-negative least squares on the media block, which is the least-squares analogue of this framework's positive-only media priors. ROAS for channel \( c \) is \( \hat\beta_c \sum_t f_c(x_{ct}) / \sum_t x_{ct} \), and every transform used satisfies \( f_c(0) = 0 \), so the estimate is the same zero-out contrast the truth is defined by. Price is excluded from every design. In the confounding and reverse-causality worlds it is exactly \( 12 + 0.5\cos(2\pi t/52) \), a perfect linear combination of the intercept and the yearly cosine, and dropping it in the other six keeps the design identical across worlds. Eight worlds by sixteen specs is 128 fits, summarised into 32 channel-world pairs. The whole thing runs in under two seconds and uses no MCMC, which is the point. This is a claim about specification sets.

Three results, in order of how uncomfortable they are.

Mild first. The truth fell inside the sixteen-spec range in 24 of the 32 pairs, which is a respectable-sounding number until you ask what generated it. Nothing did. The range has no coverage property, and a differently-sized declared set on the same worlds would have produced a different figure with equal authority.

Second, and this is the one worth carrying: the width of a channel's spec-curve range tells you nothing about how far the centre of that range sits from the truth. Across the 32 pairs the Spearman correlation between range width and absolute error of the centre is \( \rho = -0.268 \) (\( p = 0.138 \)), and on relative error it is \( \rho = -0.369 \) (\( p = 0.038 \)). Both are the wrong sign for the reading everyone applies. Three pairs are degenerate, with all sixteen specs pinned at zero, and dropping them collapses the association to \( \rho = -0.068 \) (\( p = 0.726 \)) and \( \rho = -0.152 \) (\( p = 0.431 \)). Take the conservative reading and the width is simply uninformative. Median range width was 0.527 in ROAS units, median absolute error of the centre was 0.159, and in 6 of the 32 pairs the error of the centre exceeded the entire width of the range.

Third is the one I would put on the first slide.

⚠️ Sixteen specifications agreed exactly, and were exactly wrong

In the carryover-misspecification world, where true carryover is a delayed Weibull peaking six to eight weeks out with mass past 26 weeks and every declared spec truncates at eight, all sixteen specifications returned a TV ROAS of 0.000. Not close to each other. Identical, to four decimal places. The true causal ROAS was 0.713. Social behaved the same way: sixteen specs at 0.000 against a truth of 0.846. Display spanned 0.000 to 0.052 against 0.649, and Search spanned 0.031 to 0.252 against 0.917. A zero-width spec curve is what unanimity looks like, and unanimity among models that share a wrong assumption is the least informative thing a spec curve can produce.

Nothing exotic drives that. An eight-week geometric window cannot represent the channel's response, the non-negativity constraint pins the coefficient at its boundary, and every member of the set hits the same boundary for the same reason. A real MMM with a positive media prior should land in the same place, with posterior mass piled near zero rather than a point at it. The specs had nothing to disagree about, because they were wrong in the identical direction, and spread measures disagreement.

Deep diveWhy 24 of 32 is not a coverage number

A credible interval has a construction that makes its stated mass mean something: it is a set of parameter values carrying a declared share of posterior probability under a declared model. A confidence interval has a different construction with a different guarantee, over repeated sampling. The range spanned by a spec set has neither. It is the extreme order statistics of a finite collection of point estimates from models chosen by a human, and its width is a function of how spread out that human's imagination was.

Widen the declared set with two more exotic saturation forms and the range grows, with no corresponding change in what the data can support. Narrow it to the two forms your team actually uses and the range shrinks, again with no change in the evidence. Nothing in the construction ties the width to the sampling distribution of anything, which is why the correlation with error measured above has no reason to be positive and, in this sample, was not.

What the range does measure is worth reporting under its own name: how sensitive your number is to the modelling choices you were prepared to enumerate. _robustness (validation/spec_curve.py:529) already reports it as range, spread_pct and sign_stable, alongside the fraction of the declared weight clearing the ROI reference. Sign stability across a declared set makes a useful headline. Labelling that field “total uncertainty” does not.

A Trajectory, Not a Leaderboard

If the range is the wrong summary, something has to replace it, and the canon has a proposal that this framework's machinery already computes and does not display. Section 7.4 argues that a model class is a partial ordering rather than a list. “AR(1) is simpler than AR(2) which is simpler than ARMA(2,1), and MA(1) is also simpler than ARMA(2,1), but AR(1) and MA(1) are not themselves ordered.” Some pairs of models are comparable and some are siblings, and the useful structure is the ordering itself. The purpose of that structure is stated in one sentence: “Our interest here is not in averaging over models but in navigating among them.”

What makes the move from one model to the next meaningful is that parameters in neighbouring models, in the authors' phrase, “can ‘talk with each other’ in the sense of having a shared, observable meaning outside the confines of the model itself.” Causal effect estimation is one of the two examples they name. A channel's ROAS means the same thing in the model with a linear trend and in the model with a spline, which is exactly what licenses you to plot the two on the same axis. It is also the condition The Table 2 Fallacy in Media Mix Models shows can silently fail: change the conditioning set and you may have changed the question, in which case the trajectory is comparing answers to different questions and the line through them means nothing.

Section 8.1's own figure is the display that follows. One quantity of interest, plotted across the ladder from the simplest comparison to the final model. Here is that ladder on the unobserved-confounding world, where latent demand pushes both the Search and Social budgets and the sales baseline, using geometric carryover at 0.4 and exponential saturation as the primary forms.

ModelTVSearchSocialDisplay
Spend only0.2650.7960.8440.570
+ linear trend0.2640.7370.7890.508
+ yearly seasonality0.2590.6070.6120.378
+ carryover0.4270.8720.6760.435
+ saturation0.4961.0070.7650.535
+ demand control0.4550.8910.6420.626
True causal ROAS0.6820.6300.4690.668

Read Search down the column. Seasonality is the story: adding two Fourier terms drops Search's ROAS from 0.737 to 0.607, a fall of 18%, because latent demand in this world is partly seasonal and the Fourier basis absorbs the seasonal share of the back door for free. Carryover then pushes it back up to 0.872 and saturation to 1.007, because both transforms concentrate the same spend into fewer effective exposures and the coefficient rises to compensate. Nearest the truth, at 0.607 against 0.630, sits the seasonality rung. It gets there by carrying two errors of opposite sign that have not finished cancelling. A leaderboard would have crowned it. The trajectory shows why that would be a coincidence.

Now look at the last row before the truth. Adding one control the declared sixteen never varied, a noisy proxy for latent demand, moves the whole spec curve rather than one member of it. Averaged across all sixteen specs, Social's centre moves from 0.620 to 0.508 against a truth of 0.469, closing three quarters of its bias, and Search's moves from 0.913 to 0.776 against 0.630. TV moves the wrong way, from 24.0% low to 29.1% low. That is the texture of a real adjustment rather than a magic fix. The shift is 0.19 of the declared range's own width for Social and 0.14 for Search, and no diagnostic computed inside the declared set could have pointed at it, because every member shared the omission.

Vary the conditioning set, then. SpecVariant already exposes controls as one of its seven axes and run_spec_curve will happily fit it. The default grid, and most spec curves I have seen in practice, move the transforms and hold the causal structure fixed, which is backwards: transforms are where models disagree loudly and the answer moves little, and the conditioning set is where they agree quietly and the answer moves a lot.

Why I Can Say It Missed

Everything above rests on knowing the true ROAS, which in a real engagement you do not. The stage of the workflow that supplies it is §4.3, experimentation using constructed data, and it is a different stage from the two this blog has already covered. Prior predictive checks and simulation-based calibration both draw parameters from the prior, simulate data from the model being graded, and ask whether the machinery recovers what it planted. That is a correctness check on the code and the geometry. Does Your Model Work Before You Fit It? is the reference for it. By construction it cannot tell you what happens when the world is not the model, because in that exercise the world is always the model.

Section 4.3 simulates a world your model cannot represent, then prices the failure. Its worked example plants a treatment effect of +10 in a 500-student midterm-and-final simulation and then breaks the design on purpose. Make the assignment mechanism unbalanced, with treatment probability falling in the pre-test score so the weaker students are the ones treated, and the raw comparison returns −13.8 ± 1.5. Linear adjustment for the pre-test rescues it to 9.7 ± 0.8. Then make the true conditional expectation nonlinear and re-run the same linear adjustment under the same unbalanced design, and it degrades to 7.3 ± 0.9 while the balanced design still recovers 10.5 ± 0.8. Nothing there is a coding error. The estimator is fine, the sampler is fine, and the answer is 27% low because the design and the functional form conspired. Media spend that responds to demand is the same shape of problem, drawn from marketing instead of a gradebook.

The authors are explicit that this stage is where you go when calibration is guaranteed to fail. Fitting a normal likelihood to \( t_4 \) data will fail SBC by construction, and their comment on it is the sentence that best describes what §4.3 is for: “the interesting question here is, how bad will these inferences be?”

This repository has been set up for that question for a while without quite claiming it. Every factory in synth/dgp.py is labelled with the assumption it breaks, and the module carries a flag for the case that matters most here: a scenario “sets representable=False when the truth lies outside the model's hypothesis space,” so a large error is expected by construction and reported as such rather than counted as a failure. Missing is the stage itself, as a named thing a user runs: pick the assumption you are least sure of, generate the world where it is false, run your declared spec set against it, and report the size of the error rather than its presence. What comes back is a scale, and having a scale is what separates “this could be biased” from “this is biased by about a third, upward, on the demand-chasing channels.”

The Section Nobody Quotes

Section 9 of the book is titled “Modeling as software development,” runs four subsections, and has zero coverage on this blog. The phrase “version control” appears in no post in the series. That absence is what happens when a field treats its statistical content as the interesting part and its engineering as plumbing.

Version control comes before the first prior, per §9.1. Section 9.2 makes the case for testing as you go, bottom-up, and its example is deliberately unglamorous: a one-line standard-deviation helper that divides by \( n \) instead of \( n-1 \), which the authors note “can produce mysterious biases during inference.” The bug is invisible, the fit converges, every diagnostic is green, and the posterior is quietly wrong. Nothing in a Bayesian workflow catches that. A unit test does.

Section 9.3 has the teeth, because it separates two words the replication-crisis literature has fused. Reproducibility is one analysis, consistently done, by anyone with the code and the data. Replicability is a new study reaching the same conclusion. The book's argument about the first is the sharpest sentence in the section: if re-running the analysis changes the paper's outcome, then “the original results were not particularly robust.” That is a claim about the model rather than about the plumbing. It is also why the vendor-MMM reproducibility problem carries a statistical consequence, even though its cause is that the model, the preprocessing and the deck live in three places and only one of them is versioned. Fragility under re-execution and fragility under respecification are the same fragility wearing different clothes.

What this repository already does, and what it does not

Present: a formatting gate and a lint gate in CI, contract tests that pin the public API and assert the package version matches pyproject.toml (tests/test_api_contracts.py), a lean-import test that fails if a web or LLM dependency reaches a core module, and a fit-time health snapshot stamped into every agent-driven run. Field-level spec locking is there too, so a user's manual edit survives the agent: “Anything they touch becomes locked: the LLM may no longer silently overwrite it” (agents/spec_locks.py).

Missing: a single artifact that pins the git commit, the data checksum, the resolved spec after locks and overrides, the estimand registry state and the random seed, such that one command reproduces a published number. Every ingredient exists. Nothing assembles them, and the “can you re-run last quarter's model?” question is currently answered by a person rather than a hash.

Version control before the first prior means the first thing you can answer about any number is which code produced it, which is also the only thing that makes a spec curve auditable. A declared spec set that nobody can re-fit is a paragraph, and a paragraph is exactly the thing pre-specification was invented to replace.

Takeaways

  • Iterating over many models is honest when two things are declared before any fit (the estimand, on a named scale over a named window, and the space of specifications) and one thing is forbidden (any step that moves between models by reading the estimand). Predictive criteria may look at held-out \( y \); nothing may look at TV's ROI and choose.
  • The range spanned by a declared spec set is not an uncertainty interval and has no coverage property. Across 32 channel-world pairs on eight synthetic worlds with known causal truth, it contained the truth in 24, and the correlation between its width and the absolute error of its centre was \( \rho = -0.268 \) overall and \( \rho = -0.068 \) after dropping three degenerate pairs. In 6 of 32 pairs the error of the centre exceeded the whole width.
  • Unanimity across a spec set is the failure mode with no warning light. In the carryover-misspecification world all sixteen specifications returned a TV ROAS of 0.000 against a true 0.713, and 0.000 for Social against 0.846, because every member shared the wrong eight-week window and the positive-only constraint pinned them all at the same boundary.
  • The display that replaces the range is the §8.1 trajectory: one estimand plotted across a declared complexity ladder, with unordered variants drawn as siblings rather than as further rungs (§7.4's partial ordering). On the confounding world, Search's ROAS falls 18% when yearly seasonality enters, which names seasonality as the story in a way no summary interval can.
  • Vary the conditioning set, which is where the answer moves most. Adding one demand proxy moved Social's spec-curve centre from 0.620 to 0.508 against a truth of 0.469, a shift of 0.19 of the declared set's own width, and no diagnostic inside the declared set could have found it. SpecVariant already exposes controls as one of its seven axes; the default grid does not use it.
  • Two buildable next things follow directly. A named bias-pricing stage (§4.3) that runs a declared spec set against a world where a chosen assumption is false and reports the size of the error, using the representable flag synth/dgp.py already carries. And a single provenance artifact binding commit, data checksum, resolved spec, estimand state and seed, so that a published number can be re-derived by a command instead of by a person.

References