Bayesian Workflow

A principled approach to building, checking, and refining Bayesian models. This guide follows the framework established by Gelman et al. (2020) and Gabry et al. (2019), adapted for Marketing Mix Models with practical PyMC implementations.

Introduction

Bayesian data analysis is about more than computing a posterior distribution. It is an iterative process of model building, inference, model checking, and model expansion. Each step informs the next, and visualization is indispensable at every stage.

For Marketing Mix Models, this workflow is especially important because:

• Identification challenges: MMMs often have weak identification, making prior choice critical
• Domain knowledge is rich: We have strong beliefs about plausible effect sizes and shapes
• Stakes are high: Budget allocation decisions depend on getting the model right
• Overfitting is tempting: With many channels and controls, specification shopping is easy

Beneath the mechanics, the seven phases share one purpose: protecting the causal interpretation of the model's attribution. Prior predictive checks encode marketing knowledge before the data can speak; sensitivity analysis asks whether the attribution survives perturbation. A model that passes these gates has earned the right to say "this is what your media caused," not merely "this is what fit best."

The Full Workflow

The workflow is not strictly linear. At each stage, we may discover problems that send us back to an earlier stage. This is normal and expected—the goal is to build understanding, not just to produce a single "correct" answer.

Philosophy

The workflow emphasizes several key principles:

1. Model Building

Model building starts with a generative story: how do we believe the data were generated? For an MMM, this story includes the baseline sales process, the effect of marketing activities, and the noise structure.

Prior Elicitation

Priors encode domain knowledge before seeing the data. Think of them as your media plan's institutional memory—what decades of marketing science says about plausible effect sizes, carryover, and saturation before this particular dataset speaks. In MMM, we often have strong beliefs:

Encoding Domain Knowledge

The prior is the formal mechanism for encoding domain knowledge. Consider the saturation parameter \(\lambda\) in the framework's default saturation function:

This is exponential saturation: \(f(0) = 0\), concave, asymptoting to 1. (The config name SaturationConfig.logistic() is historical — the shipped function is the curve above, not the logistic sigmoid; see transforms/saturation.py. Hill, Michaelis–Menten, and tanh shapes are available per channel via SaturationConfig.) If we know that diminishing returns are well established around 1M USD spend, we can set a prior on \(\lambda\) such that \(f(1M)\) has most of its mass in a plausible range (e.g., 60–90% of maximum effect).

# Prior on saturation parameter
# If typical_spend = 1M and we want ~70% of max effect at typical spend:
# solve: 0.7 = 1 - exp(-lambda * 1M)
# lambda = -ln(0.3) / 1M ≈ 1.20 / 1M

import pymc as pm

with pm.Model():
    # Lognormal prior on lambda, centered on expected value
    lam = pm.Lognormal(
        "lambda",
        mu=np.log(1.20 / typical_spend),
        sigma=0.5  # allows factor of ~3x variation
    )

# Note: the shipped model applies saturation to adstocked spend normalized
# to roughly [0, 1] per channel, so its default prior sat_lam ~ Exponential(0.5)
# lives on that normalized scale, not on raw dollars.

2. Prior Predictive Checks

Before fitting the model to data, we sample from the prior predictive distribution:

This generates fake datasets from the model using only the priors—no observed data involved. The question is: could the observed data plausibly have come from this generative process?

Prior Pushforward

Weakly Informative Priors

A prior is weakly informative if the prior predictive distribution:

• Has mass on plausible datasets (including extreme but possible cases)
• Has no mass on impossible datasets (negative sales, infinite values)
• Is broader than the observed data (not cherry-picked to match)

# Prior predictive check in PyMC
with model:
    prior_pred = pm.sample_prior_predictive(samples=500)

# Visualize: do these look like plausible sales data?
fig, axes = plt.subplots(2, 5, figsize=(15, 6))
for i, ax in enumerate(axes.flat):
    ax.plot(prior_pred.prior_predictive["y_obs"][0, i])
    ax.set_title(f"Draw {i+1}")
plt.suptitle("Prior Predictive Samples")
plt.tight_layout()

3. Computational Fitting

Once satisfied with the prior predictive, we fit the model using MCMC (typically NUTS in PyMC). The goal is to obtain samples from the posterior:

MCMC Diagnostics

MCMC produces samples from the posterior, not the posterior itself. We must verify that the samples are reliable:

Divergences and Pathologies

Solutions to common MCMC problems:

# Standard fitting with diagnostics
with model:
    trace = pm.sample(
        draws=2000,
        tune=2000,
        chains=4,
        target_accept=0.95,
        random_seed=42,
        return_inferencedata=True
    )

# Check diagnostics
print(az.summary(trace, var_names=["beta", "alpha", "lambda"]))
# Look for: R-hat < 1.01, ESS > 400, no divergences

4. Posterior Predictive Checks

After fitting, we generate data from the posterior predictive distribution:

If the model is good, replicated data \(y^{\text{rep}}\) should look like the observed data \(y\). Systematic differences indicate model misspecification.

Calibration

A model is calibrated if its stated uncertainty is accurate. For example, 90% prediction intervals should contain the true value 90% of the time.

Leave-One-Out Cross-Validation

LOO-CV estimates out-of-sample predictive performance without actually refitting the model for each held-out observation. PSIS-LOO uses importance sampling to approximate LOO efficiently:

# Compute LOO
loo = az.loo(trace, pointwise=True)
print(loo)

# Check Pareto k diagnostics
# k > 0.7: that observation is highly influential, LOO approximation unreliable
high_k = np.where(loo.pareto_k > 0.7)[0]
if len(high_k) > 0:
    print(f"Warning: {len(high_k)} observations with k > 0.7")

5. Model Comparison

When we have multiple candidate models, we compare them using their expected predictive performance. The key quantity is the expected log pointwise predictive density (ELPD):

ELPD and LOO-CV

Higher ELPD is better—it means the model makes better predictions on held-out data. The standard error of the ELPD difference tells us whether the difference is meaningful.

# Compare multiple models
comparison = az.compare({
    "standard": trace_standard,
    "nested": trace_nested,
    "hierarchical": trace_hierarchical
})
print(comparison)

# Visualize comparison
az.plot_compare(comparison)

Model Topology

Gelman et al. (2020) recommend thinking of models as a topology—a network of related models that can be expanded or simplified. Each model comparison teaches us something about the structure of the problem. Below is the full model topology for the MMM Framework:

Model Selection Guide

6. Sensitivity Analysis

How much do our conclusions depend on our modeling choices? Sensitivity analysis systematically varies assumptions to assess robustness. For an MMM, this is the causal stress test: if a channel's attribution survives reasonable perturbations of priors, trend, and influential weeks, you can treat it as a finding rather than an artifact of one specification.

Prior Sensitivity

Refit the model with different priors and compare posteriors. If conclusions change dramatically with reasonable prior variations, the data are not strongly informative about those parameters.

Data Influence

The Pareto k diagnostics from LOO-CV identify influential observations. High-influence observations warrant investigation:

• Are they outliers? Consider robust likelihoods (Student-t)
• Are they unusual periods? (Promotions, COVID, competitive activity)
• Are they data errors? Check the raw data

7. Model Expansion

Based on posterior predictive checks and domain understanding, we iteratively expand the model to address deficiencies. The MMM Framework provides structured expansion paths:

Expansion Decision Tree

When to Stop

Model expansion has diminishing returns. Stop when:

• Posterior predictive checks show no systematic problems
• Additional complexity doesn't improve ELPD beyond standard error
• The model answers the scientific questions with adequate precision
• You've started "shopping" for specifications that give desired answers

MMM-Specific Workflow

Applying the Bayesian workflow to Marketing Mix Models involves specific considerations for the MMM Framework's model classes:

Workflow by Model Type

MMM Prior Guidance

For each major parameter class in the framework, here is guidance on setting informative priors:

# Example: configure priors per channel with the shipped config objects
from mmm_framework.config import (
    AdstockConfig, MediaChannelConfig, PriorConfig, SaturationConfig,
)

tv = MediaChannelConfig(
    name="TV",
    adstock=AdstockConfig.geometric(          # parametric path (the default;
        l_max=8,                              # ModelConfig.use_parametric_adstock=True)
        alpha_prior=PriorConfig.beta(alpha=2, beta=2),  # longer-memory belief for TV
    ),
    saturation=SaturationConfig.logistic(     # exponential-saturation form 1 - e^(-λx)
        lam_prior=PriorConfig.gamma(alpha=2.0, beta=1.0),
    ),
)

Validation Strategies

MMM validation goes beyond statistical checks. The framework supports multiple validation approaches:

Framework Validation Code

# Holdout validation: rolling-origin backtest (refits on expanding
# training windows, grades genuine out-of-time forecasts)
from mmm_framework.validation import BacktestConfig, run_backtest

result = run_backtest(
    BayesianMMM(panel, config),
    BacktestConfig(min_train_size=104, horizon=13, step=13),
)
print(result.summary())    # MAPE/MASE vs naive + seasonal-naive baselines
print(result.coverage())   # nominal vs empirical interval coverage (~90% at 0.9)

# Compare to a geo experiment: the model's implied lift for the channel
contrib = model.compute_counterfactual_contributions(
    compute_uncertainty=True, hdi_prob=0.90
)
mmm_lift = contrib.total_contributions["TV"] / baseline_sales
print(f"MMM lift: {mmm_lift:.1%} "
      f"[{contrib.contribution_hdi_low['TV'] / baseline_sales:.1%}, "
      f"{contrib.contribution_hdi_high['TV'] / baseline_sales:.1%}]")
print(f"Geo test lift: {0.12:.1%}")  # from the experiment readout
# If the geo result falls outside the MMM interval, calibrate:
# feed the experiment in as a likelihood term via add_experiment_calibration()
# or compare systematically with ModelValidator + ValidationConfigBuilder
# .with_calibration([LiftTestResult(...)]).

References

• Gelman, A., Vehtari, A., Simpson, D., et al. (2020). Bayesian Workflow. arXiv:2011.01808.
• Gabry, J., Simpson, D., Vehtari, A., Betancourt, M., & Gelman, A. (2019). Visualization in Bayesian workflow. J. R. Stat. Soc. A, 182(2), 389-402.
• Vehtari, A., Gelman, A., & Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing, 27(5), 1413-1432.
• Betancourt, M. (2017). A Conceptual Introduction to Hamiltonian Monte Carlo. arXiv:1701.02434.