Prior Sensitivity Without Refitting Anything

Matthew Reda Bayesian workflow

Every MMM deliverable carries a paragraph about priors. It says the team tried a tighter media prior and a looser one, and nothing moved much. Writing it takes ninety seconds. The draws from the fit you already ran contain a far more demanding version of that check, and extracting it costs about three milliseconds.

What One Fit Already Contains

The posterior is proportional to prior times likelihood. Raise the prior to a power \( \alpha \) and the perturbed posterior factors into something you already have:

$$ p_\alpha(\theta \mid y) \;\propto\; p(\theta)^\alpha\, L(\theta) \;=\; p(\theta \mid y)\, \cdot\, p(\theta)^{\alpha - 1} $$

Which means the posterior under a prior raised to the power \( \alpha \) is the posterior you sampled, reweighted by \( p(\theta)^{\alpha-1} \). No new chains. The weight for draw \( s \) is \( w_s \propto \exp\!\big((\alpha-1)\log p(\theta_s)\big) \), and every term in that expression is a number you can compute from the draws in the trace object on your disk.

Kallioinen, Paananen, Bürkner and Vehtari (2024) turned that identity into a diagnostic by adding the two things it needs to be trustworthy. Pareto-smoothed importance sampling replaces the extreme tail of the weight distribution with fitted quantiles from a generalized Pareto fit, which stops a handful of draws from carrying the whole estimate. The shape parameter \( \hat{k} \) of that fit reports when the reweighting has stopped working. They also power-scale the likelihood on its own, which turns out to matter for an MMM in a way I will come back to.

I wanted to know how much of that promise survives contact with a real media model, so I measured it. The world is synth/dgp.py:438 (make_unobserved_confounding), the same generator behind the ridge-versus-Bayes comparison: 156 weeks, four channels, latent demand pushing both spend and baseline sales. One BayesianMMM with a linear trend, NUTS through NumPyro, four chains of a thousand draws. The fit took 10.8 seconds. Then I power-scaled the joint prior at eight values of \( \alpha \) from 0.5 to 2.0, once by reweighting and once by actually refitting, and compared the per-channel incremental contribution the two paths report.

\( \alpha \)TV, reweightedTV, refitLargest gap, 4 channels\( \hat{k} \)Weight ESS
0.55,5885,5501.39%0.39700
0.85,3365,2950.78%0.103,132
1.15,1025,1071.01%−0.043,832
1.54,8404,8350.58%0.171,855
2.04,5614,5431.11%0.43401

Across all thirty-two comparisons (eight values of \( \alpha \), four channels) the worst disagreement between the reweighted estimate and its refit was 1.39% of the refit's value, and no comparison moved more than 0.071 posterior standard deviations. The typical disagreement is 0.28%, against Monte Carlo error of 0.20% to 0.39% on the reference fit's own channel means. That 1.39% sits at \( \alpha = 0.5 \), where the reweighted estimate runs on a weight ESS of 700 and carries roughly 1.1% of Monte Carlo error on its own, so these gaps are draws landing in slightly different places rather than the method going wrong.

The eight refits cost 44.9 seconds of sampling. The eight reweights cost 2.59 milliseconds, on top of a 0.65-second pass to evaluate the log prior and log likelihood at the existing draws.

Beyond the size, the direction carries its own reading. Strengthening the prior pulls TV's credited contribution down by roughly a thousand units of KPI across that range, which is the same mechanism the confounding post described from the other side: on this world the media prior is quietly pushing media estimates down, and the back-door path is pushing them up. Here the pull is a number rather than an argument.

⚠️ The reweighting tells you when it has failed, and you have to look

At \( \alpha = 2 \) the effective sample size of the weights had fallen from 4,000 to 401. Push further and the estimate becomes four draws in a trench coat. The cutoff documented by this framework's own PSIS helper is 0.7, and utils/arviz_compat.py:255 says plainly that weights above it are unreliable. My sweep peaked at 0.43. Gelman et al. make the same point about importance sampling between two posteriors in §6.3, and add the part people skip: the bridge works “provided the two posteriors are similar enough for importance sampling to bridge,” and when they are not, “this is also valuable information in itself.” A prior you cannot reach by reweighting is a prior far enough from yours to deserve its own fit.

The whole recipe is about fifteen lines against PyMC and ArviZ.

# illustrative
import numpy as np
import pymc as pm
from arviz_stats.base import array_stats

pm.stats.compute_log_prior(idata, model=mmm.model)
n_draws = idata.posterior.sizes["chain"] * idata.posterior.sizes["draw"]
log_prior = sum(
    idata["log_prior"][v].values.reshape(n_draws, -1).sum(axis=-1)
    for v in idata["log_prior"].data_vars
)

def power_scale(alpha, quantity):
    """Posterior mean of `quantity` under a prior raised to `alpha`."""
    # array_stats.psislw negates its argument, so pass the ratio negated.
    log_w, khat = array_stats.psislw(-(alpha - 1.0) * log_prior, r_eff=1.0, axis=-1)
    w = np.exp(np.asarray(log_w) - np.max(log_w))
    w /= w.sum()
    return float(w @ quantity), float(khat)

compute_log_prior evaluates the prior with the value transforms removed, so the density is the one on the constrained scale and no Jacobian sneaks into the power. That detail is load-bearing. Scale the transformed-space density instead and you are power-scaling a change of variables along with the prior, which is a different and less meaningful perturbation.

Width Against Location

This framework already ships a prior-versus-posterior diagnostic. diagnostics/learning.py:288 (parameter_learning) reports contraction \( c = 1 - \mathrm{Var}_{\text{post}} / \mathrm{Var}_{\text{prior}} \), a prior-posterior overlap coefficient, and a location shift in prior standard deviations, then grades each parameter as strong, moderate, weak, relocated or prior-dominated. It is good, and honest about being heuristic. Its own docstring names the boundary of what it measures: “contraction tells you the data spoke; the posterior location tells you what it said.”

The confounding post put the complaint more bluntly. Contraction measures how much the data moved your uncertainty, not how much the prior moved your point estimate, and that post proposed the missing sibling in two sentences: refit with the prior progressively widened and watch the coefficient walk. Power-scaling is that instrument, without the refits, and with a summary statistic instead of an eyeball.

Prior elasticity

For a scalar quantity \( q \) (a coefficient, a channel's contribution, a ROI), define

$$ e_{\text{prior}}(q) \;=\; \frac{1}{\mathrm{sd}[q]} \cdot \frac{\partial\, \mathbb{E}_\alpha[q]}{\partial \log \alpha}\bigg|_{\alpha = 1} $$

estimated by a central difference between \( \alpha = 0.8 \) and \( \alpha = 1.25 \). Units are posterior standard deviations of \( q \) per unit of log prior sharpness. This is my summary, not the one priorsense publishes; power-scaling is the shared machinery underneath.

Run both diagnostics over the nineteen scalar parameters of that fit and they disagree almost completely. Rank-correlate contraction against the magnitude of prior elasticity and you get −0.121. Pearson gives −0.007. Two parameters both graded strong by contraction sit thirteenfold apart on how much their location leans on the prior.

ParameterContractionVerdictPrior elasticity
sat_lam_Display0.080prior-dominated−0.613
adstock_alpha_Display0.115moderate−0.618
beta_TV0.719strong−0.486
adstock_alpha_TV0.846strong−0.652
beta_Display0.934strong−0.049
trend_slope0.952strong−0.629

Look at the last two rows against each other. beta_Display contracted 0.934 and barely notices the prior being sharpened or flattened. trend_slope contracted 0.952, an even tighter posterior by the width reading, and its location moves 0.63 posterior standard deviations per unit of log prior sharpness. Both are strong. Only one of them is a number the data placed.

None of that makes contraction wrong. It answers the question it was built to answer, which is whether the data narrowed the parameter, and the docstring is explicit that this is informativeness rather than importance. The two diagnostics are close to orthogonal on this fit, which is the useful discovery: a delivery gate that wants to stamp a headline number as prior-dominated needs both, because a parameter can be pinned tight and still be sitting exactly where its prior put it.

Scale the Likelihood Too

Power-scaling the likelihood by the same trick costs nothing extra once you have the pointwise log likelihood, and it turns one reading into two. Kallioinen et al. describe the joint pattern. A parameter sensitive to both prior and likelihood scaling is showing prior-data conflict, two sources of information pulling against each other. A parameter sensitive to prior scaling and insensitive to likelihood scaling is showing something else entirely: the likelihood had little to say, so the posterior tracks the prior because there is nothing else for it to track.

That second pattern is the MMM's native failure mode, and it sorts this fit's parameters cleanly.

ParameterPrior elasticityLikelihood elasticityRatio
sat_lam_Social−0.3880.0785.0
sat_lam_Display−0.6130.2053.0
sat_lam_TV−0.238−0.0862.8
sat_lam_Search−0.149−0.0662.2
beta_TV−0.4860.4681.0
adstock_alpha_TV−0.6520.8940.7

All four saturation parameters land above 2. Every adstock rate and every media coefficient with an elasticity worth reading lands below 2. beta_Display posts 11.1 on this scale by dividing a prior elasticity of −0.049 by a likelihood elasticity of 0.004, which is a ratio of two quantities near zero rather than a coefficient the prior is placing. On this world the saturation block is being placed by its prior because the likelihood is close to flat along it, and the diagnostic says so from one fit, per parameter, with a number attached. That these transform parameters are weakly identified in an MMM is the carryover and shape post's result, established on a simulation with a planted truth to check against. The two-axis reading is how you detect the same thing on a client's data, where there is none.

The Robustness Paragraph Checks the Wrong Block

Nothing forces you to scale the whole prior at once. Split the log prior into blocks, scale each block on its own, and you get a decomposition of where the headline number's dependence on assumption actually lives. For TV's incremental contribution on this fit, the joint prior elasticity is −0.835 and it splits like this.

Prior blockElasticity of TV contributionShare
baseline (intercept, trend, seasonality)−0.41549%
adstock rates−0.21425%
saturation−0.13716%
media coefficients−0.0739%
control coefficients−0.0071%
observation noise0.0040%

The six blocks sum to −0.842 against the joint −0.835, which is what local additivity is supposed to look like and a decent check that the arithmetic is sound. I verified two of the rows the slow way: power-scaling only the media coefficient priors at \( \alpha = 2 \) and refitting gives 5,090.7 for TV against the reweighted 5,119.8, and doing the same to the baseline block gives 4,844.4 against 4,856.0. Both inside 0.6%.

Now read the table as an audit of standard practice. The robustness paragraph in the deliverable varies the media ROI prior, because that is the prior anyone argues about in a kickoff meeting. On this fit that block accounts for 9% of the sensitivity of the number being defended. Half of it sits in the baseline block, which is the trend and seasonality prior nobody discusses, and a further quarter sits in the adstock rates. An analyst who ran the conventional check and reported stability would have been telling the truth about the axis they varied.

Four Points Are Not a Neighbourhood

There is a second problem with the robustness paragraph, and it is geometric. The check varies one prior hyperparameter at a time, holding the others at their defaults, then reports the range of the answer over those points. That procedure has no claim on the worst case, because the worst case does not have to lie on an axis.

Odnoblyudova, Dellaporta and Briol (2026) show how little a finite sweep guarantees even along one axis. On a telecom model the sensitivity is non-convex in a prior correlation parameter, so evaluating four candidate values and taking the worst selects a local maximum at 0.05 while the true supremum sits at 0.2 and is never visited. Their contribution is the machinery that makes the global question tractable. A Fisher divergence between the reference posterior's score and a candidate's never requires sampling from the candidate, only evaluating that candidate's score at the reference draws, and it costs \( O(md) \) in draws and dimension. For exponential-family priors with a loss linear in the hyperparameters it reduces further to a convex quadratic, so the supremum over a hypercube neighbourhood can be read off a handful of vertices.

Reweighting gets at a smaller version of the same point without the machinery, provided the neighbourhood stays close enough for \( \hat{k} \) to allow it. The media coefficient prior in this framework is \( \mathrm{Gamma}(\mu = 1.5, \sigma = 1) \) on the standardized coefficient. I swept a rectangle around it, \( \mu \) from 1.0 to 2.0 and \( \sigma \) from 0.6 to 1.5, at 41 by 37 resolution. That is 1,517 candidate priors, evaluated in 1.42 seconds, with a maximum \( \hat{k} \) of 0.547 across the whole grid. Reweighting back to the reference prior reproduced the fit's own channel means to within \( 3 \times 10^{-13} \) of a KPI unit, which is the check that the weights do nothing when they should do nothing.

The conventional sweep of five points (the default, a tighter \( \sigma \), a looser \( \sigma \), a lower \( \mu \), a higher \( \mu \)) moves total media contribution over a range of 207 units. Over the full rectangle it moves 361, with both extremes at corners the axis-wise sweep never visits: 15,275.6 at \( (\mu = 1.0,\, \sigma = 1.5) \) and 15,636.6 at \( (\mu = 2.0,\, \sigma = 0.6) \). So the informal check found 57% of the movement available inside the same neighbourhood it thought it was searching.

For one channel it gets much worse. Search moves 5.1 units across the five axis points and 50.1 across the rectangle, so the conventional sweep reported a tenth of the sensitivity that a joint search over the same hyperparameters exposes. I refit at both corners to confirm the reweighting was not inventing it. At \( (\mu = 1.0,\, \sigma = 1.5) \) the largest per-channel disagreement was 0.58%, and at \( (\mu = 2.0,\, \sigma = 0.6) \) it was 1.15%.

These are small moves. Search's fifty units sit against a 90% credible interval 2,192 units wide. I chose a defensible neighbourhood and the answer inside it is stable, which is a real finding about this fit and worth reporting as one. My claim is about the procedure. A five-point axis-wise sweep understated the movement inside its own neighbourhood by a factor of ten on one channel, and the only reason I know that is that scanning fifteen hundred alternatives cost a quarter of one refit.

What This Does Not Establish

Power-scaling is a local sensitivity analysis. It perturbs infinitesimally around the prior you actually used and reports the gradient, so it answers whether the posterior is leaning on that prior. It does not answer how bad a defensible alternative prior could make things, and Odnoblyudova et al. classify it that way explicitly when they motivate the global version. A parameter with zero prior elasticity at \( \alpha = 1 \) can still move a long way under a prior centred somewhere else entirely.

Somebody still has to pick the neighbourhood. My rectangle was mine, and a different analyst's rectangle would produce a different worst case. The global method removes the search but not the specification.

None of this touches confounding. The back-door path on this world is open, and the fit is badly wrong because of it: it credits Search with 5,126 against a planted truth of 3,654. A prior sensitivity analysis reports how much of that wrong answer came from the prior. It does not make the answer right. The prior post made the stronger version of that argument: a prior tight enough to offset a back-door path is also tight enough to override real evidence, and the route to closing the path is an experiment.

And a low \( \hat{k} \) is necessary rather than sufficient. It certifies that the weights have a finite-variance tail, not that the perturbed posterior is unimodal or that the sampler explored it in the first place. If the reference fit had the pathologies in the sampling diagnostics post, reweighting inherits every one of them and reports a clean \( \hat{k} \) while doing it.

This is one fit on one synthetic world, with one seed. The mechanism generalizes because it is an identity about densities. The particular numbers do not.

What to Build

Nothing in this framework computes any of it. The closest thing is tests/test_continuous_learning.py:978, a @pytest.mark.slow audit that fits the bandit model twice at two values of one prior scale and asserts the posterior spread grows. That test is correct and it costs two fits to interrogate one hyperparameter. The 1,517-cell sweep above cost 1.42 seconds and no fits at all.

A diagnostics/prior_sensitivity.py sitting beside learning.py starts with one table: per-parameter elasticity on both axes, prior and likelihood, with \( \hat{k} \) and weight ESS beside every row so a reader can see which cells to distrust. The block decomposition comes next, keyed to a named quantity rather than to a raw coefficient, because what anyone asks an MMM about is a channel's contribution or its ROI. Then the neighbourhood scan, which takes hyperparameter ranges and evaluates the grid. Where \( \hat{k} \) exceeds 0.7 it refuses to report the cell and prints the refit command for it instead.

Reporting is the part that changes what a client sees. The acceptance gate demands an evidence tier on every headline number, and prior-dominated is one of the tiers it can assign. Contraction cannot assign it, for the reason the middle of this post measures. Prior elasticity against a stated threshold can, and it comes with a sentence a client can act on: strengthening or flattening your assumptions by a factor of two moves this channel's credited contribution by a third of its reported interval, and here is which assumption did it.

Takeaways

  • The posterior under a power-scaled prior is the posterior you already sampled, reweighted by \( p(\theta)^{\alpha-1} \). On a 156-week four-channel MMM it reproduced eight refits to within 1.39% on per-channel contribution, at 2.59 milliseconds against 44.9 seconds.
  • The Pareto \( \hat{k} \) of the weights is the part that makes it safe to publish. It reached 0.43 at \( \alpha = 2 \) here, and a value past 0.7 is an instruction to refit rather than a number to interpret.
  • Contraction and prior elasticity measure different things and barely correlate (Spearman −0.121 over nineteen parameters). Two parameters graded strong by contraction differed thirteenfold in how much their location leaned on the prior.
  • Scaling the likelihood on its own tells prior-data conflict apart from a likelihood that had nothing to say. All four saturation parameters on this fit were more than twice as prior-sensitive as likelihood-sensitive, and no adstock rate or media coefficient with an elasticity worth reading was.
  • Blocking the prior showed the baseline block carrying 49% of TV's sensitivity against the media coefficient block's 9%, which inverts the priority of the robustness paragraph everyone writes.
  • A five-point axis-wise sweep of the media prior captured 57% of the movement available in its own rectangle, and 10% of it for one channel. Scanning 1,517 alternatives instead cost 1.42 seconds.

References