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

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. 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 a logistic saturation function:

If we know that saturation typically occurs 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% saturation at typical spend:
# solve: 0.7 = 1/(1 + exp(-lambda * 1M))
# lambda ≈ 0.85 / 1M

import pymc as pm

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

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.

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 in the framework
from mmm_framework import ModelConfigBuilder, AdstockConfigBuilder, SaturationConfigBuilder

config = (
    ModelConfigBuilder()
    .with_adstock(
        AdstockConfigBuilder()
        .with_type("geometric")
        .with_alpha_prior("Beta", alpha=3, beta=3)  # Centered at 0.5
        .build()
    )
    .with_saturation(
        SaturationConfigBuilder()
        .with_type("logistic")
        .with_lam_prior("Gamma", alpha=2, beta=1/mean_spend)  # Scale to data
        .build()
    )
    .with_media_prior("HalfNormal", sigma=0.5)  # Positive effects
    .build()
)

Validation Strategies

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

Framework Validation Code

# Holdout validation
train_data, test_data = panel.split(holdout_periods=8)
model = BayesianMMM(train_data, config)
model.fit()

# Predict on holdout
predictions = model.predict(test_data)
coverage = model.compute_coverage(test_data, ci=0.9)
print(f"90% CI Coverage: {coverage:.1%}")  # Should be ~90%

# Compare to geo experiment
geo_test_lift = 0.12  # 12% lift from experiment
mmm_lift = model.compute_channel_lift("tv", scenario="actual_vs_zero")
print(f"MMM Lift: {mmm_lift.mean():.1%} [{mmm_lift.hdi(0.9)}]")
print(f"Geo Test Lift: {geo_test_lift:.1%}")
# Check if geo test result within MMM credible interval

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.