When Your MMM Won't Mix: Sweep the Priors Before You Refit Everything

A PyMC-Marketing 1.0 case study with 60 SKU × customer cells, realistic commercial confounding and one expensive validation fit

A hierarchical marketing mix model can fail in a particularly frustrating way: the code runs, the sampler finishes, the notebook contains posterior summaries, and none of that means the model has actually learned the posterior you intended to fit.

That is the situation in this case study.

We are modelling a realistic soft-drinks portfolio with three brands, twelve SKUs, five retail customers, four media channels, price, promotion, distribution and three years of weekly data. The business question is not merely, "does TV work?" We want a model that can estimate media contribution while also telling us whether a 500 ml zero-sugar SKU at an online retailer is more price sensitive than the same brand in a national grocer.

That richer model is useful. It is also exactly the sort of model where weak identification, correlated commercial drivers and hierarchical geometry can make NUTS miserable.

The first fit in this example has 187 divergent transitions, a maximum rank R-hat of about 2.50, a minimum bulk ESS below 5, and two SKU/customer price elasticities with the wrong sign. More draws would give us a larger NetCDF file, not a more trustworthy answer.

The repair workflow is:

1. understand the data before blaming the sampler;
2. make the broad model sample correctly using non-centred parameterisation;
3. use that trustworthy reference posterior to screen plausible prior changes with pymc-marketing-sweep;
4. reject candidate priors that the reference posterior cannot support;
5. protect business quantities such as channel contribution and price elasticity from excessive movement;
6. run one real PyMC-Marketing validation refit with the recommended prior configuration.

The result is a model with clean chain mixing and price elasticities that line up across SKU, customer and brand without forcing every cell to be identical.

Reproducibility note: all business data are synthetic. The ZIP includes deterministic posterior fixtures so every chart and sweep result renders immediately. The fixture is deliberately constructed to mimic realistic convergence failure, but it is not being presented as a real MCMC run. The full PyMC-Marketing 1.x scripts are included so you can replace the fixture with actual NUTS fits. The PSIS prior sweep itself is genuinely executed by the bundled sweep engine against the well-sampled synthetic reference posterior.

First, what business are we actually modelling?

Before talking about R-hat, it is worth looking at the thing the model is trying to explain.

The synthetic portfolio contains:

• 3 brands: Aster, Brio and Cinder
• 4 variants per brand: Original 330 ml, Zero 330 ml, Original 500 ml and Zero 500 ml
• 5 customers: NorthMart, ValueHub, CityExpress, FreshBasket and eGrocer
• 156 weeks from January 2023 to December 2025
• 12 SKUs and 60 SKU/customer cells
• 9,360 panel rows
• 23.1 million units over the full sample
• roughly 33.7 million in synthetic revenue
• around 148,000 units per week across the portfolio
• 4 media channels: TV, paid social, paid search and retailer media
• 3 commercial controls: price, promotion and weighted distribution

The footprint is deliberately uneven. NorthMart is larger than CityExpress. Aster is larger than Brio. That matters because a model that pretends every customer and brand combination has the same amount of information is already starting from the wrong place.

Figure 1. Total unit volume by customer and brand. The model has to learn across cells with different commercial scale and information content.

The weekly portfolio is not stationary either. Demand moves over time, promotions arrive in bursts, national media pulses across brands, list price trends upwards and Zero variants gradually reach broader distribution.

Figure 2. Weekly portfolio volume and the share of SKU/customer cells on promotion. This is intentionally closer to retail data than a set of independent Gaussian regressors.

Zoom into one cell

Take Aster Zero 500 at eGrocer. This one time series already contains most of the identification problem in miniature.

Figure 3. Aster Zero 500 at eGrocer. Promotion weeks coincide with lower price and large movements in demand.

Across the complete panel, roughly 24.7% of cell-weeks are on an active promotion. The median active promotion depth is about 14.7%. Zero variants ramp from average weighted distribution of roughly 57% in 2023 to 84% in 2025.

Those are useful business features. They are also sources of statistical confounding.

Price is almost never changed in isolation

The model uses:

price_discount_signal = -log(actual_price / list_price)

so a larger value means a deeper effective discount. Because the target is log_units, a positive coefficient on this signal maps naturally to a negative conventional price elasticity.

But in retail data, price changes rarely happen alone. A retailer lowers price because a promotion is running. Retail media is often increased for the same event. Search demand reacts. Distribution may also change around launches and promotional windows.

The synthetic data make that problem explicit. The correlation between the engineered discount signal and promotion depth is about 0.993. The correlation between discount signal and retailer media is about 0.70.

Figure 4. Left: promotion depth and the price signal are nearly the same commercial event. Right: retailer media also rises when products are discounted.

This is not a bug in the data generator. It is the point of the exercise. A useful MMM has to decide how much demand movement belongs to price, promotion, retailer media, distribution and national media when those things were planned together.

We genuinely need hierarchical variation

The answer is not to collapse everything into one global price coefficient.

The synthetic truth deliberately contains economically plausible heterogeneity. eGrocer is generally more price sensitive than NorthMart, and some 500 ml or zero-sugar variants are more elastic than their 330 ml counterparts. The true elasticities span roughly -1.24 to -2.48.

Figure 5. The true price-elasticity surface across 60 SKU/customer cells. We want partial pooling, not complete pooling.

A good hierarchical model should borrow information across those cells while preserving meaningful differences. It should not let every cell wander independently, and it should not force every cell to be equal.

That is the modelling tension behind the rest of the article.

Figure 6. Correlation across the commercial drivers. Price is the clearest problem, but it is not the only one.

The model we actually want

The model is a multidimensional PyMC-Marketing 1.x MMM with dimensions:

brand × customer × variant

Media has adstock and saturation. Commercial controls include price, promotion and distribution. Seasonality is included. The control coefficients and media-response parameters are partially pooled through hierarchical priors.

A simplified version of the broad control hierarchy is:

from pymc_extras.prior import Prior


gamma_control = Prior(
    "Normal",
    mu=0.25,
    sigma=Prior(
        "HalfNormal",
        sigma=1.0,
        dims=("control", "brand", "variant"),
    ),
    dims=("control", "brand", "customer", "variant"),
    centered=True,
)

The media saturation hierarchy is deliberately flexible too:

from pymc_extras.prior import Prior
from pymc_marketing.special_priors import LogNormalPrior


saturation_beta = LogNormalPrior(
    mean=Prior(
        "Gamma",
        mu=0.22,
        sigma=0.18,
        dims=("channel", "brand"),
    ),
    std=Prior(
        "HalfNormal",
        sigma=0.80,
        dims=("channel", "brand"),
    ),
    dims=("channel", "brand", "customer", "variant"),
    centered=True,
)

There is nothing inherently absurd about either prior. The difficulty comes from the combination of broad hierarchical scales, dozens of weakly identified local effects, correlated controls and a centred parameterisation.

What non-convergence actually means

This term is often used too casually, so it is worth spelling out.

When we run four MCMC chains, we are asking four independent simulations to explore the same posterior distribution. They may start in different places, but after warm-up they should forget those starting points and spend their time moving through the same probability landscape.

A converged fit is not one where every chain draws the same values. It is one where the chains look like interchangeable samples from the same target distribution.

A non-converged fit is one where that has not happened.

For practical MMM work, I look at four things together:

What bad mixing looks like

Here is one price elasticity from the initial centred hierarchy.

Figure 7. Four chains for the exact same SKU/customer elasticity. They are not noisy versions of one answer. They spend long periods in different regions.

The fixture has a maximum rank R-hat of about 2.50. That is not a marginal diagnostic warning. It is evidence that the chains disagree dramatically about the distribution they are supposed to be sampling.

Another way to see the same problem is to ignore draw order and overlay each chain's marginal posterior density.

Figure 8. Per-chain posterior densities for the same elasticity. If the chains had mixed, these curves would largely lie on top of one another.

This is why posterior overlay plots are so useful. A single combined density can hide the fact that it is actually an average of several incompatible chain-specific answers.

The ESS tells the same story numerically. The weakest parameter in the bad fixture has fewer than 5 effective bulk draws. You cannot rescue that by reporting more decimal places.

Most importantly, the failure reaches the business quantities. Two SKU/customer elasticities even flip to the wrong sign, and overall elasticity RMSE against the known synthetic truth is about 0.73.

That makes this much more than a sampler-aesthetics problem.

Step 1: make the broad model sample correctly

The first intervention is not to tighten the prior.

If the reference posterior itself is poorly explored, an importance-reweighting sweep has no solid foundation. Reweighting a bad sample does not turn it into a good sample.

Instead, we keep the intended broad prior distribution and switch the hierarchical parameterisation from centred to non-centred:

gamma_control = Prior(
    "Normal",
    mu=0.25,
    sigma=Prior(
        "HalfNormal",
        sigma=1.0,
        dims=("control", "brand", "variant"),
    ),
    dims=("control", "brand", "customer", "variant"),
    centered=False,
)

Conceptually, the prior is still saying the same thing about plausible control effects. We are changing the geometry exposed to NUTS.

In the reference fixture, maximum R-hat falls to roughly 1.006 and minimum bulk ESS rises above 1,360. The chain density overlay now looks like this:

Figure 9. Same business quantity, same broad-prior target, better parameterisation. The chains now largely overlay.

The aggregate diagnostics make the contrast clear.

Figure 10. The centred fit has serious disagreement across several parts of the hierarchy. The reference and validation fits are close to one.

Figure 11. The weakest bulk ESS moves from roughly five effective draws to thousands.

At this point we have something valuable: a trustworthy broad-prior posterior that contains information about where plausible alternative priors might place their mass.

Now the sweep becomes meaningful.

Step 2: use the sweep to screen priors before paying for another full fit

A conventional iteration loop might be:

1. guess that one prior is too loose;
2. tighten it;
3. run NUTS for hours;
4. inspect convergence;
5. discover that the prior either did not help or changed the model too much;
6. repeat.

With several hierarchical scales, that quickly becomes an expensive Cartesian search.

pymc-marketing-sweep asks a cheaper question first:

Given a well-sampled posterior under the fitted prior, which candidate prior changes are still supported by those posterior draws, and which ones alter the business quantities as little as possible?

The methodology is intentionally simple at a high level.

Figure 12. The sweep is a screening layer between one trustworthy reference fit and one real validation refit.

For a fitted prior (p_0(\theta)) and candidate prior (p_1(\theta)), an existing posterior draw receives an importance ratio proportional to:

p1(theta) / p0(theta)

The package then uses Pareto-smoothed importance sampling (PSIS) to stabilise those ratios and reports diagnostics that tell us whether the approximation is trustworthy.

The key outputs are:

• Pareto-k, which diagnoses whether a few extreme weights dominate;
• reweighting ESS, which says how much effective information survives the prior change;
• reweighting ESS fraction, which makes that quantity comparable across posterior sizes;
• posterior distortion, measured on variables you explicitly want to protect;
• prior distance, used to prefer the least invasive feasible change.

It is important to keep the language precise: reweighting ESS is not predicted NUTS ESS. The sweep estimates whether the old posterior can safely approximate the candidate posterior. Only a real validation refit tells us whether NUTS now samples the candidate model well.

The grid

We target three places where the broad reference model has more flexibility than we think the data can reliably support:

1. the hierarchical scale controlling customer variation in commercial controls, including price;
2. the cross-customer/variant scale on media saturation beta;
3. the cell-level intercept prior.

The YAML is deliberately readable:

sweep:
  mode: product
  max_candidates: 100

  priors:
    price_control_pooling:
      prior_path: gamma_control.sigma
      posterior_variable: gamma_control_sigma
      parameters:
        sigma:
          values: [0.50, 0.35, 0.25]

    media_response_pooling:
      prior_path: saturation_beta.std
      posterior_variable: saturation_beta_std
      parameters:
        sigma:
          values: [0.50, 0.35, 0.25]

    baseline_regularisation:
      prior_path: intercept
      posterior_variable: intercept_contribution
      parameters:
        sigma:
          values: [1.5, 1.0]

constraints:
  max_pareto_k: auto
  min_reweight_ess: 250
  min_reweight_ess_fraction: 0.08
  max_relative_rmse: 0.05

outputs:
  comparison_variables:
    - gamma_control
    - channel_contribution_original_scale

That gives us 18 candidate prior configurations from one reference posterior.

The guardrails matter as much as the grid. We are explicitly telling the sweep that a numerically attractive candidate is not acceptable if it materially rewrites price effects or channel contribution.

The most useful candidates are sometimes the ones the sweep rejects

Tighter priors can make a model look calmer simply because they remove parameter space. That does not mean the candidate is supported by the data.

Several aggressive candidates in this example push Pareto-k above 0.8 and collapse reweighting ESS to about 40 effective draws. The worst candidate reaches Pareto-k around 0.95 with only about 0.7% of the original posterior information remaining.

The sweep rejects those candidates.

Figure 13. Feasible prior changes retain substantial posterior support. Aggressive candidates move into regions the reference draws cannot safely represent.

That refusal is a feature. Without it, a prior-search tool can become a very efficient way to manufacture overconfident answers.

The recommended change is deliberately boring

The selected fixture candidate changes the control hierarchy scale from 1.0 to 0.50, the media saturation-beta scale from 0.80 to 0.50, and the intercept standard deviation from 2.5 to 1.5.

It has:

The point is not that 0.50 is a universal magic number. The point is that this is the smallest supported changed prior that clears the configured guardrails.

The package emits a native PyMC-Marketing 1.x configuration patch so the recommendation can be validated directly in the real model:

from pymc_extras.prior import Prior
from pymc_marketing.special_priors import LogNormalPrior


recommended_model_config_patch = {
    "gamma_control": Prior(
        "Normal",
        mu=0.25,
        sigma=Prior(
            "HalfNormal",
            sigma=0.50,
            dims=("control", "brand", "variant"),
        ),
        dims=("control", "brand", "customer", "variant"),
        centered=False,
    ),
    "intercept": Prior(
        "Normal",
        mu=7.25,
        sigma=1.50,
        dims=("brand", "customer", "variant"),
    ),
    "saturation_beta": LogNormalPrior(
        mean=Prior(
            "Gamma",
            mu=0.22,
            sigma=0.18,
            dims=("channel", "brand"),
        ),
        std=Prior(
            "HalfNormal",
            sigma=0.50,
            dims=("channel", "brand"),
        ),
        dims=("channel", "brand", "customer", "variant"),
        centered=False,
    ),
}

What happens to price elasticity?

This is where the example becomes economically interesting.

The bad centred fit has price-elasticity RMSE of about 0.73 against the known synthetic truth. Some cells are simply far away, and two have the wrong sign.

Figure 14. Poor chain mixing produces bad business estimates, not merely ugly diagnostics.

The non-centred reference improves elasticity RMSE to about 0.176 because it actually samples the intended broad-prior posterior. The sweep approximation improves it slightly further to about 0.172.

That modest change is exactly what we want from a screening step. PSIS should not invent a completely different model. It should identify a plausible prior change that is worth validating.

Then we perform the only expensive refit that matters.

python scripts/fit_fixed_model.py

The final validation fixture has:

• 0 divergences
• maximum R-hat around 1.004
• minimum bulk ESS above 2,200
• price-elasticity RMSE around 0.036
• 0 wrong-sign price elasticities

Figure 15. The PSIS approximation is useful for screening, but the real validation refit is what confirms the final posterior.

The important bit is what did not happen: the hierarchy did not squash every SKU and customer onto the same number.

Figure 16. Aster across customers and variants. The final model keeps real retailer and pack-size differences while removing unsupported wandering.

You can see the whole repair path in one chart.

Figure 17. Reparameterisation fixes the broad model's geometry. The sweep chooses a supported regularisation change. The final NUTS refit confirms the result.

This is not "tight priors fix convergence"

That would be the wrong lesson.

The workflow is:

understand the commercial data
        ↓
identify genuine convergence failure
        ↓
repair parameterisation or model geometry first
        ↓
obtain a trustworthy broad-prior reference posterior
        ↓
screen plausible prior changes with PSIS
        ↓
reject poor-overlap candidates
        ↓
protect business quantities from excessive movement
        ↓
choose the smallest supported changed prior
        ↓
run one real validation fit

There are many problems a prior sweep should not be expected to repair:

• omitted-variable bias;
• incorrect causal structure;
• a misspecified likelihood;
• deterministic or near-deterministic collinearity;
• a reference posterior that never explored its own target;
• a candidate prior so different that PSIS says the old posterior has no useful overlap.

The last case is especially important. The correct output is "refit required", not "extrapolate harder".

The repository is agent-ready, including custom PyMC-Marketing models

The bundled pymc-marketing-sweep repository is not hard-coded to Aster, eGrocer, price or these exact dimensions.

It includes agent instruction files:

pymc-marketing-sweep/
├── AGENTS.md
├── CLAUDE.md
└── .agents/
    └── skills/
        └── pymc-marketing-sweep/
            └── SKILL.md

The skill tells Codex, Claude or Cursor to:

1. inspect the saved PyMC-Marketing 1.x model first;
2. discover the actual model configuration, prior paths, posterior variable names and dimensions;
3. modify the YAML sweep grid rather than rewriting the statistical engine;
4. preserve the fitted distribution and dimensions unless the user explicitly asks to change them;
5. protect model-specific posterior quantities chosen by the user;
6. emit native PyMC-Marketing Prior(...) or supported special-prior configuration for the final refit.

That is what makes the workflow useful for custom PyMC-Marketing models.

Your dimensions could be:

brand × customer × SKU

or:

country × channel

or:

DMA × product × customer_segment

The workflow is the same. The agent should inspect the actual saved model rather than assume that a variable is called saturation_beta or that a dimension is called customer.

A practical prompt can be as simple as:

Read AGENTS.md and the bundled pymc-marketing-sweep skill first.

Inspect outputs/runtime/reference_model.nc with pmm-sweep inspect.
Do not guess prior paths or posterior variable names.

The model has unstable SKU/customer price elasticities.
Find the hierarchy controlling those coefficients and propose a conservative
sweep around its fitted scale. Also inspect the media-response hierarchy.

Keep the fitted distributions and dimensions unchanged.
Reject candidates with poor PSIS overlap or low reweighting ESS.
Protect price coefficients and channel_contribution_original_scale from
material movement.

Run the sweep, explain the rejected candidates, and give me the exact
PyMC-Marketing 1.x prior patch for the least-invasive feasible candidate.
Do not claim the issue is fixed until the final real NUTS refit has acceptable
divergences, R-hat and ESS.

This is a much safer task for an agent than "keep changing priors until the diagnostics look better".

Approximate future predictions are possible too

Once a candidate passes the overlap guardrails, the package can importance-resample the joint posterior using the candidate weights and temporarily put those resampled draws into a PyMC-Marketing 1.x DataTree.

Future prediction is then delegated to PyMC-Marketing's own predict_posterior implementation.

The example config includes:

prediction:
  data: ../data/future_media_price_plan.csv
  draws: 1000

That means you can ask questions such as:

• does the candidate prior materially change next-quarter demand predictions?
• do contribution shares move enough to alter a budget recommendation?
• does the price response become more economically coherent without rewriting the media story?

Those are still importance-reweighted approximate predictions. Production answers should come from the final validation fit.

Running the case study

This repo targets PyMC-Marketing 1.x and Python 3.12+.

Create an environment:

python -m venv .venv
source .venv/bin/activate
pip install -r requirements.txt
pip install -e ./pymc-marketing-sweep

Render the article immediately

The fast path uses the deterministic posterior fixture:

python scripts/generate_data.py
python scripts/summarise_data.py
PYTHONPATH=pymc-marketing-sweep/src python scripts/generate_story_artifacts.py
python scripts/make_figures.py
python scripts/render_html.py

Then open:

BLOG.html

Run the full PyMC-Marketing workflow

Fit the deliberately difficult centred model:

python scripts/fit_bad_model.py
python scripts/diagnose_fit.py \
  --model outputs/runtime/bad_model.nc \
  --label bad

Fit the non-centred broad-prior reference:

python scripts/fit_reference_model.py
python scripts/diagnose_fit.py \
  --model outputs/runtime/reference_model.nc \
  --label reference

Inspect the saved model before writing or changing a sweep:

pmm-sweep inspect outputs/runtime/reference_model.nc

Run the configured sweep:

pmm-sweep run configs/sku_customer_brand_sweep.yaml

The sweep writes:

outputs/runtime/sweep/
├── sweep_results.csv
├── recommended_scenario.yaml
├── recommended_model_config.py
├── recommended_importance_weights.npz
├── recommended_weighted_posterior_summary.nc
└── recommended_future_prediction_summary.csv

Finally run the recommended real refit:

python scripts/fit_fixed_model.py
python scripts/diagnose_fit.py \
  --model outputs/runtime/fixed_model.nc \
  --label fixed
python scripts/extract_price_elasticities.py

Or run the complete sequence with:

./run_demo.sh

Why this workflow is useful

A complex MMM has two requirements that pull in opposite directions.

It needs enough flexibility to represent the business. Customers differ. Pack sizes differ. Media response differs. Price sensitivity differs. Complete pooling would be statistically convenient and commercially wrong.

But every extra dimension creates another place where weak identification can become difficult posterior geometry.

The answer is not to avoid hierarchical MMMs. It is to make model iteration more disciplined.

A prior sweep turns:

"try a few priors and see"

into a set of explicit questions:

• Is this candidate supported by the posterior we already sampled well?
• How much effective information survives the reweighting?
• Which business quantities move, and by how much?
• Is the more aggressive prior really worth that movement?
• What exact PyMC-Marketing configuration should we validate next?

That changes the workflow from repeated expensive guessing to one reference fit, a cheap screening stage, and one targeted validation fit.

The final lesson is not that NUTS needed more patience.

It is that a model can be statistically sophisticated and still need a better-shaped problem.

Technical notes

This repository targets PyMC-Marketing 1.x, uses the public pymc_marketing.mmm.MMM API, native PyMC-Extras Prior objects, supported PyMC-Marketing special priors and DataTree-backed saved model artefacts. The bundled sweep package is the PyMC-Marketing 1.x-native version and includes its own unit, stress and integration tests.

See SOURCES.md for the official references used to align the example with PyMC-Marketing 1.0 and the associated PyMC/ArviZ stack.