Mmm Modeling

Media Mix Modeling (MMM) Skill

Bayesian Media Mix Modeling using PyMC-Marketing (primary) with a custom scipy-based fallback for restricted environments. Includes Meridian-inspired scenario planning patterns for forward-looking budget optimization.

Architecture

Detection logic — try import pymc_marketing first; if ImportError, fall back to custom scripts.

Agency Growth OS Integration

This skill operates within the 17-skill Agency Growth OS. Its position in the three cycles:

Execution Cycle

media-routing-planner → mmm-modeling → measurement-incrementality
                                     ↘ budget_optimization → media-routing-planner (re-optimize)

Intelligence Cycle

weekly-control-tower → performance-diagnosis → mmm-modeling (if root cause is channel mix)
mmm-modeling → qbr-generator (quarterly decomposition + scenario slides)
mmm-modeling → client-memory-synthesizer (store model parameters + integrity score)

Chain Routing

After mmm-modeling completes, suggest the next skill via sendPrompt():

Data Ingestion

Digital media data → Adspirer MCP (automated)

When the model needs digital campaign spend data (Google Ads, Meta, LinkedIn, TikTok), call the Adspirer MCP tools already connected:

# Claude calls Adspirer tools to pull spend data:
# - get_campaign_performance → spend by channel × date
# - get_campaign_structure → channel taxonomy
# - analyze_search_terms → search query volume (control variable)

Use `tool_search` for Adspirer tools when the user mentions digital spend data. The skill does NOT embed Adspirer tool calls directly — Claude resolves them at runtime via the connected MCP.

Non-digital data → Human batch upload (CSV/Excel)

Data the human must provide (no MCP available):

Workflow: Human uploads → data_validator.py validates → context layer resolves priors → model fits.

Visualization Strategy

This skill ALWAYS produces visual deliverables. The strategy controls WHERE they render to avoid crashing the app.

In-chat: Tables + lightweight SVG (always)

Every MMM run produces formatted tables in chat as immediate output:

• Channel ROAS with HDI ranges (table)
• Integrity scorecard pass/warn/fail (table)
• Contribution share % breakdown (table)
• Scenario comparison (table)

For visual charts in chat, use the Visualizer with lightweight SVG only:

• Max 6 channels × 50 data points per chart
• No animations, no JS interactivity
• Keep total SVG under 50KB
• One chart per Visualizer call (never stack multiple)

Suitable for inline SVG:

• Horizontal bar: contribution waterfall, ROAS comparison
• Line chart: saturation curves (≤6 lines), efficient frontier
• Donut: four-way decomposition (base/media/trade/interaction)
• Status indicators: integrity scorecard with color-coded pass/warn/fail

File deliverables: Full charts + dashboards (always generated)

After in-chat summary, ALWAYS generate file artifacts in /mnt/user-data/outputs/:

Chart generation in HTML files

HTML dashboards use inline <svg> or lightweight Chart.js (CDN) — NOT matplotlib, NOT plotly, NOT heavy React state. This ensures they open fast in browser without crashing.

<!-- Pattern for HTML dashboard charts -->
<script src="https://cdn.jsdelivr.net/npm/chart.js@4"></script>
<canvas id="saturation_chart"></canvas>
<script>
  new Chart(document.getElementById('saturation_chart'), {
    type: 'line',
    data: { /* from mmm results JSON */ },
    options: { responsive: true, animation: false }
  });
</script>

PPTX charts

Use the pptx skill to embed static chart images in slides. Generate chart as PNG via matplotlib in script → embed in slide. Charts render server-side, no app crash risk.

Quick Start (PyMC-Marketing)

import arviz as az
import numpy as np
import pandas as pd
from pymc_extras.prior import Prior
from pymc_marketing.mmm import GeometricAdstock, LogisticSaturation
from pymc_marketing.mmm.multidimensional import MMM

data_df = pd.read_csv("data.csv", parse_dates=["date"])
X = data_df.drop(columns=["y"])
y = data_df["y"]

mmm = MMM(
    date_column="date",
    channel_columns=["tv", "radio", "social"],
    target_column="y",
    adstock=GeometricAdstock(l_max=6),
    saturation=LogisticSaturation(),
    yearly_seasonality=5,
)

# CRITICAL: Always fit on FULL dataset. Train/test splits are ONLY for stability assessment.
mmm.build_model(X, y)
mmm.fit(X=X, y=y, nuts_sampler="nutpie", target_accept=0.9, random_seed=42)
mmm.sample_posterior_predictive(X=X, random_seed=42)

Quick Start (Custom Fallback)

# Validate data
python scripts/data_validator.py data.csv

# Generate report from results JSON
python scripts/report_generator.py model.json --roi roi.json --html

Reference Architecture

Scripts

Model Specification (PyMC-Marketing)

Transformations

GeometricAdstock — exponential decay carryover:

from pymc_marketing.mmm import GeometricAdstock
adstock = GeometricAdstock(l_max=6, normalize=True)
# Default prior: alpha ~ Beta(1, 3) — favors fast decay
# Custom: priors={"alpha": Prior("Beta", alpha=2, beta=5, dims="channel")}

LogisticSaturation — S-shaped diminishing returns:

from pymc_marketing.mmm import LogisticSaturation
saturation = LogisticSaturation()
# Default: lam ~ Gamma(3, 1), beta ~ HalfNormal(2)
# Informed: priors={"beta": Prior("HalfNormal", sigma=spend_shares, dims="channel")}

Full Model Config

from pymc_extras.prior import Prior

model_config = {
    "intercept": Prior("Normal", mu=0.2, sigma=0.05),
    "saturation_beta": Prior("HalfNormal", sigma=spend_shares, dims="channel"),
    "gamma_control": Prior("Normal", mu=0, sigma=1, dims="control"),
    "gamma_fourier": Prior("Laplace", mu=0, b=1, dims="fourier_mode"),
    "likelihood": Prior("TruncatedNormal", lower=0, sigma=Prior("HalfNormal", sigma=1)),
}

See references/model_specification.md for constructor reference, all saturation alternatives (Hill, MichaelisMenten, Tanh, Root, etc.), hierarchical prior patterns, and scaling config.

Multidimensional (Geo-Hierarchical)

Activate with dims=("geo",). Use partial pooling (default recommendation):

from pymc_marketing.special_priors import LogNormalPrior

model_config = {
    "saturation_beta": LogNormalPrior(
        mean=Prior("Gamma", mu=1.0, sigma=1.0),
        std=Prior("HalfNormal", sigma=1.0),
        dims=("channel", "geo"), centered=False,
    ),
}

Workflow

Client Brief (industry, category, market, channels, KPI)
    ↓
Context Resolution (context_resolver.py → calibrated priors)
    ↓
EDA & Data Prep (data_validator.py)
    ↓
Acid-Test Pre-Model (acid_test.py → integrity baseline)
    ↓
Model Specification (priors from context layer)
    ↓
Build Model (mmm.build_model)
    ↓
Prior Predictive Checks
    ↓
[Optional] Add Lift Test Calibration
    ↓
[Optional] Add Trade Marketing Variables (see trade_marketing_decomposition.md)
    ↓
Fit on FULL Dataset (mmm.fit with nutpie)
    ↓
Diagnostics (divergences=0, R-hat<1.01, ESS>400)
    ↓
Acid-Test Post-Model (absorption check, ROI plausibility, permutation)
    ↓
Media Deep Dive (contributions, ROAS, saturation, sensitivity)
    ↓
Four-Way Decomposition (base + media + trade + interaction)
    ↓
Budget Optimization + Scenario Planning
    ↓
Report Generation (HTML dashboard + PPTX deck + JSON API)

Key APIs

Incrementality (preferred for ROAS/CAC)

roas = mmm.incrementality.contribution_over_spend(frequency="all_time")
marginal_roas = mmm.incrementality.marginal_contribution_over_spend(frequency="all_time", spend_increase_pct=0.01)
cac = mmm.incrementality.spend_over_contribution(frequency="quarterly")

Summary DataFrames

mmm.summary.posterior_predictive()    # mean, median, HDI, observed
mmm.summary.contributions()           # per-channel contributions
mmm.summary.roas()                     # ROAS with HDI
mmm.summary.saturation_curves()        # saturation response
mmm.summary.adstock_curves()           # decay profiles

Budget Optimization

from pymc_marketing.mmm.multidimensional import MultiDimensionalBudgetOptimizerWrapper

optimizer = MultiDimensionalBudgetOptimizerWrapper(model=mmm, start_date=..., end_date=...)
allocation, result = optimizer.optimize_budget(budget=1_000_000, budget_bounds=bounds)
response = optimizer.sample_response_distribution(allocation_strategy=allocation, include_carryover=True)

See references/budget_optimization.md for bounds setup, channel fixing, custom constraints, and budget sweeps.

Scenario Planning

See references/scenario_planning.md for Meridian-inspired patterns adapted for PyMC-Marketing:

• Fixed budget optimization (maximize ROI at given budget)
• Flexible budget optimization (find max budget at target ROI)
• Budget sweeps with efficient frontier
• Flighting / temporal distribution
• Multi-scenario comparison (conservative/moderate/aggressive)
• Cost-per-media-unit sensitivity

Save/Load

mmm.save("mmm_model.nc", engine="h5netcdf")
loaded = MMM.load("mmm_model.nc")

YAML Specification

from pymc_marketing.mmm.builders.yaml import build_mmm_from_yaml
mmm = build_mmm_from_yaml("model_spec.yaml", X=X, y=y)

Layer 1: Acid-Test Validation

Runs pre/post-model integrity checks to verify MMM truthfulness. This is how you expose holdco manipulation.

from scripts.acid_test import AcidTestValidator

# Pre-model (before fitting)
validator = AcidTestValidator(df, "date", "sales", ["tv", "social", "search"])
pre_report = validator.run_pre_model_tests()

# Post-model (after fitting, with MMM results)
post_report = validator.run_post_model_tests({
    "r2": 0.82,
    "channel_contributions": {"tv": 500000, "social": 200000, "search": 300000},
    "channel_roas": {"tv": 2.1, "social": 1.5, "search": 3.2},
    "baseline_pct": 0.55,
})

print(validator.to_summary())
validator.to_json("integrity_report.json")

Key tests: Holt-Winters baseline, Granger causality, VIF, baseline absorption check, ROI plausibility. See references/acid_test_validation.md.

Layer 2: Context Layer

Resolves industry/category/market benchmarks into calibrated priors:

from scripts.context_resolver import ContextResolver

resolver = ContextResolver(
    industry="CPG", category="Beverages", market="Mexico",
    channels=["tv", "social", "search", "ooh"],
    distribution_model="indirect", trade_marketing_share=0.45,
)
brief = resolver.to_json("context_brief.json")
print(resolver.to_summary())

# Use resolved priors in model config
# brief.suggested_model_config → ready-to-use Prior strings
# brief.roi_benchmarks → feeds acid-test ROI plausibility check
# brief.adstock_priors → calibrated Beta distributions per channel

See references/context_layer.md for full benchmark databases and adaptation rules.

Layer 3: Trade Marketing Decomposition

Four-way decomposition: base demand + media-driven + trade-driven + interaction.

Add trade variables as controls with domain-informed priors:

control_columns = ["promo_depth", "price_index", "distribution_acv", "feature_flag"]

model_config = {
    "saturation_beta": Prior("HalfNormal", sigma=spend_shares, dims="channel"),
    "gamma_control": Prior("Normal", mu=trade_prior_means, sigma=trade_prior_sds, dims="control"),
}

For sophisticated decomposition with custom trade response curves, use a custom PyMC model. See references/trade_marketing_decomposition.md for flexible data schema (Tier 1-4), functional forms, and the holdco exposure play.

Critical Rules

1. Always fit on FULL dataset — train/test splits are ONLY for TimeSliceCrossValidator stability checks 2. Zero divergences required — any divergences invalidate the posterior 3. R-hat < 1.01 for all parameters before proceeding 4. Call `add_original_scale_contribution_variable` before sample_posterior_predictive to get *_original_scale variables 5. Use `mmm.incrementality` for ROAS (accounts for adstock carryover), not element-wise division 6. Use `MultiDimensionalBudgetOptimizerWrapper` not BudgetOptimizer directly (handles geo allocation) 7. Never present platform ROAS as incremental without caveating (Agency Growth OS rule) 8. Search web for current industry benchmarks before comparing ROI results

Diagnostics Quick Reference

See references/diagnostics_benchmarks.md for full checklist and industry ROI ranges.

Acid-Test Validation Layer (Model Integrity Report)

The acid-test layer runs before and after the MMM to answer one question: "Is the MMM telling the truth, or is it fitting noise to media variables?"

This is the weapon against holdco MMMs that overattribute to channels where the agency earns margin.

Architecture

Raw Data
  ↓
┌─────────────────────────────────────────────────┐
│ ACID-TEST BATTERY (pre-model)                   │
│                                                  │
│  1. Baseline Forecast (Holt-Winters / ETS)      │
│  2. Granger Causality (per channel)             │
│  3. VIF Multicollinearity                        │
│  4. Stationarity (ADF)                           │
│  5. Spend Concentration Index                    │
│                                                  │
│  → Produces: pre_model_integrity.json            │
└─────────────────────────────────────────────────┘
  ↓
MMM Fit (PyMC-Marketing or fallback)
  ↓
┌─────────────────────────────────────────────────┐
│ ACID-TEST BATTERY (post-model)                  │
│                                                  │
│  6. Permutation Test (shuffle media, refit)     │
│  7. Baseline Absorption Check                    │
│  8. Leave-One-Channel-Out (LOCO)                │
│  9. ROI Plausibility vs Industry Benchmarks     │
│ 10. Contribution Stability (bootstrap)          │
│                                                  │
│  → Produces: post_model_integrity.json           │
└─────────────────────────────────────────────────┘
  ↓
Model Integrity Report (HTML + PPTX + JSON)

Test Descriptions

Pre-Model Tests

1. Baseline Forecast (Holt-Winters / ETS)

Fits an exponential smoothing model using ONLY the target variable (no media). If this model explains >85% of variance, the MMM's marginal value is questionable.

from statsmodels.tsa.holtwinters import ExponentialSmoothing

ets_model = ExponentialSmoothing(
    y, trend="add", seasonal="add", seasonal_periods=52
).fit(optimized=True)
ets_pred = ets_model.fittedvalues
ets_r2 = 1 - np.sum((y - ets_pred)**2) / np.sum((y - y.mean())**2)

Key metric: media_marginal_r2 = mmm_r2 - ets_r2. If < 0.05, the MMM adds almost nothing beyond seasonality.

2. Granger Causality

Tests whether media spend temporally precedes KPI movement, or just correlates.

from statsmodels.tsa.stattools import grangercausalitytests

for channel in media_cols:
    test_data = pd.DataFrame({"kpi": y, "media": X[channel]})
    result = grangercausalitytests(test_data, maxlag=4, verbose=False)
    # p-value < 0.05 → channel Granger-causes KPI

Channels that fail Granger should be flagged — their MMM coefficients may be spurious.

3. VIF Multicollinearity

from statsmodels.stats.outliers_influence import variance_inflation_factor

vif_data = X[media_cols]
vif = pd.DataFrame({
    "channel": media_cols,
    "VIF": [variance_inflation_factor(vif_data.values, i) for i in range(len(media_cols))]
})

Holdco red flag: If TV and digital have VIF > 10, the agency can shift attribution between them at will.

4. Stationarity (Augmented Dickey-Fuller)

from statsmodels.tsa.stattools import adfuller

adf_stat, adf_p, _, _, _, _ = adfuller(y, maxlag=12)
# p < 0.05 → stationary (good for MMM)
# p ≥ 0.05 → unit root, needs differencing or detrending

5. Spend Concentration Index

Measures how evenly spend is distributed across time. Highly concentrated (pulsed) spend makes identification harder.

def spend_concentration(x):
    """Herfindahl-like index. 1/n = perfectly even, 1.0 = all in one period."""
    shares = x / x.sum()
    return (shares ** 2).sum()

for channel in media_cols:
    hhi = spend_concentration(X[channel])
    # hhi > 0.1 → concentrated spend, harder to identify

Post-Model Tests

6. Permutation Test

Shuffle each media variable independently, refit, measure R² drop. If R² doesn't drop, the channel's coefficient is spurious.

def permutation_test(mmm, X, y, channel, n_permutations=100):
    """Shuffle one channel, refit, compare R²."""
    original_r2 = compute_r2(mmm, X, y)
    perm_r2s = []
    for _ in range(n_permutations):
        X_perm = X.copy()
        X_perm[channel] = np.random.permutation(X_perm[channel].values)
        # Refit or use posterior predictive with permuted data
        perm_r2 = compute_r2_with_permuted(mmm, X_perm, y)
        perm_r2s.append(perm_r2)
    r2_drop = original_r2 - np.mean(perm_r2s)
    p_value = np.mean([pr2 >= original_r2 for pr2 in perm_r2s])
    return r2_drop, p_value

7. Baseline Absorption Check

Compares baseline % between the MMM and the Holt-Winters forecast.

mmm_baseline_pct = intercept_contribution / total_predicted
ets_baseline_pct = ets_pred.mean() / y.mean()

absorption_gap = ets_baseline_pct - mmm_baseline_pct
# If gap > 0.15 → MMM is absorbing baseline demand into media channels

Holdco red flag: If the MMM shows 35% baseline but Holt-Winters shows 70%, the MMM is attributing 35% of natural demand to media. This inflates media ROI and justifies the agency's fee.

8. Leave-One-Channel-Out (LOCO)

Remove each channel, refit, measure how much total predicted KPI drops.

for channel in media_cols:
    X_loco = X.drop(columns=[channel])
    # Refit MMM without this channel
    # Compare: if total R² barely changes, channel is redundant

9. ROI Plausibility

Compare MMM-derived ROI against industry benchmarks (loaded from context layer).

for channel, roi in mmm_roas.items():
    benchmark_range = industry_benchmarks.get(channel, (0.5, 5.0))
    if roi > benchmark_range[1] * 2:
        flag = "IMPLAUSIBLE_HIGH"
    elif roi < benchmark_range[0] * 0.5:
        flag = "IMPLAUSIBLE_LOW"

10. Contribution Stability (Bootstrap)

Resample data with replacement, refit N times, measure coefficient of variation of channel contributions.

# If CV of a channel's contribution share > 0.50 across bootstraps,
# that channel's attribution is unstable and should not drive budget decisions.

Integrity Score

Each test produces a pass/warn/fail. The composite score:

integrity_score = {
    "baseline_independence": ets_r2 < 0.85,          # weight: 20%
    "causal_direction": granger_pass_rate > 0.60,     # weight: 15%
    "multicollinearity": max_vif < 10,                # weight: 15%
    "permutation_validity": min_r2_drop > 0.01,       # weight: 20%
    "baseline_absorption": absorption_gap < 0.15,     # weight: 15%
    "roi_plausibility": all_within_benchmarks,        # weight: 15%
}

total_score = sum(weight * passed for (_, passed), weight
                  in zip(integrity_score.items(), weights))

Output Format

The acid-test produces three artifacts: 1. integrity_report.json — machine-readable results for API consumers 2. integrity_report.html — interactive dashboard with test-by-test results 3. integrity_slide.pptx — single executive slide for C-suite decks

--- See also: context_layer.md, diagnostics_benchmarks.md

Budget Optimization

Setup

from pymc_marketing.mmm.multidimensional import MultiDimensionalBudgetOptimizerWrapper

last_date = pd.Timestamp(X[date_column].max())
start_date = last_date + pd.Timedelta(weeks=1)
end_date = start_date + pd.Timedelta(weeks=12)

optimizer = MultiDimensionalBudgetOptimizerWrapper(
    model=mmm, start_date=str(start_date), end_date=str(end_date),
)

Budget Bounds

import xarray as xr

# Single-geo
budget_bounds = xr.DataArray(
    data=np.array([[0.5, 1.5], [0.3, 2.0], [0.5, 1.5]]) * equal_share,
    dims=["channel", "bound"],
    coords={"channel": channel_columns, "bound": ["lower", "upper"]},
)

# Multi-geo
budget_bounds = xr.DataArray(
    data=np.stack([np.full((n_ch, n_geos), 0.0),
                   np.full((n_ch, n_geos), max_budget)], axis=-1),
    dims=["channel", "geo", "bound"],
    coords={"channel": channel_columns, "geo": geos, "bound": ["lower", "upper"]},
)

Running Optimization

allocation, result = optimizer.optimize_budget(
    budget=budget_per_period,
    budget_bounds=budget_bounds,
    minimize_kwargs={"method": "SLSQP", "options": {"ftol": 1e-4, "maxiter": 10_000}},
)

Response Sampling

response = optimizer.sample_response_distribution(
    allocation_strategy=allocation,
    additional_var_names=["channel_contribution_original_scale"],
    include_last_observations=True, include_carryover=True, noise_level=0.05,
)

Visualization

optimizer.plot.budget_allocation(samples=response)
mmm.plot.allocated_contribution_by_channel_over_time(response)

Advanced Patterns

Fix Certain Channels

budgets_to_optimize = xr.DataArray(
    data=[True, True, False], dims=["channel"],
    coords={"channel": channel_columns},
)
allocation, result = optimizer.optimize_budget(
    budget=budget_per_period, budget_bounds=budget_bounds,
    budgets_to_optimize=budgets_to_optimize, ...
)

Custom Temporal Distribution (Flighting)

time_weights = np.linspace(1.5, 0.5, n_periods)
time_weights /= time_weights.sum()

budget_distribution = xr.DataArray(
    data=np.tile(time_weights[:, None], (1, len(channel_columns))),
    dims=["date", "channel"], coords={"channel": channel_columns},
)

allocation, result = optimizer.optimize_budget(
    budget=budget_per_period, budget_bounds=budget_bounds,
    budget_distribution_over_period=budget_distribution, ...
)

Custom Constraints

from pymc_marketing.mmm.constraints import Constraint

def tv_ge_2x_radio(budgets_sym, total_budget_sym, optimizer):
    tv_idx = list(channel_columns).index("tv")
    radio_idx = list(channel_columns).index("radio")
    return budgets_sym[tv_idx] - 2 * budgets_sym[radio_idx]

constraint = Constraint(key="tv_ge_2x_radio", constraint_type="ineq",
                        constraint_fun=tv_ge_2x_radio)

allocation, result = optimizer.optimize_budget(
    budget=budget_per_period, budget_bounds=budget_bounds,
    constraints=[constraint], ...
)

Budget Sweep (Efficient Frontier)

budget_levels = np.linspace(50_000, 500_000, 10)
sweep_results = []
for budget in budget_levels:
    alloc, res = optimizer.optimize_budget(budget=budget, budget_bounds=bounds, ...)
    sweep_results.append({"budget": budget, "allocation": alloc, "result": res})

In-Sample Uplift Estimation

X_cf = X.copy()
X_cf["tv"] *= 1.2
X_cf["radio"] *= 0.8
response_cf = mmm.sample_posterior_predictive(X_cf, extend_idata=False, random_seed=42)
uplift = (response_cf["posterior_predictive"]["y"].mean(dim=("chain", "draw"))
          - mmm.idata["posterior_predictive"]["y"].mean(dim=("chain", "draw")))

--- See also: scenario_planning.md, media_deep_dive.md

Context Layer (Industry-Calibrated Priors)

The context layer resolves industry, category, market, and SKU-level benchmarks before model specification, so the MMM encodes domain knowledge rather than relying on flat priors.

This is how you break the holdco "proprietary benchmarks" moat — all the data is publicly available; the value is in the systematic resolution and encoding.

Architecture

Client Brief (industry, category, market, channels)
  ↓
┌─────────────────────────────────────────────────────┐
│ CONTEXT RESOLVER                                     │
│                                                      │
│  1. Industry Classification (GICS/NAICS)            │
│  2. Web Search: benchmarks for {industry} × {market}│
│  3. Channel Benchmark Resolution                     │
│  4. Prior Calibration (map benchmarks → priors)     │
│  5. Baseline Expectation Setting                     │
│                                                      │
│  → Produces: context_brief.json                      │
│  → Produces: calibrated_priors.py (model_config)    │
└─────────────────────────────────────────────────────┘
  ↓
MMM Specification (priors informed by context)

Context Schema

The context layer expects a structured brief:

context = {
    "client": "Acme Corp",
    "industry": "CPG",                    # CPG, Automotive, Financial, Retail, Pharma, Tech, Telco, QSR
    "category": "Beverages",              # Category within industry
    "subcategory": "Carbonated Soft Drinks",  # Optional
    "market": "Mexico",                   # Country or region
    "market_tier": "emerging",            # developed, emerging
    "channels": ["tv", "digital_display", "social", "search", "ooh"],
    "kpi": "sales_volume",               # revenue, sales_volume, leads, app_installs, store_visits
    "currency": "MXN",
    "annual_media_budget": 50_000_000,
    "distribution_model": "indirect",     # direct, indirect, hybrid
    "trade_marketing_share": 0.45,        # % of total marketing in trade
    "seasonality_pattern": "holiday_heavy", # flat, holiday_heavy, back_to_school, weather_driven
    "competitive_intensity": "high",      # low, medium, high
    "brand_maturity": "established",      # launch, growth, established, decline
}

Benchmark Sources (Web-Searchable)

The LLM should search for current data from these sources:

Prior Calibration Rules

Adstock Priors by Channel × Industry

ADSTOCK_PRIORS = {
    # (industry, channel) → (alpha_mean, alpha_concentration)
    # Higher alpha = longer carryover
    ("CPG", "tv"):          {"alpha": Prior("Beta", alpha=4, beta=3)},      # ~0.57, long
    ("CPG", "digital"):     {"alpha": Prior("Beta", alpha=2, beta=5)},      # ~0.29, short
    ("CPG", "social"):      {"alpha": Prior("Beta", alpha=2, beta=4)},      # ~0.33
    ("CPG", "search"):      {"alpha": Prior("Beta", alpha=1, beta=6)},      # ~0.14, very short
    ("CPG", "ooh"):         {"alpha": Prior("Beta", alpha=3, beta=3)},      # ~0.50
    ("Automotive", "tv"):   {"alpha": Prior("Beta", alpha=5, beta=2)},      # ~0.71, very long
    ("Automotive", "digital"): {"alpha": Prior("Beta", alpha=3, beta=4)},   # ~0.43
    ("Financial", "tv"):    {"alpha": Prior("Beta", alpha=4, beta=3)},      # ~0.57
    ("Financial", "search"):{"alpha": Prior("Beta", alpha=2, beta=5)},      # ~0.29
    ("Retail", "digital"):  {"alpha": Prior("Beta", alpha=2, beta=5)},      # ~0.29
    ("Pharma", "tv"):       {"alpha": Prior("Beta", alpha=5, beta=2)},      # ~0.71, very long
    ("Tech", "digital"):    {"alpha": Prior("Beta", alpha=2, beta=4)},      # ~0.33
    ("QSR", "tv"):          {"alpha": Prior("Beta", alpha=3, beta=3)},      # ~0.50
    ("Telco", "digital"):   {"alpha": Prior("Beta", alpha=2, beta=4)},      # ~0.33
}

ROI Benchmark Ranges by Industry × Channel

ROI_BENCHMARKS = {
    ("CPG", "tv"):          (1.0, 3.5),
    ("CPG", "digital"):     (1.5, 5.0),
    ("CPG", "social"):      (0.8, 3.0),
    ("CPG", "search"):      (2.0, 6.0),
    ("CPG", "ooh"):         (0.5, 2.0),
    ("Automotive", "tv"):   (0.5, 2.0),
    ("Automotive", "digital"): (1.0, 4.0),
    ("Automotive", "search"):  (2.0, 8.0),
    ("Financial", "tv"):    (0.8, 2.5),
    ("Financial", "digital"): (1.5, 5.0),
    ("Retail", "digital"):  (2.0, 8.0),
    ("Retail", "search"):   (3.0, 10.0),
    ("Pharma", "tv"):       (0.3, 1.5),
    ("Tech", "digital"):    (1.5, 6.0),
    ("QSR", "tv"):          (1.5, 4.0),
    ("QSR", "digital"):     (2.0, 7.0),
    ("Telco", "digital"):   (1.0, 4.0),
}

Baseline Expectations by Industry

BASELINE_EXPECTATIONS = {
    "CPG":        {"range": (0.55, 0.80), "reason": "Strong brand equity, repeat purchase"},
    "Automotive": {"range": (0.30, 0.55), "reason": "Considered purchase, media-influenced"},
    "Financial":  {"range": (0.50, 0.75), "reason": "Regulated, trust-driven"},
    "Retail":     {"range": (0.40, 0.65), "reason": "Mix of brand and promotion-driven"},
    "Pharma":     {"range": (0.60, 0.85), "reason": "Prescription-driven, long cycles"},
    "Tech":       {"range": (0.25, 0.50), "reason": "Performance-driven, short cycles"},
    "QSR":        {"range": (0.45, 0.65), "reason": "Habitual + promotional"},
    "Telco":      {"range": (0.50, 0.70), "reason": "Subscription-based, churn-driven"},
}

Saturation Speed by Category

SATURATION_SPEED = {
    # How quickly channels saturate (lam prior for LogisticSaturation)
    "FMCG_food":        Prior("Gamma", alpha=4, beta=1),    # Fast saturation
    "FMCG_beverages":   Prior("Gamma", alpha=3.5, beta=1),
    "FMCG_personal":    Prior("Gamma", alpha=3, beta=1),
    "Automotive":       Prior("Gamma", alpha=2, beta=1),     # Slow saturation
    "Luxury":           Prior("Gamma", alpha=1.5, beta=1),   # Very slow
    "Financial":        Prior("Gamma", alpha=2.5, beta=1),
    "Tech_SaaS":        Prior("Gamma", alpha=3, beta=1),
    "Retail_ecommerce": Prior("Gamma", alpha=3.5, beta=1),  # Fast
    "QSR":              Prior("Gamma", alpha=3, beta=1),
}

Context Brief Output

The resolver produces a structured brief:

{
  "context_id": "acme_beverages_mexico_2026",
  "industry": "CPG",
  "category": "Beverages",
  "market": "Mexico",
  "resolved_benchmarks": {
    "tv_roi_range": [1.0, 3.5],
    "social_roi_range": [0.8, 3.0],
    "baseline_expected": [0.55, 0.80],
    "adstock_tv_alpha_mean": 0.57
  },
  "calibrated_model_config": {
    "intercept": "Prior('Normal', mu=0.65, sigma=0.10)",
    "saturation_beta": "Prior('HalfNormal', sigma=spend_shares, dims='channel')"
  },
  "benchmark_sources": [
    {"source": "Kantar", "metric": "CPG TV ROI", "value": "2.1x", "year": 2025},
    {"source": "Nielsen", "metric": "Beverages baseline", "value": "68%", "year": 2025}
  ],
  "flags": [
    "Mexico is an emerging market — adjust for lower media saturation thresholds",
    "Indirect distribution — trade marketing decomposition recommended"
  ]
}

How the LLM Uses Context

1. Before model specification: Read context_brief.json → set informative priors 2. During diagnostics: Compare MMM ROI against resolved benchmarks 3. In acid-test: Use baseline expectations for absorption check 4. In reporting: Include benchmark comparisons in executive output 5. In scenario planning: Use industry saturation curves for budget sweep interpretation

Adaptive Resolution

The context layer adapts depth based on available information:

--- See also: acid_test_validation.md, model_specification.md, diagnostics_benchmarks.md

Custom Models

When the MMM class cannot express your model structure (non-standard hierarchies, spline baselines, custom likelihoods), use PyMC-Marketing components with plain PyMC.

Standalone Component Usage

import pymc as pm
from pymc_marketing.mmm import GeometricAdstock, LogisticSaturation

adstock = GeometricAdstock(l_max=6)
saturation = LogisticSaturation()

with pm.Model(coords=coords) as custom_mmm:
    channel_data_ = pm.Data("channel_data", channel_scaled, dims=("date", "geo", "channel"))
    adstocked = adstock.apply(channel_data_, dims=("geo", "channel"))
    channel_contribution = saturation.apply(adstocked, dims=("geo", "channel"))
    # ... add intercept, controls, seasonality, likelihood

Trade-offs

When to Go Custom

• Non-standard hierarchical structures (e.g., brand × channel × geo)
• Spline-based intercepts instead of Fourier seasonality
• Custom likelihoods (Student-t, NegativeBinomial)
• Time-series specific components (AR terms, changepoints)
• Integration with non-media models (demand forecasting, pricing)

--- See also: model_specification.md

Data Analysis & Preparation

Data Format Requirements

National Model (single geo)

Multidimensional (geo-level) — Long Format

Critical: Geo data must be in long format (one row per date × geo).

Minimum Data Requirements

EDA Checklist

1. Temporal structure: Weekly frequency, no gaps, consistent geo coverage 2. Target distribution: Check skewness → apply log/Box-Cox if highly skewed 3. Spend patterns: Continuous vs. pulsed/flighted → informs adstock_first 4. Correlation matrix: High multi-collinearity (VIF > 5) → combine or drop channels 5. Seasonality: Visual inspection of yearly patterns → set yearly_seasonality Fourier terms 6. Outliers: Winsorize extreme values (> 3σ) or model with robust likelihood 7. Structural breaks: Regime changes (COVID, brand events) → add control dummies

Spend Share Computation

Spend shares inform saturation_beta priors:

spend_shares = X[channel_columns].sum() / X[channel_columns].sum().sum()
# Normalize to use as prior sigma for HalfNormal
spend_shares = spend_shares.values

Pre-Modeling Validation Script

Use scripts/data_validator.py for automated checks:

python scripts/data_validator.py data.csv

Returns quality score (0-100). Target ≥ 75 before modeling.

Transformation Order Decision

--- See also: model_specification.md, diagnostics_benchmarks.md

Diagnostics & Industry Benchmarks

Convergence Diagnostics

Critical Checks

R-hat Interpretation

Model Validity

Posterior Predictive Metrics

Business Logic Guards

ROI Benchmarks by Channel

Sources: Nielsen, Kantar, Triple Whale, industry meta-analyses. Always search web for current benchmarks specific to country/industry/channel.

Adstock Typical Values

Half-life: ln(0.5) / ln(alpha)

Saturation Sanity Checks

Pre/Post Modeling Checklist

## Pre-Modeling
- [ ] Data quality score ≥ 75 (data_validator.py)
- [ ] ≥ 104 weeks of data
- [ ] Media variation CV > 0.3
- [ ] No high multicollinearity (VIF < 5)
- [ ] Prior predictive checks passed
- [ ] Lift tests added (if available)

## Post-Modeling
- [ ] 0 divergences
- [ ] R-hat < 1.01 for all parameters
- [ ] ESS > 400 for all parameters
- [ ] Posterior predictive R² > 0.70
- [ ] Baseline 20-80%
- [ ] Posterior predictive visually reasonable

## Business Validation
- [ ] ROI within industry benchmarks
- [ ] Adstock half-lives reasonable
- [ ] No single channel > 50% contribution
- [ ] Saturation points sensible
- [ ] Prior vs posterior contribution shares differ (data is informative)

Mathematical Formulas

Adstock (Geometric Decay)

AdStock(x_t; α, L) = Σ(s=0 to L)[α^s * x_{t-s}] / Σ(s=0 to L)[α^s]

LogisticSaturation

f(x; λ) = 1 - exp(-λx)

Hill Saturation

Hill(x; ec, slope) = 1 / (1 + (x/ec)^(-slope))

Full Model

KPI_g,t = intercept_g + trend_t + seasonality_t
        + Σ_m[β_g,m * Saturation(Adstock(x_g,m,t))]
        + Σ_c[γ_c * control_c,t] + ε_g,t

ROI

ROI_m = incremental_contribution_m / total_spend_m
mROI_m = ∂contribution/∂spend_m |_{current_spend}

Optimality Condition

At optimal allocation: mROI_1 = mROI_2 = ... = mROI_M (marginal ROIs equalized)

--- See also: model_fit.md, media_deep_dive.md

Lift Test Calibration

Lift tests resolve causal identification when channels are correlated.

Adding Lift Tests

# After build_model, before fit
mmm.build_model(X, y)
mmm.add_lift_test_measurements(df_lift_test)
mmm.fit(X=X, y=y, nuts_sampler="nutpie", ...)

Data Format

Example

df_lift = pd.DataFrame({
    "channel": ["tv", "social"],
    "x": [50000, 30000],
    "delta_x": [10000, 5000],
    "delta_y": [25000, 8000],
    "sigma": [5000, 2000],
})

Sigma Estimation

If sigma is unavailable from the experiment:

• Use standard error from the experiment's statistical test
• Conservative: sigma = 0.3 * abs(delta_y) (30% relative uncertainty)
• From confidence interval: sigma = (CI_upper - CI_lower) / (2 * 1.96)

When to Add Lift Tests

• Channels are highly correlated (VIF > 5)
• Prior and posterior contribution shares are nearly identical (data not informative)
• ROAS estimates seem implausible
• Stakeholders require causal identification

Calibrated vs Uncalibrated

Compare models with and without lift tests:

# Uncalibrated
mmm_uncal = MMM(...)
mmm_uncal.build_model(X, y)
mmm_uncal.fit(X=X, y=y, ...)

# Calibrated
mmm_cal = MMM(...)
mmm_cal.build_model(X, y)
mmm_cal.add_lift_test_measurements(df_lift)
mmm_cal.fit(X=X, y=y, ...)

Lift tests typically narrow posterior uncertainty and resolve channel attribution conflicts.

--- See also: model_specification.md, model_fit.md

Media Deep Dive: Contributions, ROAS, and Diagnostics

Channel Contributions Over Time

mmm.plot.contributions_over_time(
    var=["channel_contribution_original_scale",
         "control_contribution_original_scale",
         "intercept_contribution_original_scale"],
    combine_dims=True, hdi_prob=0.94,
)

Original Scale Contributions

CRITICAL: Call add_original_scale_contribution_variable BEFORE sample_posterior_predictive:

mmm.add_original_scale_contribution_variable(
    var=["channel_contribution", "control_contribution",
         "intercept_contribution", "yearly_seasonality_contribution", "y"]
)
mmm.sample_posterior_predictive(X=X, random_seed=rng)

Without this, *_original_scale variables will NOT exist.

Waterfall Decomposition

mmm.plot.waterfall_components_decomposition()

ROAS Computation

Element-wise (quick approximation)

channel_contrib = mmm.idata["posterior"]["channel_contribution_original_scale"]
sum_dims = ["date"] + (["geo"] if "geo" in channel_contrib.dims else [])
contrib_total = channel_contrib.sum(dim=sum_dims)

spend_xr = xr.DataArray(X[channel_columns].sum().values,
                         dims=["channel"], coords={"channel": channel_columns})
roas_samples = contrib_total / spend_xr
az.plot_forest(roas_samples, combined=True, hdi_prob=0.94)

Incremental (preferred — accounts for adstock carryover)

roas = mmm.incrementality.contribution_over_spend(frequency="all_time")
az.plot_forest(roas, combined=True)

Incremental Analysis (mmm.incrementality)

Counterfactual analysis that correctly handles adstock carryover. Preferred approach.

# Total vs Marginal ROAS
total_roas = mmm.incrementality.contribution_over_spend(frequency="all_time")
marginal_roas = mmm.incrementality.marginal_contribution_over_spend(
    frequency="all_time", spend_increase_pct=0.01)

Total ROAS: Overall return per dollar. Marginal ROAS: Return on the next dollar — directly informs reallocation.

Saturation Curves

Scatterplot (direct/marginal)

mmm.plot.saturation_scatterplot(original_scale=True)

Shows actual spend vs. saturated effect at each observed time point.

Smooth posterior curves

curve = mmm.saturation.sample_curve(mmm.idata.posterior, max_value=2)
mmm.plot.saturation_curves(curve, original_scale=True)

Sensitivity Analysis

Counterfactual response under spend scaling:

sweeps = np.linspace(0, 1.5, 16)
mmm.sensitivity.run_sweep(
    sweep_values=sweeps, var_input="channel_data",
    var_names="channel_contribution_original_scale", extend_idata=True,
)
mmm.plot.sensitivity_analysis(hue_dim="channel", x_sweep_axis="relative")

Interpretation: Slope at x=1.0 (current spend) = marginal efficiency. Flat = saturated.

Channel Contribution Share

mmm.plot.channel_contribution_share_hdi(hdi_prob=0.94)

Compare prior vs posterior shares to assess data informativeness. If identical → channels not identifiable → add lift tests.

Summary DataFrames

mmm.summary.posterior_predictive()
mmm.summary.contributions()
mmm.summary.roas()
mmm.summary.channel_spend()
mmm.summary.saturation_curves()
mmm.summary.adstock_curves()
mmm.summary.total_contribution()
mmm.summary.change_over_time()

--- See also: budget_optimization.md, lift_test_calibration.md, plot_api.md

Model Fit & Diagnostics

Fitting

Always fit on the full dataset:

mmm.fit(X=X, y=y, target_accept=0.9, chains=6, draws=800,
        tune=1_500, nuts_sampler="nutpie", random_seed=rng)
mmm.sample_posterior_predictive(X=X, random_seed=rng)

Diagnostics Checklist

1. Divergences (must be 0)

mmm.idata["sample_stats"]["diverging"].sum().item()

2. R-hat (must be < 1.01)

az.summary(data=mmm.idata, var_names=[...])["r_hat"].describe()

3. Trace plots

az.plot_trace(data=mmm.fit_result, var_names=[...], compact=True)

4. ESS

az.summary(data=mmm.idata, var_names=[...])["ess_bulk"].describe()

5. Posterior predictive

mmm.plot.posterior_predictive(hdi_prob=0.94)
mmm.plot.residuals_over_time()

Diagnostic Thresholds

Time-Slice Cross-Validation

For stability assessment before the final fit:

from pymc_marketing.mmm.time_slice_cross_validation import TimeSliceCrossValidator

cv = TimeSliceCrossValidator(
    n_init=163, forecast_horizon=12,
    date_column="date", step_size=1,
)
results = cv.run(X, y, sampler_config={...}, yaml_path=...)

cv.plot.param_stability(results, parameter=["adstock_alpha"], dims={...})
cv.plot.cv_predictions(results)
cv.plot.cv_crps(results)

Key: CV is for assessing stability, NOT for selecting the model. The final model is always fit on all data.

Common Issues

High divergences

• Increase target_accept (0.9 → 0.95 → 0.99)
• Use more informative priors (especially on saturation params)
• Check for near-zero spend channels → remove or combine

R-hat > 1.01

• Increase draws and tune
• Check for multimodality (bimodal trace plots)
• Reparameterize: use non-centered for hierarchical (set centered=False)

Poor in-sample fit

• Add control variables (promotions, competitors, holidays)
• Check yearly_seasonality setting (try 3-7 Fourier pairs)
• Consider time_varying_intercept=True for irregular trends
• Check target scaling (log transform if heavily skewed)

Implausible ROAS

• Add lift test calibration (see lift_test_calibration.md)
• Check for multicollinearity between channels
• Verify spend data units (impressions vs. dollars)
• Review saturation_beta priors (too wide → unstable)

--- See also: model_specification.md, media_deep_dive.md, diagnostics_benchmarks.md

Model Specification

Table of Contents

• MMM Constructor Reference
• Adstock Transformations
• Saturation Transformations
• Prior Specification
• Model Config Dictionary
• Multidimensional / Hierarchical Models
• Scaling Configuration
• Prior Predictive Checks
• Building the Model

MMM Constructor Reference

from pymc_marketing.mmm.multidimensional import MMM

mmm = MMM(
    date_column="date",
    channel_columns=["tv", "radio", "social"],
    target_column="y",
    adstock=GeometricAdstock(l_max=6),
    saturation=LogisticSaturation(),
    dims=None,                              # tuple[str, ...] for multidimensional
    scaling=None,                           # dict or Scaling object
    model_config=None,                      # dict of Prior objects
    sampler_config=None,                    # dict of sampler kwargs
    control_columns=None,                   # list[str] for control variables
    yearly_seasonality=None,                # int: number of Fourier terms
    adstock_first=True,                     # adstock before saturation
    time_varying_intercept=False,           # bool or HSGPBase
    time_varying_media=False,               # bool or HSGPBase
)

Adstock Transformations

GeometricAdstock

from pymc_marketing.mmm import GeometricAdstock
from pymc_extras.prior import Prior

adstock = GeometricAdstock(l_max=6, normalize=True)

# Custom priors per channel
adstock = GeometricAdstock(
    l_max=6,
    priors={"alpha": Prior("Beta", alpha=2, beta=5, dims="channel")},
)

# Per-channel-per-geo
adstock = GeometricAdstock(
    l_max=6,
    priors={"alpha": Prior("Beta", alpha=2, beta=5, dims=("channel", "geo"))},
)

Default prior: alpha ~ Beta(1, 3) — favors fast decay.

l_max rule of thumb: Max plausible carryover. Weekly: l_max=6. Daily: l_max=14+. Too large wastes computation; too small truncates real effects.

Saturation Transformations

LogisticSaturation (most common)

from pymc_marketing.mmm import LogisticSaturation

saturation = LogisticSaturation(
    priors={
        "lam": Prior("Gamma", alpha=3, beta=1, dims="channel"),
        "beta": Prior("HalfNormal", sigma=spend_shares, dims="channel"),
    },
)

beta controls max reachable effect per channel. Setting prior proportional to spend shares encodes that higher-spend channels should have proportionally larger effects.

All Saturation Functions

Prior Specification

All priors use pymc_extras.prior.Prior:

from pymc_extras.prior import Prior

Prior("Normal", mu=0, sigma=1)                              # Simple
Prior("HalfNormal", sigma=0.5, dims="channel")              # With dims
Prior("Beta", alpha=2, beta=5, dims=("channel", "geo"))     # Multiple dims

# Hierarchical (non-centered) via LogNormalPrior
from pymc_marketing.special_priors import LogNormalPrior

LogNormalPrior(
    mean=Prior("Gamma", mu=1.0, sigma=1.0),
    std=Prior("HalfNormal", sigma=1.0),
    dims=("channel", "geo"),
    centered=False,
)

Model Config Dictionary

model_config = {
    "intercept": Prior("Normal", mu=0.2, sigma=0.05),
    "saturation_beta": Prior("HalfNormal", sigma=spend_shares, dims="channel"),
    "gamma_control": Prior("Normal", mu=0, sigma=1, dims="control"),
    "gamma_fourier": Prior("Laplace", mu=0, b=1, dims="fourier_mode"),
    "likelihood": Prior("TruncatedNormal", lower=0, sigma=Prior("HalfNormal", sigma=1)),
}

Multidimensional / Hierarchical Models

Pooling Strategies

Partial pooling (recommended):

model_config = {
    "saturation_beta": LogNormalPrior(
        mean=Prior("Gamma", mu=1.0, sigma=1.0),
        std=Prior("HalfNormal", sigma=1.0),
        dims=("channel", "geo"), centered=False,
    ),
}

Full pooling — same params across all geos:

model_config = {"saturation_beta": Prior("HalfNormal", sigma=spend_shares, dims="channel")}

No pooling — independent per geo (needs substantial data):

model_config = {"saturation_beta": Prior("HalfNormal", sigma=0.5, dims=("geo", "channel"))}

Scaling Configuration

# Global max across geos (preferred — keeps channels comparable)
scaling_config = {
    "channel": {"method": "max", "dims": ()},
    "target": {"method": "max", "dims": ()},
}

# Per-geo scaling
scaling_config = {
    "channel": {"method": "max", "dims": ("geo",)},
    "target": {"method": "max", "dims": ("geo",)},
}

Prior Predictive Checks

Always run before fitting:

# CRITICAL: call add_original_scale_contribution_variable BEFORE prior predictive
mmm.add_original_scale_contribution_variable(
    var=["channel_contribution", "control_contribution",
         "intercept_contribution", "yearly_seasonality_contribution", "y"]
)
mmm.sample_prior_predictive(X, y, samples=4_000, random_seed=42)

Check that prior contributions are not concentrated on a single channel. If one channel dominates, widen saturation_beta priors.

Building the Model

mmm.build_model(X, y)
mmm.graphviz()    # DAG visualization
mmm.table()       # Rich summary table of all variables

After building, optionally add lift tests before fitting:

mmm.add_lift_test_measurements(df_lift_test)

--- See also: data_analysis.md, model_fit.md, lift_test_calibration.md

Plot API Reference

All via mmm.plot namespace.

--- See also: media_deep_dive.md, budget_optimization.md

Scenario Planning (Meridian-Inspired)

Forward-looking budget optimization patterns adapted from Google Meridian's Scenario Planner for use with PyMC-Marketing. These patterns go beyond single-point optimization to enable interactive, multi-scenario planning.

Core Concepts (from Meridian)

Key Meridian insight: MMM forecasts incremental outcome, not total outcome. Control variables and temporal effects cancel in the counterfactual difference. This means scenario planning should focus on incremental contribution, not total KPI prediction.

Pattern 1: Fixed Budget Optimization

Find optimal allocation for a given budget:

allocation, result = optimizer.optimize_budget(
    budget=1_000_000,
    budget_bounds=bounds,
    minimize_kwargs={"method": "SLSQP", "options": {"ftol": 1e-4, "maxiter": 10_000}},
)
response = optimizer.sample_response_distribution(
    allocation_strategy=allocation, include_carryover=True)

Pattern 2: Flexible Budget (Target ROI)

Find maximum budget that maintains a target ROAS. PyMC-Marketing doesn't have a direct target_roi parameter, so implement via binary search:

def find_max_budget_at_target_roas(optimizer, bounds, target_roas,
                                     budget_range=(100_000, 5_000_000), tol=10_000):
    """Binary search for max budget maintaining target ROAS."""
    low, high = budget_range

    while high - low > tol:
        mid = (low + high) / 2
        alloc, res = optimizer.optimize_budget(budget=mid, budget_bounds=bounds,
            minimize_kwargs={"method": "SLSQP", "options": {"ftol": 1e-4}})

        # Sample response and compute ROAS
        resp = optimizer.sample_response_distribution(
            allocation_strategy=alloc, include_carryover=True)
        total_contrib = resp["posterior_predictive"]["y"].sum(dim="date").mean(dim=("chain", "draw"))
        roas = total_contrib.item() / mid

        if roas >= target_roas:
            low = mid  # Can afford more budget
        else:
            high = mid  # Budget too high, ROAS drops below target

    return low

Pattern 3: Budget Sweep (Efficient Frontier)

Evaluate optimal response across budget levels to identify diminishing returns:

import pandas as pd

budget_levels = np.linspace(200_000, 2_000_000, 15)
frontier = []

for budget in budget_levels:
    alloc, res = optimizer.optimize_budget(budget=budget, budget_bounds=bounds,
        minimize_kwargs={"method": "SLSQP", "options": {"ftol": 1e-4, "maxiter": 10_000}})
    resp = optimizer.sample_response_distribution(
        allocation_strategy=alloc, include_carryover=True)

    total_response = resp["posterior_predictive"]["y"].sum(dim="date")
    mean_resp = total_response.mean(dim=("chain", "draw")).item()
    hdi = az.hdi(total_response, hdi_prob=0.90)

    # Channel allocation breakdown
    ch_alloc = {ch: alloc.sel(channel=ch).sum().item() for ch in channel_columns}

    frontier.append({
        "budget": budget,
        "expected_response": mean_resp,
        "roas": mean_resp / budget,
        "hdi_low": hdi.values[0],
        "hdi_high": hdi.values[1],
        **ch_alloc,
    })

frontier_df = pd.DataFrame(frontier)

Interpretation: The point where the efficient frontier curve flattens is the budget saturation point — spending beyond this yields minimal incremental return.

Pattern 4: Multi-Scenario Comparison

Create named scenarios for stakeholder comparison:

scenarios = {
    "Conservative": {
        "budget": 800_000,
        "bounds_multiplier": (0.8, 1.2),  # ±20% from current
    },
    "Moderate": {
        "budget": 1_000_000,
        "bounds_multiplier": (0.5, 1.5),  # ±50%
    },
    "Aggressive": {
        "budget": 1_500_000,
        "bounds_multiplier": (0.3, 2.0),  # wide reallocation
    },
}

results = {}
for name, config in scenarios.items():
    lo, hi = config["bounds_multiplier"]
    scenario_bounds = xr.DataArray(
        data=np.array([[lo, hi]] * len(channel_columns)) * equal_share,
        dims=["channel", "bound"],
        coords={"channel": channel_columns, "bound": ["lower", "upper"]},
    )
    alloc, res = optimizer.optimize_budget(
        budget=config["budget"], budget_bounds=scenario_bounds,
        minimize_kwargs={"method": "SLSQP", "options": {"ftol": 1e-4}},
    )
    resp = optimizer.sample_response_distribution(
        allocation_strategy=alloc, include_carryover=True)
    results[name] = {"allocation": alloc, "response": resp, "config": config}

Pattern 5: Cost-Per-Media-Unit Sensitivity

When CPM/CPC is expected to change, re-evaluate without refitting:

def evaluate_cpm_change(mmm, X, channel, cpm_multiplier):
    """What happens if CPM for a channel increases/decreases?"""
    X_modified = X.copy()
    # Same budget buys fewer impressions → lower spend-equivalent
    X_modified[channel] = X_modified[channel] / cpm_multiplier
    resp = mmm.sample_posterior_predictive(X_modified, extend_idata=False, random_seed=42)
    return resp

# CPM increases 30% on social → same budget buys fewer units
resp_high_cpm = evaluate_cpm_change(mmm, X, "social", 1.3)

Pattern 6: Quarterly/Monthly Planning Horizon

Adapt Meridian's time-breakdown approach:

import pandas as pd

planning_start = pd.Timestamp("2026-01-01")
quarters = pd.date_range(planning_start, periods=4, freq="QS")

quarterly_plans = []
for q_start in quarters:
    q_end = q_start + pd.offsets.QuarterEnd()
    q_optimizer = MultiDimensionalBudgetOptimizerWrapper(
        model=mmm, start_date=str(q_start), end_date=str(q_end))

    alloc, res = q_optimizer.optimize_budget(
        budget=quarterly_budget, budget_bounds=bounds,
        minimize_kwargs={"method": "SLSQP", "options": {"ftol": 1e-4}})

    quarterly_plans.append({
        "quarter": f"Q{(q_start.month-1)//3 + 1} {q_start.year}",
        "allocation": alloc, "start": q_start, "end": q_end,
    })

Pattern 7: Spend Shift Ratio Bounds (Meridian Style)

Express bounds as % shift from historical spend rather than absolute values:

def meridian_style_bounds(historical_spend, min_shift=0.3, max_shift=0.3):
    """
    Meridian convention: min_shift=0.3 means spend can decrease by 30%.
    max_shift=0.3 means spend can increase by 30%.
    """
    lower = (1 - min_shift) * historical_spend
    upper = (1 + max_shift) * historical_spend
    return xr.DataArray(
        data=np.stack([lower.values, upper.values], axis=-1),
        dims=["channel", "bound"],
        coords={"channel": historical_spend.index.tolist(), "bound": ["lower", "upper"]},
    )

historical = X[channel_columns].sum()
bounds = meridian_style_bounds(historical, min_shift=0.5, max_shift=1.0)

Reporting Template

Generate a scenario comparison table for stakeholder decks:

Use scripts/report_generator.py to produce HTML reports or mmm.plot methods for inline visualizations.

--- See also: budget_optimization.md, media_deep_dive.md

Time-Varying Parameters

GP-based time-varying intercept and media multiplier via Hilbert Space Gaussian Processes (HSGP).

When to Use TVP

Use when residuals show irregular, non-repeating temporal variation NOT explained by seasonality, trend, or controls. The GP is for in-sample decomposition; it reverts to prior mean out of sample.

Configuration

from pymc_marketing.hsgp_kwargs import HSGPKwargs

mmm = MMM(
    ...,
    time_varying_intercept=True,
    time_varying_media=True,
    model_config={
        "intercept_tvp_config": HSGPKwargs(
            m=500, L=188, eta_lam=5.0, ls_mu=5.0, ls_sigma=10.0),
        "media_tvp_config": HSGPKwargs(
            ls_mu=11.0, ls_sigma=5.0),
    },
)

HSGPKwargs Parameters

Diagnostics

• Check that TVP doesn't absorb media effects (compare contributions with/without TVP)
• Length scale > data length → TVP is essentially constant (unnecessary)
• Length scale < 4 weeks → TVP may be overfitting to noise

Lift Tests with TVP

When time_varying_media=True, include date in lift test DataFrame so each measurement maps to the correct media_temporal_latent_multiplier time coordinate.

--- See also: model_specification.md, model_fit.md

Trade Marketing Decomposition Framework

Separates media-driven sales from trade-driven sales and base demand. This is the layer that exposes the attribution gap holdco MMMs hide.

The Problem

Standard MMMs model: KPI = baseline + Σ(media) + ε

Reality for most brands: KPI = base_demand + media_driven + trade_driven + interaction + unexplained

Trade marketing (promotions, distributor incentives, shelf placement, co-op programs) drives 40-60% of CPG/FMCG sales but is either omitted or reduced to a binary control variable in agency MMMs. This systematically overattributes to media.

Flexible Data Schema

Since client data varies wildly, the schema adapts to what's available:

TRADE_VARIABLES = {
    # Tier 1: Almost always available
    "promo_flag": {
        "type": "binary",
        "description": "Was a promotion active?",
        "functional_form": "step",
        "availability": "high",
    },
    "promo_depth": {
        "type": "continuous",
        "description": "Discount depth (% off or absolute)",
        "functional_form": "linear_with_ceiling",
        "availability": "high",
    },
    "price_index": {
        "type": "continuous",
        "description": "Price relative to category average or competitor",
        "functional_form": "log_linear",
        "availability": "high",
    },

    # Tier 2: Often available for CPG/FMCG
    "distribution_acv": {
        "type": "continuous",
        "description": "All Commodity Volume — % of stores carrying the product",
        "functional_form": "log_linear",
        "availability": "medium",
    },
    "weighted_distribution": {
        "type": "continuous",
        "description": "Distribution weighted by store revenue",
        "functional_form": "log_linear",
        "availability": "medium",
    },
    "feature_flag": {
        "type": "binary",
        "description": "Featured in retailer circular/flyer",
        "functional_form": "step_with_decay",
        "availability": "medium",
    },
    "display_flag": {
        "type": "binary",
        "description": "End-cap or special display placement",
        "functional_form": "step",
        "availability": "medium",
    },

    # Tier 3: Available for sophisticated clients
    "sell_in_volume": {
        "type": "continuous",
        "description": "Volume shipped to distributors/retailers",
        "functional_form": "linear",
        "availability": "low",
    },
    "shelf_share": {
        "type": "continuous",
        "description": "Share of shelf space vs competitors",
        "functional_form": "log_linear",
        "availability": "low",
    },
    "pos_sell_out": {
        "type": "continuous",
        "description": "Point-of-sale actual consumer purchases",
        "functional_form": "target_variable",
        "availability": "low",
    },
    "trade_spend": {
        "type": "continuous",
        "description": "Total trade marketing investment",
        "functional_form": "hill_saturation",
        "availability": "low",
    },
    "ppa_price_per_agreement": {
        "type": "continuous",
        "description": "Price-per-agreement with retailer/distributor",
        "functional_form": "linear",
        "availability": "low",
    },
    "co_op_spend": {
        "type": "continuous",
        "description": "Co-operative advertising with retailers",
        "functional_form": "linear_with_ceiling",
        "availability": "low",
    },

    # Tier 4: Best-in-class data infrastructure
    "store_visits": {
        "type": "continuous",
        "description": "Foot traffic or store visit data",
        "functional_form": "funnel_input",
        "availability": "rare",
    },
    "digital_to_store": {
        "type": "continuous",
        "description": "Online-influenced offline purchases",
        "functional_form": "funnel_input",
        "availability": "rare",
    },
    "shopper_marketing_spend": {
        "type": "continuous",
        "description": "In-store marketing investment",
        "functional_form": "linear_with_ceiling",
        "availability": "rare",
    },
}

Functional Forms

Trade variables need different response functions than media:

Step Function (promotions)

def step_response(promo_flag, lift_magnitude):
    """Promotion drives instant volume lift, no carryover."""
    return promo_flag * lift_magnitude

Linear with Ceiling (trade spend)

def linear_with_ceiling(x, slope, ceiling):
    """Linear response that caps at a maximum effect."""
    return np.minimum(slope * x, ceiling)

Step with Decay (features/displays)

def step_with_decay(flag, lift, decay, max_lag=4):
    """Promotion lifts then decays — like adstock but typically faster."""
    result = np.zeros_like(flag, dtype=float)
    for t in range(len(flag)):
        if flag[t]:
            for s in range(min(max_lag, len(flag) - t)):
                result[t + s] += lift * (decay ** s)
    return result

Log-Linear (distribution/price)

def log_linear_response(x, elasticity):
    """Log-linear price/distribution elasticity."""
    return elasticity * np.log(x)

Model Specification with Trade Variables

In PyMC-Marketing (as control variables with informed priors)

# Trade variables as controls with domain-informed priors
model_config = {
    # Media priors (from context layer)
    "saturation_beta": Prior("HalfNormal", sigma=spend_shares, dims="channel"),

    # Trade control priors — these are NOT saturated, they're linear/step
    "gamma_control": Prior("Normal", mu=trade_prior_means, sigma=trade_prior_sds, dims="control"),
}

# Control columns include trade variables
control_columns = ["promo_depth", "price_index", "distribution_acv", "feature_flag"]

Custom Model (when trade needs its own response curves)

For sophisticated decomposition where trade variables need non-linear treatment beyond what control_columns supports:

import pymc as pm
from pymc_marketing.mmm import GeometricAdstock, LogisticSaturation

with pm.Model(coords=coords) as trade_mmm:
    # --- Media component (standard) ---
    channel_data_ = pm.Data("channel_data", channel_scaled, dims=("date", "channel"))
    adstocked = adstock.apply(channel_data_, dims=("channel",))
    media_contribution = saturation.apply(adstocked, dims=("channel",))
    total_media = pm.math.sum(media_contribution, axis=-1)

    # --- Trade component (custom response curves) ---
    promo_depth = pm.Data("promo_depth", promo_data, dims=("date",))
    promo_beta = pm.HalfNormal("promo_beta", sigma=0.5)
    promo_ceiling = pm.HalfNormal("promo_ceiling", sigma=2.0)
    trade_promo = pm.math.minimum(promo_beta * promo_depth, promo_ceiling)

    price_idx = pm.Data("price_index", price_data, dims=("date",))
    price_elasticity = pm.Normal("price_elasticity", mu=-1.5, sigma=0.5)  # Negative = higher price → lower sales
    trade_price = price_elasticity * pm.math.log(price_idx)

    dist_acv = pm.Data("distribution_acv", acv_data, dims=("date",))
    dist_elasticity = pm.HalfNormal("dist_elasticity", sigma=0.5)
    trade_dist = dist_elasticity * pm.math.log(dist_acv)

    total_trade = trade_promo + trade_price + trade_dist

    # --- Interaction: media × trade ---
    interaction_beta = pm.Normal("interaction_beta", mu=0, sigma=0.1)
    interaction = interaction_beta * total_media * trade_promo

    # --- Baseline ---
    intercept = pm.Normal("intercept", mu=0.5, sigma=0.1)
    # + trend + seasonality as needed

    # --- Likelihood ---
    mu = intercept + total_media + total_trade + interaction
    sigma = pm.HalfNormal("sigma", sigma=0.5)
    y_obs = pm.Normal("y", mu=mu, sigma=sigma, observed=y_scaled, dims=("date",))

Decomposition Output

The framework produces a four-way decomposition:

┌──────────────────────────────────────────────┐
│              TOTAL KPI (100%)                │
├──────────────────────────────────────────────┤
│                                              │
│  ┌────────────┐  ┌────────────┐              │
│  │ Base Demand│  │ Media-     │              │
│  │ (Intercept │  │ Driven     │              │
│  │ + Trend +  │  │ Sales      │              │
│  │ Seasonality│  │            │              │
│  │ )          │  │ TV: 8%     │              │
│  │            │  │ Digital: 5%│              │
│  │ 55%        │  │ Social: 3% │              │
│  │            │  │ Search: 4% │              │
│  │            │  │            │              │
│  │            │  │ Total: 20% │              │
│  └────────────┘  └────────────┘              │
│                                              │
│  ┌────────────┐  ┌────────────┐              │
│  │ Trade-     │  │ Interaction│              │
│  │ Driven     │  │ (Media ×   │              │
│  │ Sales      │  │ Trade)     │              │
│  │            │  │            │              │
│  │ Promo: 10% │  │ 3%         │              │
│  │ Price: 5%  │  │            │              │
│  │ Dist: 7%   │  │            │              │
│  │            │  │            │              │
│  │ Total: 22% │  │            │              │
│  └────────────┘  └────────────┘              │
│                                              │
└──────────────────────────────────────────────┘

KPI Adaptation by Industry

Funnel Decomposition (Store Visits + Digital Sales)

When store visit or digital-to-store data is available, add a funnel layer:

Media Spend → Awareness (impressions/reach)
                 ↓
             Consideration (search, site visits)
                 ↓
             ┌───┴───┐
         Store Visit  Digital Visit
             ↓            ↓
         In-Store     Online
         Purchase     Purchase
             ↑            ↑
        Trade Marketing   Digital Trade
        (shelf, promo)    (coupons, delivery promos)

This doesn't require a separate model — it's a post-hoc allocation using observed funnel ratios:

# Media-attributable store visits
media_driven_visits = total_visits * (media_contribution_share)
trade_driven_visits = total_visits * (trade_contribution_share)
organic_visits = total_visits * (baseline_share)

# Per-visit conversion rate
visit_to_purchase_rate = total_purchases / total_visits

# Final attribution
media_driven_sales = media_driven_visits * visit_to_purchase_rate * avg_ticket
trade_driven_sales = trade_driven_visits * visit_to_purchase_rate * avg_ticket

The Holdco Exposure Play

When presenting to a CMO whose agency runs the MMM:

1. Show the decomposition gap: "Your agency's MMM shows 40% media contribution. Our analysis shows 20% media, 22% trade, 58% baseline." 2. Quantify the overattribution: "This means ~$X million of your media budget is being justified by sales that would have happened anyway." 3. Show the interaction: "Media and trade work together — cutting trade by 20% would reduce media ROI by 15% because in-store presence amplifies ad recall." 4. Recommend rebalancing: "Shifting $Y from underperforming media channels to trade co-op programs would increase total sales by Z%."

--- See also: context_layer.md, acid_test_validation.md, scenario_planning.md