MMM Case Study II: Funnel-Aware Modeling, ROAS Calibration, and Budget Planning#
This notebook is the second iteration of the MMM End-to-End Case Study notebook. The first iteration left us with a working baseline model and a first set of insights and decisions: the baseline drives about \(60\%\) of sales, Online Display and Insert show the highest return on ad spend (ROAS) but with wide uncertainty bands, and Direct Mail takes the largest budget while showing the lowest ROAS. A counterfactual re-allocation suggested a meaningful uplift in media-driven sales.
In many situations, these initial results are not immediately used for decision-making. Instead, the team has to present them to stakeholders, who come back with questions. For example, the stakeholders might ask:
Why is search engine marketing (SEM) not in the model? SEM is the second largest line in the media budget.
How stable are these ROAS numbers over time, and can we tighten the wide ones?
This iteration answers both.
The analysis proceeds in nine main steps (model specification, diagnostics and the appendix sit between them):
The SEM Funnel
ROAS Calibration from Experiments
Time-Slice Cross-Validation
Final Model, Diagnostics and Posterior Predictive Checks
Media Deep Dive
Contribution Decomposition
Direct and Mediated Contributions
Return on Ad Spend (ROAS)
Budget Optimization
In-Sample Uplift
Comparison with the First Iteration
Conclusion and Next Steps
Note
This notebook focuses on the business application. For the model building blocks see Media Mix Modeling (MMM) in Python: End-to-End Example, and for the funnel machinery in depth see Funnel-Aware MMM: Upper to Lower-Funnel Mediation via a custom MuEffect and Advanced Funnel-Aware MMM: Geo-Level Mediation Through Lower-Funnel Spend. The appendix shows the same model built from a YAML registry and checks that it builds the same graph.
The SEM Funnel#
Why did the first iteration deliberately leave the SEM spend column out? For an initial (flat) baseline model, that was the right call. Paid search spend reacts to search demand, and search demand reacts to upper funnel advertising. So if we add SEM as an eighth flat channel, the model conditions on a variable that sits in the middle of the causal path. It blocks the indirect path and credits SEM with sales that TV, Radio, and Newspaper created. This is the first failure mode analyzed in Funnel-Aware MMM: Upper to Lower-Funnel Mediation via a custom MuEffect.
This is discussed in the paper “Challenges And Opportunities In Media Mix Modeling”, in section \(4.3.2\) on funnel effects:
Selection bias arises when there are funnel effects in the media and the model is mis-specified. When an ad channel also impacts the level of another ad channel, using a model like in Equation (3), which simultaneously estimate the impact of all ad channels in one equation, will lead to biased estimates. An example is a TV campaign driving more related queries, which in turn increase volume of paid search ads. For a derivation of the size of the bias, see Angrist and Krueger (1999). The underlying reason for the bias in the mis-specified model is that the downstream ads were affected by the TV ads. In assessing the ROAS of TV, the linear regression model does not account for the changes in paid search ads caused by TV.
Downstream ads should not be included with exogenously-determined ads in a single regression equation. Alternatives include graphical models (Pearl (2009)) and structural equation models. However, the problem isn’t only a matter of model form. Both these alternatives require estimating the causal effect of the upstream ads on the downstream ads, which has just as stringent data requirements as estimating the effect of an ad channel on sales.
Now let’s come back to our specific use case:
Business Context: SEM is \(26.9\%\) of the modeled media spend, the second largest line after Direct Mail. A mediator-blind model undervalues TV, Radio, and Newspaper exactly where they earn part of their return, and the resulting budget shift can destroy the demand SEM lives on.
The dataset gives us two SEM columns to work with. sem_spend is the money the search team spends, an exogenous budget lever. sem_impressions is how often the ads showed, which scales with search volume, so it is a noisy measurement of harvested demand. We use both: the spend as a decision variable, the impressions as evidence.
Model Structure#
The model adds one latent demand pool \(D_t\) and one mediator \(M^{*}_t\) on top of the first iteration structure:
and the sales mean gains the term \(\text{sat}^{lf}\left(\text{adstock}^{lf}(M^{*}_t)\right)\). The two fixed scaling constants are the maxima of the two SEM series, \(s_{\text{spd}} = \max_t(\text{SEM spend}_t)\) and \(s_{\text{imp}} = \max_t(\text{SEM impressions}_t)\).
How should you read the mediator? Because the measurement coefficient is fixed at one times \(s_{\text{imp}}\), it is expressed in units of the peak week of observed impressions. That is, \(M^{*}_t = 1\) means “as much search demand as the busiest week in the data”.
\(D_t\) is the pool of people looking for the product in week \(t\), in those units. Its first part, the baseline, is ambient interest that exists whatever the brand does. Its second part is the interest that TV, Radio, and Newspaper added, arriving with a lag (the adstock) and with diminishing returns (the saturation), so a doubled TV budget does not double the pool.
\(M^{*}_t\) is then the demand that paid search actually touches: the pool \(D_t\) plus what the search team’s own spend adds on top of it, because bidding on more queries and higher positions reaches people the organic pool would have missed.
Later in the notebook we measure how the mediator splits between these three sources (ambient baseline, upper funnel media, and SEM’s own spend). That split is what decides how much credit the upper funnel channels can claim for search sales.
Note
Why this parametrization?
\(\lambda\) is the exchange rate between search money and search demand: it says how much of a peak-impressions week one peak-spend week buys. Dividing the spend by \(s_{\text{spd}}\) keeps \(\lambda\) unit-free, so a printed value like \(0.7\) reads directly as “seven tenths of a peak week per peak week of spend” and the prior does not depend on the currency scale.
The spend term is linear because a search auction turns money into impressions almost mechanically over the observed range. The exploratory plot of spend against impressions below shows the two series moving together week by week. We do not lose the diminishing returns to SEM in sales by this choice: they enter downstream, where \(\text{sat}^{lf}\) acts on the whole mediator before it reaches the sales mean.
Impressions enter as a measurement. They scale with query volume, so they tell the model where the pool is, while sales are what the pool converts into.
Like Funnel-Aware MMM: Upper to Lower-Funnel Mediation via a custom MuEffect, this parametrization keeps organic demand out of the funnel except through \(D_t\). In reality, an organic spike in brand increase could increase SEM clicks and, in turn inflate sem_spend, without that necessarily being incremental to sales. See {ref}`mmm_funnel_mueffect_advanced for a more advanced DAG representation.
One more modeling decision worth highlighting: the measurement coefficient on \(M^{*}\) is fixed at one. This anchors the scale of the latent mediator. With a free coefficient, the model could double the demand and halve the coefficient and fit the impressions equally well, and the sampler would wander along that ridge.
Finally, all seven channels from the first iteration keep their direct paths. Only TV, Radio, and Newspaper feed the demand pool, while Direct Mail, Insert, Online Display, and Social Media stay direct-only. This matters for the calibration later on: the ROAS calibration constrains the direct path, so we calibrate direct-only channels where the experiment measures exactly what the model constrains.
Here is the resulting causal structure.
import graphviz as gr
dag = gr.Digraph()
dag.attr(rankdir="LR")
dag.node("upper", "TV, Radio,\nNewspaper", style="filled", fillcolor="lightpink")
dag.node(
"direct",
"Direct Mail, Insert,\nOnline Display, Social Media",
style="filled",
fillcolor="lightpink",
)
dag.node("demand", "latent demand D", style="dashed", color="goldenrod")
dag.node("sem_spend", "sem_spend B", style="filled", fillcolor="lightgray")
dag.node("mediator", "mediator M*", style="dashed", color="steelblue")
dag.node("impressions", "sem_impressions", style="filled", fillcolor="skyblue")
dag.node("sales", "sales Y", style="filled", fillcolor="lightgreen")
dag.edge("upper", "sales", label="direct")
dag.edge("direct", "sales", label="direct")
dag.edge("upper", "demand")
dag.edge("demand", "mediator")
dag.edge("sem_spend", "mediator", label="lambda")
dag.edge("mediator", "impressions", label="measurement", style="dashed")
dag.edge("mediator", "sales", label="conversion")
dag
Tip
Filled nodes are observed and dashed nodes are latent. The dashed edge into sem_impressions is a measurement: impressions do not cause anything in this graph, they simply tell us where the mediator is.
Prepare Notebook#
We load the necessary libraries and set the notebook’s configuration.
Note
One convention for the whole notebook: every interval we report or draw is a \(94\%\) highest density interval (HDI). The two ArviZ settings below make that the default for every summary table and every forest or density plot, so a number in the prose, a row in a table, and a bar in a figure always mean the same thing. ArviZ 1.x defaults to an \(89\%\) equal-tailed interval, which is why we set this explicitly.
The one exception is the experiment input in the calibration section, which is a reported \(95\%\) range. We say so where it appears.
import tempfile
import warnings
from pathlib import Path
import arviz as az
import arviz_plots as azp
import matplotlib.pyplot as plt
import matplotlib.ticker as mtick
import numpy as np
import pandas as pd
import pymc as pm
import pymc.dims as pmd
import seaborn as sns
import xarray as xr
import yaml
from pydantic import InstanceOf
from pymc_extras.prior import Prior
from pytensor.compile.mode import Mode
from pytensor.graph.traversal import ancestors
from pymc_marketing.mmm import GeometricAdstock, LogisticSaturation
from pymc_marketing.mmm.additive_effect import DataVarMuEffect, IncrementalitySpec
from pymc_marketing.mmm.builders.factories import REGISTRY
from pymc_marketing.mmm.builders.yaml import build_mmm_from_yaml
from pymc_marketing.mmm.constraints import Constraint, build_default_sum_constraint
from pymc_marketing.mmm.data_conversion import to_mmm_dataset
from pymc_marketing.mmm.media_transformation import MediaTransformation
from pymc_marketing.mmm.mmm import MMM
from pymc_marketing.mmm.time_slice_cross_validation import TimeSliceCrossValidator
from pymc_marketing.paths import data_dir
warnings.filterwarnings("ignore", category=FutureWarning)
# pytensor falls back to object mode for the truncated normal sampler: correct, only slower.
warnings.filterwarnings("ignore", message="Numba will use object mode")
# pytensor formats a diagnostic message with a 0 / 0 relative error on all-zero tensors.
warnings.filterwarnings(
"ignore",
message="invalid value encountered in divide",
module="pytensor.tensor.type",
)
az.style.use("arviz-darkgrid")
HDI_PROB = 0.94
az.rcParams["stats.ci_kind"] = "hdi"
az.rcParams["stats.ci_prob"] = HDI_PROB
plt.rcParams["figure.figsize"] = [12, 7]
plt.rcParams["figure.dpi"] = 100
%config InlineBackend.figure_format = "retina"
seed: int = sum(map(ord, "mmm_case_study2"))
rng: np.random.Generator = np.random.default_rng(seed=seed)
Load Data#
We use the same public dataset as the first iteration, with the same cleaning and channel names. There are just two changes: we keep the two SEM columns, renamed to sem_spend and sem_impressions, and we rename the date column to date (the name the cross-validation plotting utilities expect). As before, mdsp_viddig and mdsp_auddig stay excluded.
Note that there is no train and test split this time. Time-slice cross-validation replaces it, and the final model uses all \(209\) weeks.
data_path = data_dir / "case_study_data.csv"
raw_df = pd.read_csv(data_path)
date_column = "date"
target_column = "sales"
channel_mapping = {
"mdsp_dm": "Direct Mail",
"mdsp_inst": "Insert",
"mdsp_nsp": "Newspaper",
"mdsp_audtr": "Radio",
"mdsp_vidtr": "TV",
"mdsp_so": "Social Media",
"mdsp_on": "Online Display",
}
channel_columns = sorted(channel_mapping.values())
funnel_mapping = {"mdsp_sem": "sem_spend", "mdip_sem": "sem_impressions"}
funnel_vars = list(funnel_mapping.values())
upper_channels = ["Newspaper", "Radio", "TV"]
control_columns = [col for col in raw_df.columns if "hldy_" in col]
data_df = raw_df[
["wk_strt_dt", target_column, *channel_mapping, *funnel_mapping, *control_columns]
].rename(columns=channel_mapping | funnel_mapping | {"wk_strt_dt": date_column})
data_df[date_column] = pd.to_datetime(data_df[date_column])
X = data_df.drop(columns=target_column)
y = data_df[target_column]
data_df.head()
| date | sales | Direct Mail | Insert | Newspaper | Radio | TV | Social Media | Online Display | sem_spend | ... | hldy_Memorial Day | hldy_Mother's Day | hldy_NYE | hldy_New Year's Day | hldy_Pre Thanksgiving | hldy_Presidents Day | hldy_Prime Day | hldy_Thanksgiving | hldy_Valentine's Day | hldy_Veterans Day | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2014-08-03 | 72051457.64 | 678410.26 | 129759.66 | 507635.16 | 254134.14 | 216725.59 | 0.0 | 61364.75 | 355954.51 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 1 | 2014-08-10 | 78794770.54 | 661097.71 | 70355.33 | 625719.37 | 106602.98 | 108919.26 | 0.0 | 132643.90 | 336306.72 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2 | 2014-08-17 | 70071185.56 | 1874338.43 | 107720.37 | 277532.61 | 59011.42 | 74032.19 | 0.0 | 76741.27 | 344787.96 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 3 | 2014-08-24 | 68642464.59 | 306583.37 | 141708.61 | 186028.61 | 178010.46 | 110995.66 | 0.0 | 193043.55 | 365495.11 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 4 | 2014-08-31 | 86190784.65 | 2366274.16 | 160691.30 | 706694.09 | 144097.56 | 225140.84 | 0.0 | 249612.89 | 401789.35 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5 rows × 33 columns
Exploratory Data Analysis#
The first iteration covers the full exploratory analysis of channels, controls, and sales. Here we look only at what is new: the SEM pair and what it implies for the funnel.
Let’s start with the spend share, to see how big this line item really is.
spend_share = (
pd.concat(
[X[channel_columns], X[["sem_spend"]].rename(columns={"sem_spend": "SEM"})],
axis=1,
)
.sum()
.sort_values()
)
spend_share = spend_share / spend_share.sum()
fig, ax = plt.subplots()
colors = ["C1" if name == "SEM" else "C0" for name in spend_share.index]
spend_share.plot.barh(color=colors, ax=ax)
ax.xaxis.set_major_formatter(mtick.PercentFormatter(xmax=1, decimals=0))
ax.set(title="Share of Total Media Spend (including SEM)", xlabel="spend share")
print(f"SEM share of the modeled media spend: {spend_share['SEM']:.1%}")
SEM share of the modeled media spend: 26.9%
SEM (orange) is \(26.9\%\) of the modeled media spend, second only to Direct Mail. It is not surprising that the stakeholders are curious about this channel. Still, adding it to the flat baseline model from the first iteration would have biased the upper funnel channels.
Next, let’s look at the two SEM series over time.
fig, ax = plt.subplots()
sns.lineplot(
data=data_df, x=date_column, y="sem_spend", color="C1", label="sem_spend", ax=ax
)
ax2 = ax.twinx()
sns.lineplot(
data=data_df,
x=date_column,
y="sem_impressions",
color="C2",
label="sem_impressions (right axis)",
ax=ax2,
)
ax2.grid(False)
ax2.set(ylabel="sem_impressions")
handles, labels = ax.get_legend_handles_labels()
handles2, labels2 = ax2.get_legend_handles_labels()
ax2.get_legend().remove()
ax.legend(handles + handles2, labels + labels2, loc="upper left")
ax.set(title="SEM Spend and SEM Impressions", ylabel="sem_spend");
The two series clearly move together. We can quantify this with the correlations against sales.
corr_sales = data_df[["sem_spend", "sem_impressions", target_column]].corr()[
target_column
]
corr_pair = data_df["sem_spend"].corr(data_df["sem_impressions"])
print(f"corr(sem_spend, sales) = {corr_sales['sem_spend']:.2f}")
print(f"corr(sem_impressions, sales) = {corr_sales['sem_impressions']:.2f}")
print(f"corr(sem_spend, sem_impressions) = {corr_pair:.2f}")
corr(sem_spend, sales) = 0.59
corr(sem_impressions, sales) = 0.59
corr(sem_spend, sem_impressions) = 0.79
SEM spend and impressions move together (correlation \(0.79\)) and both track sales (correlations \(0.59\)). Note, however, that this pattern is compatible with two very different stories: SEM drives sales, or SEM follows demand that something else created and sales follow the same demand.
Warning
Critical Assessment: correlation cannot separate “SEM causes sales” from “SEM harvests demand”. The practical question is causal: if we doubled SEM spend next quarter, would sales move? Only structure plus experiments can answer it, which is why this iteration adds both.
fig, axes = plt.subplots(
nrows=len(channel_columns),
ncols=1,
figsize=(12, 16),
sharex=True,
layout="constrained",
)
for ax, channel in zip(axes, channel_columns, strict=True):
sns.lineplot(data=data_df, x=date_column, y=channel, color="C0", ax=ax)
ax.set(title=channel, ylabel="spend")
fig.suptitle("Weekly Spend per Channel", fontsize=18, fontweight="bold");
Here is an important observation: Newspaper spend collapses to almost nothing in the recent period. The panel shows a steady decline from its early level to a weekly average of about \(33K\) over the last year, and to effectively zero over the last quarter.
Let’s keep this in mind. It is why the budget planning section later anchors on the last \(52\) weeks: last-quarter averages would freeze Newspaper at a level that no longer represents a real decision.
Model Specification#
The direct part of the model is the first iteration model unchanged: geometric adstock, logistic saturation, holiday controls, yearly seasonality, and spend-share priors on the saturation betas. See MMM End-to-End Case Study for the reasoning behind these priors.
The configuration helper below computes the spend shares from whatever data it receives, namely the full \(209\) weeks for the final model, and only the fold’s training window inside the cross-validation.
Tip
We pass the per-channel prior width as a labelled xr.DataArray so that the channel dimension is explicit instead of inferred from position.
def make_model_config(X: pd.DataFrame) -> dict:
"""Build the first iteration model configuration with spend-share prior widths from X."""
spend_shares = (
X[channel_columns].sum().reindex(channel_columns)
/ X[channel_columns].sum().sum()
).to_numpy()
return {
"intercept": Prior("Normal", mu=0.2, sigma=0.05),
"saturation_beta": Prior(
"HalfNormal",
sigma=xr.DataArray(spend_shares, dims="channel"),
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)
),
}
The Funnel Effect#
The funnel enters the model as a custom MuEffect: a component that adds a term to the sales mean and can declare its own likelihoods. We subclass DataVarMuEffect, which registers sem_spend and sem_impressions as pm.Data containers in the model, so prediction calls can swap in new values without rebuilding the model.
Four design points are worth explaining, because each one prevents a silent failure:
The effect reads the model’s own channel data (
mmm.channel_data_scaled) and selects the upper channels by name. This is what makes budget optimization work later: an intervention on the channel budgets flows through the demand pool automatically. An effect that reads its own copy of the spend would be invisible to the optimizer.The upper channel subset gets its own named dimension (
channel_upper), so its adstock and saturation priors cannot collide with the seven-channelchanneldimension.incrementality_specdeclares the extra carryover the second adstock adds. This is how the ROAS machinery and the optimizer window arithmetic learn that the mediated path reaches further into the future than the direct one.funnel_baselinecompetes with the sales intercept, as a constant shift in latent demand produces a near-constant mediated contribution. The tight intercept prior inherited from the first iteration resolves that ridge in this fit. With a wider intercept prior, this confound is the first thing you should check.
upper_dim = "channel_upper"
class FunnelEffect(DataVarMuEffect):
"""Upper funnel spend -> latent demand -> SEM-harvested mediator -> sales.
The measurement coefficient on the mediator is fixed at one (times the
impressions scale). This anchors the mediator scale and removes the
likelihood-invariant ridge a free coefficient would create.
"""
upper_channels: list[str]
sem_scale: float
impressions_scale: float
upper_transform: InstanceOf[MediaTransformation]
demand_transform: InstanceOf[MediaTransformation]
model_config = {"arbitrary_types_allowed": True}
def to_dict(self) -> dict:
"""Serialize the effect."""
return {
"data_vars": self.data_vars,
"prefix": self.prefix,
"upper_channels": self.upper_channels,
"sem_scale": self.sem_scale,
"impressions_scale": self.impressions_scale,
"upper_transform": self.upper_transform.to_dict(),
"demand_transform": self.demand_transform.to_dict(),
}
def incrementality_spec(self) -> IncrementalitySpec:
"""Opt in to incrementality, declaring the second adstock's reach."""
return IncrementalitySpec(
additional_carryover_lags=self.demand_transform.adstock.l_max
)
def create_data(self, mmm) -> None:
"""Register the funnel data variables and the upper channel coordinate."""
mmm.model.add_coord(upper_dim, self.upper_channels)
super().create_data(mmm)
def create_effect(self, mmm):
"""Build the mediator equation and the term it adds to the sales mean."""
model = mmm.model
channel_index = [mmm.channel_columns.index(c) for c in self.upper_channels]
upper = mmm.channel_data_scaled.isel(channel=channel_index).rename(
{"channel": upper_dim}
)
upper_on_demand = self.upper_transform(upper, dim="date").sum(dim=upper_dim)
baseline = pmd.HalfNormal(f"{self.prefix}_baseline", sigma=1.0)
demand = pmd.Deterministic(f"{self.prefix}_demand", baseline + upper_on_demand)
lam = pmd.HalfNormal(f"{self.prefix}_lambda", sigma=0.5)
mediator = pmd.Deterministic(
f"{self.prefix}_mediator",
demand + lam * model["sem_spend"] / self.sem_scale,
)
pmd.TruncatedNormal(
f"{self.prefix}_sem_impressions_likelihood",
mu=self.impressions_scale * mediator,
sigma=pmd.HalfNormal(
f"{self.prefix}_sigma_s", sigma=0.2 * self.impressions_scale
),
lower=0.0,
observed=model["sem_impressions"],
)
return pmd.Deterministic(
f"{self.prefix}_effect_contribution",
self.demand_transform(mediator, dim="date"),
)
l_max = 6 # direct path adstock, as in the first iteration
# The mediator adstock length is a modeling declaration, not an estimate: with weekly
# data we assume harvested demand converts within about six weeks, mirroring the direct
# path. It is deliberate and conservative, because effective_carryover_lags uses it to
# widen the optimizer window.
l_max_lf = 6
def make_funnel_effect(sem_scale: float, impressions_scale: float) -> FunnelEffect:
"""Build the funnel effect with its two media transformations."""
return FunnelEffect(
data_vars=funnel_vars,
prefix="funnel",
upper_channels=upper_channels,
sem_scale=sem_scale,
impressions_scale=impressions_scale,
upper_transform=MediaTransformation(
adstock=GeometricAdstock(
l_max=l_max,
prefix="adstock_uf",
priors={"alpha": Prior("Beta", alpha=2, beta=3, dims=upper_dim)},
),
saturation=LogisticSaturation(
prefix="sat_uf",
priors={
"lam": Prior("Gamma", mu=2.0, sigma=1.0, dims=upper_dim),
"beta": Prior("HalfNormal", sigma=1.0, dims=upper_dim),
},
),
adstock_first=True,
dims=(upper_dim,),
),
demand_transform=MediaTransformation(
adstock=GeometricAdstock(
l_max=l_max_lf,
prefix="adstock_lf",
priors={"alpha": Prior("Beta", alpha=2, beta=3)},
),
saturation=LogisticSaturation(
prefix="sat_lf",
priors={
"lam": Prior("Gamma", mu=1.5, sigma=0.5),
"beta": Prior("HalfNormal", sigma=1.0),
},
),
adstock_first=True,
dims=(),
),
)
Data Plumbing#
There are two important points when building funnel models:
First, the model must be built from an xr.Dataset, because the funnel columns (sem_spend, sem_impressions) only survive the data conversion as extra variables of a Dataset. Fitting still uses the flat DataFrame.
Second, every prediction call converts DataFrames without the extra variables by default, which would silently freeze the funnel data at its training values. We subclass the model once to add them back. With this override, all posterior predictive workflows in this notebook (cross-validation, counterfactuals, uplift) can work with plain DataFrames, and a frame that lacks the funnel columns is refused.
Note
One consequence of the first point: the stored fit_data carries no funnel columns, so save and load round-trips do not reconstruct this model.
Both mechanics are candidates for the library itself. The model knows which data variables its effects declare, so it could carry them through build, prediction, cross-validation, and save and load without a subclass. That proposal is tracked in #2945.
class FunnelMMM(MMM):
"""MMM that carries the funnel data variables through DataFrame conversion."""
def _posterior_predictive_data_transformation(
self, X, y=None, include_last_observations: bool = False
):
if isinstance(X, pd.DataFrame):
missing = [var for var in funnel_vars if var not in X.columns]
if missing:
raise ValueError(
f"funnel columns missing from prediction frame: {missing}"
)
X = to_mmm_dataset(
X,
y,
date_column=self.date_column,
dims=self.dims,
channel_columns=self.channel_columns,
control_columns=self.control_columns,
target_column=self.target_column,
extra_vars=funnel_vars,
)
y = None
return super()._posterior_predictive_data_transformation(
X, y=y, include_last_observations=include_last_observations
)
With the plumbing settled, here are the two helpers that build the dataset and the model itself.
def make_funnel_dataset(X: pd.DataFrame, y: pd.Series | None = None) -> xr.Dataset:
"""Convert the flat frame to the Dataset the funnel model builds from."""
return to_mmm_dataset(
X,
y,
date_column=date_column,
channel_columns=channel_columns,
control_columns=[col for col in X.columns if "hldy_" in col],
target_column=target_column,
extra_vars=funnel_vars,
)
def make_mmm(X: pd.DataFrame) -> FunnelMMM:
"""Build the funnel MMM with scales and prior widths computed from X."""
mmm = FunnelMMM(
model_config=make_model_config(X),
target_column=target_column,
date_column=date_column,
adstock=GeometricAdstock(l_max=l_max),
saturation=LogisticSaturation(),
channel_columns=channel_columns,
control_columns=[col for col in X.columns if "hldy_" in col],
yearly_seasonality=5,
)
mmm.add_mu_effect(
make_funnel_effect(
sem_scale=float(X["sem_spend"].max()),
impressions_scale=float(X["sem_impressions"].max()),
)
)
return mmm
ROAS Calibration from Experiments#
The first iteration flagged its widest uncertainty exactly where the budget decisions were largest, and recommended experiments as the next step. So for this iteration, let’s imagine the measurement team ran two experiments in the first half of 2017:
A conversion lift study on
Online Display, the channel with the highest estimated ROAS and the widest uncertainty band. The study reported a ROAS between \(20\) and \(36\).A regional holdout test on
Direct Mail, the largest and most contested budget line. The test reported a ROAS between \(2\) and \(8\), confirming the low estimate from the model.
On a real application these numbers come from geo holdouts or conversion lift studies. See Lift Test Calibration for how to design them, and Calibrating an MMM with ROAS Estimates for the calibration API we use below.
Before we use them, let’s do a quick plausibility check on the scale. Total sales are about \(22.6B\) against \(355M\) of direct channel spend. So if media drives roughly a third of sales, the blended direct ROAS is around \(20\) (\(22.58\text{B} / 3 / 355.4\text{M} \approx 21\), printed below). Both reported ranges sit comfortably on that scale.
We treat each reported range as a \(95\%\) interval and convert it to the value and standard deviation the calibration API expects.
Note
This is the one place the notebook uses a \(95\%\) figure. It is an input reported by the measurement team, not a model interval. Every model interval from here on is a \(94\%\) HDI.
def roas_range_to_calibration(lower: float, upper: float) -> dict[str, float]:
"""Convert a reported 95% ROAS interval to a Normal(value, sigma) target."""
return {"roas": (lower + upper) / 2, "sigma": (upper - lower) / (2 * 1.96)}
roas_experiments = {
"Online Display": roas_range_to_calibration(20.0, 36.0),
"Direct Mail": roas_range_to_calibration(2.0, 8.0),
}
roas_calibration_df = pd.DataFrame(
[{"channel": channel, **target} for channel, target in roas_experiments.items()]
)
total_sales = float(y.sum())
total_direct_spend = float(X[channel_columns].to_numpy().sum())
print(
f"total sales {total_sales / 1e9:.2f}B, direct channel spend {total_direct_spend / 1e6:.1f}M, "
f"blended ROAS if media drives a third of sales: {total_sales / 3 / total_direct_spend:.1f}"
)
roas_calibration_df
total sales 22.58B, direct channel spend 355.4M, blended ROAS if media drives a third of sales: 21.2
| channel | roas | sigma | |
|---|---|---|---|
| 0 | Online Display | 28.0 | 4.081633 |
| 1 | Direct Mail | 5.0 | 1.530612 |
Warning
The calibration constrains the model’s direct channel contribution (add_cost_per_target_calibration() targets channel_contribution_original_scale). That is exactly right here because both calibrated channels are direct-only in our graph.
Calibrating an upper funnel channel this way would pin its direct path while the experiment measured its total effect, direct plus mediated, so the calibration target would have to span both paths. Always match the calibration target to the estimand of the experiment! Calibrating upper funnel channels is out of scope here and is the subject of a separate notebook (see next step 2 at the end).
Time-Slice Cross-Validation#
The first iteration validated on a single \(26\)-week holdout and flagged the limits of that design: one split, one verdict, and a test window that happened to be easier to predict than the training data. As promised in its next steps, this iteration uses TimeSliceCrossValidator instead. See Time-Slice-Cross-Validation and Parameter Stability for a deeper look at the technique.
Here is the fold design. The first training window is \(157\) weeks (about three years), each fold forecasts the next \(13\) weeks (one quarter), and the window grows by \(13\) weeks per fold. With \(209\) weeks of data that gives \(\lfloor (209 - 157 - 13) / 13 \rfloor + 1 = 4\) folds, so we test the model on the last four quarters of the data.
The validator rebuilds and refits the model on each fold. Because our model needs a Dataset at build time and a per-fold ROAS calibration, we hand the validator a small builder object instead of a model instance. Anything with a build_model(X, y) method works: the validator deep-copies it per fold and uses what it returns.
Note
Leakage control: everything data-derived is computed inside the builder from the fold’s own training window, including the spend-share prior widths, the SEM scaling constants, and the calibration spend table. The experiment ROAS values themselves date from the first half of 2017, before the first fold’s training window ends, so every fold is allowed to see them.
class CalibratedFunnelMMMBuilder:
"""Fold-level factory for the time-slice cross-validator."""
def __init__(self, calibration_data: pd.DataFrame):
self.calibration_data = calibration_data.copy()
def build_model(self, X: pd.DataFrame, y: pd.Series) -> FunnelMMM:
"""Build a fresh, ROAS-calibrated funnel model from a training fold."""
mmm = make_mmm(X)
mmm.build_model(X=make_funnel_dataset(X, y))
mmm.add_original_scale_contribution_variable(
var=["channel_contribution", "funnel_effect_contribution", "y"]
)
mmm.add_cost_per_target_calibration(
data=X,
calibration_data=self.calibration_data,
name_prefix="roas_calibration",
target_column="roas",
target_per_cost=True,
)
return mmm
Now we set up the cross-validator itself and confirm we get the four folds we expect.
cv = TimeSliceCrossValidator(
n_init=157,
forecast_horizon=13,
date_column=date_column,
step_size=13,
)
cv.plot_suite = "new"
n_splits = cv.get_n_splits(X, y)
print(f"Number of folds: {n_splits}")
if n_splits != 4:
raise RuntimeError(f"expected 4 folds, got {n_splits}")
Number of folds: 4
We are now ready to run the cross-validation. This refits the model four times, so it takes a while.
cv_results = cv.run(
X,
y,
mmm=CalibratedFunnelMMMBuilder(calibration_data=roas_calibration_df),
sampler_config={
"tune": 1_000,
"draws": 500,
"chains": 4,
"nuts_sampler": "nutpie",
"target_accept": 0.9,
"random_seed": seed,
"progressbar": False,
},
model_names=[f"Fold {i}" for i in range(n_splits)],
)
NUTS[nutpie]: [y_sigma, adstock_lf_alpha, funnel_lambda, adstock_uf_alpha, sat_uf_lam, sat_uf_beta, funnel_baseline, sat_lf_lam, sat_lf_beta, gamma_fourier, gamma_control, adstock_alpha, saturation_lam, saturation_beta, intercept_contribution, funnel_sigma_s]
Sampling: [funnel_sem_impressions_likelihood, roas_calibration, y]
NUTS[nutpie]: [y_sigma, adstock_lf_alpha, funnel_lambda, adstock_uf_alpha, sat_uf_lam, sat_uf_beta, funnel_baseline, sat_lf_lam, sat_lf_beta, gamma_fourier, gamma_control, adstock_alpha, saturation_lam, saturation_beta, intercept_contribution, funnel_sigma_s]
Sampling: [funnel_sem_impressions_likelihood, roas_calibration, y]
NUTS[nutpie]: [y_sigma, adstock_lf_alpha, funnel_lambda, adstock_uf_alpha, sat_uf_lam, sat_uf_beta, funnel_baseline, sat_lf_lam, sat_lf_beta, gamma_fourier, gamma_control, adstock_alpha, saturation_lam, saturation_beta, intercept_contribution, funnel_sigma_s]
Sampling: [funnel_sem_impressions_likelihood, roas_calibration, y]
NUTS[nutpie]: [y_sigma, adstock_lf_alpha, funnel_lambda, adstock_uf_alpha, sat_uf_lam, sat_uf_beta, funnel_baseline, sat_lf_lam, sat_lf_beta, gamma_fourier, gamma_control, adstock_alpha, saturation_lam, saturation_beta, intercept_contribution, funnel_sigma_s]
Sampling: [funnel_sem_impressions_likelihood, roas_calibration, y]
Before looking at any result, we check that the sampler was happy across all folds.
n_divergences_cv = int(
cv_results["sample_stats"].to_dataset()["diverging"].sum().item()
)
print(f"Total divergences across the {n_splits} folds: {n_divergences_cv}")
if n_divergences_cv != 0:
warnings.warn(f"{n_divergences_cv} divergences across the CV folds", stacklevel=1)
Total divergences across the 4 folds: 0
Out-of-Sample Predictions per Fold#
Each panel shows one fold: the \(94\%\) HDI of the posterior predictive over the training window, the same band over the held-out quarter (the forecast), and the observed sales line.
cv.plot.predictions(hdi_prob=HDI_PROB, figsize=(15, 18));
The forecast bands cover the observed sales in all four held-out quarters, which is good! The hardest one is Fold 1, whose test window contains the 2017 holiday peak. The model covers the spike, but with visibly wider uncertainty, and this is the fold with the highest test score in the CRPS table below.
Parameter Stability#
Here we are interested in the following question: do the channel effects and the funnel strength depend on which window we fit? Stable medians across folds mean the conclusions are not an artifact of one period.
For the two calibrated channels, part of the stability is by construction, since every fold sees the same two experiments and their saturation betas are anchored. So the interesting rows are the uncalibrated ones.
In the forests below there is one row per fold and channel, where the thick bar is the \(50\%\) HDI, the thin bar the \(94\%\) HDI, and the dot the point estimate.
cv.plot.param_stability(
var_names=["saturation_beta"],
figsize=(12, 9),
figure_kwargs={"layout": "constrained"},
);
We do the same for the adstock and saturation parameters.
cv.plot.param_stability(
var_names=["adstock_alpha", "saturation_lam"],
figsize=(12, 14),
figure_kwargs={"layout": "constrained"},
);
The funnel strength funnel_lambda is the parameter the whole SEM analysis relies on, as it measures how strongly an extra unit of SEM budget lifts the mediator. We should take a closer look.
It is a scalar, so the channel-oriented stability plot and its tabular summary do not apply. A plain forest over the folds shows the same thing.
pc = azp.plot_forest(
cv_results["posterior"].to_dataset()[["funnel_lambda"]],
combined=True,
figure_kwargs={"figsize": (10, 4)},
)
fig = pc.viz["/"]["figure"].values.item()
fig.suptitle(
"Stability of the Funnel Strength Across Folds", fontsize=16, fontweight="bold"
);
def summary_table(da: xr.DataArray) -> pd.DataFrame:
"""Posterior mean and 94% HDI over the remaining dims, as a DataFrame."""
hdi = az.hdi(da, prob=HDI_PROB)
return xr.Dataset(
{
"mean": da.mean(("chain", "draw")),
"hdi_low": hdi.sel(ci_bound="lower", drop=True),
"hdi_high": hdi.sel(ci_bound="upper", drop=True),
}
).to_dataframe()
funnel_lambda_folds = cv_results["posterior"].to_dataset()["funnel_lambda"]
funnel_lambda_by_fold = summary_table(funnel_lambda_folds)
funnel_lambda_by_fold.round(3)
| mean | hdi_low | hdi_high | |
|---|---|---|---|
| cv | |||
| Fold 0 | 0.701 | 0.602 | 0.802 |
| Fold 1 | 0.696 | 0.598 | 0.798 |
| Fold 2 | 0.688 | 0.600 | 0.761 |
| Fold 3 | 0.683 | 0.611 | 0.764 |
The media parameters are stable! The per-fold posteriors overlap heavily for every channel, with no jumps as the training window grows quarter by quarter. The two calibrated channels show the tightest saturation posteriors, which is partly by construction as noted above.
The funnel strength is stable too. The four fold means of funnel_lambda sit between \(0.68\) and \(0.70\) and their \(94\%\) HDIs overlap almost entirely (each spans roughly \(0.60\) to \(0.80\)). Each additional quarter of data tells us the same story about how strongly SEM budget moves the mediator.
Forecast Quality#
The continuous ranked probability score (CRPS) scores the whole forecast distribution against the observed value, where lower is better. We compare training and test scores per fold. The signature to watch for is a test score drifting up while the training score drifts down. That is what overfitting looks like, and it is exactly what the single split in the first iteration could not reveal.
cv.plot.crps(combine_splits=True, figsize=(10, 6));
The same scores in tabular form:
cv.summary.crps()
| split | cv | mean_crps | |
|---|---|---|---|
| 0 | train | Fold 0 | 1.741123e+07 |
| 1 | train | Fold 1 | 1.718090e+07 |
| 2 | train | Fold 2 | 1.788575e+07 |
| 3 | train | Fold 3 | 1.763065e+07 |
| 4 | test | Fold 0 | 1.532885e+07 |
| 5 | test | Fold 1 | 2.923965e+07 |
| 6 | test | Fold 2 | 1.511433e+07 |
| 7 | test | Fold 3 | 1.736371e+07 |
Training CRPS is flat across folds at \(17\) to \(18\) million. Two test scores sit clearly below the training level (Folds 0 and 2, about \(15\) million) and Fold 3 sits right at it (its test score lands within a few tenths of a million of its training score).
Fold 1 is the exception at about \(29\) million. Recall its test window is the 2017 holiday quarter, the hardest \(13\) weeks of the year to forecast.
Most importantly, there is no widening gap between train and test as the window grows. This is precisely the overfitting signature this check exists to catch.
Tip
With stability established, we can now go for what we actually want: a final model fit on all the data. We discard the cross-validation models, as their job was to certify the specification, not to make predictions.
Final Model#
The final model uses the same builder logic on the full \(209\) weeks. The API fixes the required order: build the model, register the original-scale contribution variables, add the ROAS calibration, then fit.
mmm = make_mmm(X)
mmm.build_model(X=make_funnel_dataset(X, y))
mmm.add_original_scale_contribution_variable(
var=[
"channel_contribution",
"control_contribution",
"intercept_contribution",
"yearly_seasonality_contribution",
"funnel_effect_contribution",
"y",
]
)
mmm.add_cost_per_target_calibration(
data=X,
calibration_data=roas_calibration_df,
name_prefix="roas_calibration",
target_column="roas",
target_per_cost=True,
)
pm.model_to_graphviz(mmm.model)
The graph shows what this iteration actually fits: three data sources at once. The sales likelihood, the SEM impressions likelihood inside the funnel effect, and the two experiment observations under roas_calibration.
Prior Predictive Check#
Before fitting, we confirm the priors put sales in a plausible range: weakly informative, roughly the right order of magnitude.
mmm.sample_prior_predictive(X, y, samples=2_000, random_seed=seed)
fig, ax = plt.subplots()
prior_y = mmm.idata["prior"]["y_original_scale"]
dates = X[date_column]
for prob, alpha in [(HDI_PROB, 0.2), (0.5, 0.4)]:
hdi = az.hdi(prior_y, prob=prob)
ax.fill_between(
dates,
hdi.sel(ci_bound="lower"),
hdi.sel(ci_bound="upper"),
color="C0",
alpha=alpha,
label=f"{prob:.0%} HDI",
)
ax.plot(dates, y, color="black", label="observed sales")
ax.legend(loc="upper right")
ax.set(title="Prior Predictive Check", xlabel="date", ylabel="sales");
Sampling: [adstock_alpha, adstock_lf_alpha, adstock_uf_alpha, funnel_baseline, funnel_lambda, funnel_sem_impressions_likelihood, funnel_sigma_s, gamma_control, gamma_fourier, intercept_contribution, roas_calibration, sat_lf_beta, sat_lf_lam, sat_uf_beta, sat_uf_lam, saturation_beta, saturation_lam, y, y_sigma]
The priors are not very informative: the bands reach about an order of magnitude above the observed sales level, and the observed series sits well inside them. This is what we want at this stage.
Model Fitting#
We are now ready to fit the model on all \(209\) weeks.
mmm.fit(
X=X,
y=y,
target_accept=0.95,
chains=4,
cores=4,
draws=2_000,
tune=2_000,
nuts_sampler="nutpie",
random_seed=rng,
)
mmm.sample_posterior_predictive(X, random_seed=rng)
mmm.plot_suite = "new"
NUTS[nutpie]: [y_sigma, adstock_lf_alpha, funnel_lambda, adstock_uf_alpha, sat_uf_lam, sat_uf_beta, funnel_baseline, sat_lf_lam, sat_lf_beta, gamma_fourier, gamma_control, adstock_alpha, saturation_lam, saturation_beta, intercept_contribution, funnel_sigma_s]
Sampling: [funnel_sem_impressions_likelihood, roas_calibration, y]
Model Diagnostics#
Before looking into the results, let’s check some diagnostics. We hold the model to the same bar as the first iteration: zero divergences, all \(\hat{R}\) below \(1.01\), and clean traces. The funnel block simply adds a handful of parameters to the checklist.
We start with arviz.diagnose(), which runs the standard sampler checks in one call (divergences, energy fraction of missing information, effective sample size, and \(\hat{R}\) against the thresholds it prints), and then deep dive into each item as before.
n_divergences = int(mmm.idata["sample_stats"]["diverging"].sum().item())
print(f"Number of divergences: {n_divergences}")
if n_divergences != 0:
warnings.warn(f"the sampler reported {n_divergences} divergences", stacklevel=1)
Number of divergences: 0
diagnostic_var_names = [
"intercept_contribution",
"y_sigma",
"saturation_beta",
"saturation_lam",
"adstock_alpha",
"gamma_control",
"gamma_fourier",
"funnel_baseline",
"funnel_lambda",
"funnel_sigma_s",
"adstock_uf_alpha",
"sat_uf_lam",
"sat_uf_beta",
"adstock_lf_alpha",
"sat_lf_lam",
"sat_lf_beta",
]
sampler_ok = not az.diagnose(mmm.idata, var_names=diagnostic_var_names)
r_hat_max = max(
float(values.max())
for values in az.rhat(mmm.idata, var_names=diagnostic_var_names).data_vars.values()
)
print(f"max r_hat (unrounded): {r_hat_max:.4f}")
print(f"sampler checks passed: {sampler_ok}")
Divergences
No divergent transitions found.
E-BFMI
E-BFMI satisfactory for all chains.
ESS
Effective sample size satisfactory for all parameters.
R-hat
R-hat values satisfactory for all parameters.
Processing complete, no problems detected.
max r_hat (unrounded): 1.0052
sampler checks passed: True
default_max_subplots = az.rcParams["plot.max_subplots"]
az.rcParams["plot.max_subplots"] = 100
azp.plot_trace_dist(
mmm.fit_result,
var_names=[
"saturation_beta",
"adstock_alpha",
"funnel_lambda",
"sat_uf_beta",
"sat_lf_beta",
],
figure_kwargs={"figsize": (12, 15)},
)
az.rcParams["plot.max_subplots"] = default_max_subplots
We see a good mixing in the traces for the media and funnel blocks alike. This is consistent with the summary above: no divergences, satisfactory energy fraction of missing information, effective sample size, and \(\hat{R}\) for every parameter (the largest unrounded \(\hat{R}\) is \(1.005\)).
The funnel block deserves a closer look than the trace alone, because its parameters carry very different weights of evidence.
az.summary(
mmm.idata, var_names=["funnel_lambda", "sat_lf_beta", "sat_lf_lam"], kind="stats"
)
| mean | sd | hdi94_lb | hdi94_ub | |
|---|---|---|---|---|
| funnel_lambda | 0.67 | 0.042 | 0.59 | 0.75 |
| sat_lf_beta | 0.28 | 0.2 | 0.021 | 0.74 |
| sat_lf_lam | 1.4 | 0.47 | 0.62 | 2.4 |
How to read and interpret the funnel parameters?
funnel_lambda is tight (\(0.67\), \(94\%\) HDI \(0.59\) to \(0.75\)) because the impressions data speak directly to it. It sits clearly away from zero, but keep in mind that the impressions likelihood identifies it, and search spend buying search impressions is close to tautological.
The parameter that carries the mediated path into sales is sat_lf_beta, with sat_lf_lam shaping the curve. A weak sat_lf_beta would mean the demand pool exists but does not convert. Here it is positive-leaning but wide (mean \(0.28\), \(94\%\) HDI \(0.02\) to \(0.74\), its lower edge at zero).
In other words, the sales data support conversion of the demand pool without pinning its strength. This is exactly why the mediated shares later in the notebook come out positive but imprecise.
Posterior Predictive Checks#
As our model and posterior samples are looking good, we can now look into the posterior predictive checks. This model has two likelihoods, so we get two of them.
Sales#
fig, ax = plt.subplots()
post_y = mmm.idata["posterior_predictive"]["y_original_scale"]
for prob, alpha in [(HDI_PROB, 0.2), (0.5, 0.4)]:
hdi = az.hdi(post_y, prob=prob)
ax.fill_between(
dates,
hdi.sel(ci_bound="lower"),
hdi.sel(ci_bound="upper"),
color="C0",
alpha=alpha,
label=f"{prob:.0%} HDI",
)
ax.plot(dates, y, color="black", label="observed sales")
ax.legend(loc="upper right")
ax.set(title="Posterior Predictive Check", xlabel="date", ylabel="sales");
The posterior predictive sales look quite reasonable. The model tracks the weekly level and captures all four holiday peaks, and the observed line stays inside the bands throughout.
SEM Impressions#
Here is a check that does not even exist in a flat model. The funnel model claims the latent mediator explains SEM impressions, so we can hold the model’s implied impressions against the observed series.
Business Value: this is a check the search team can audit against their own dashboards. If the model’s demand story disagreed with the impressions they see, the funnel would be wrong in a way a sales fit could never reveal.
funnel_effect = next(
effect for effect in mmm.mu_effects if getattr(effect, "prefix", None) == "funnel"
)
impressions_scale = funnel_effect.impressions_scale
implied_impressions = impressions_scale * mmm.idata["posterior"]["funnel_mediator"]
fig, ax = plt.subplots()
for prob, alpha in [(HDI_PROB, 0.2), (0.5, 0.4)]:
hdi = az.hdi(implied_impressions, prob=prob)
ax.fill_between(
dates,
hdi.sel(ci_bound="lower"),
hdi.sel(ci_bound="upper"),
color="C2",
alpha=alpha,
label=f"{prob:.0%} HDI",
)
ax.plot(dates, X["sem_impressions"], color="black", label="observed sem_impressions")
ax.legend(loc="upper right")
ax.set(
title="Measurement Check: Implied vs Observed SEM Impressions",
xlabel="date",
ylabel="impressions",
);
The implied impressions reproduce the level of the observed series and every seasonal demand peak, including their heights.
Note that the narrow bands show the implied mean of the measurement and not its noise, so the week-to-week wiggles of the observed line sit outside them by design. For an audit, the level and the peaks are what matter, and both match.
Now we add the noise back. This is the version you should use to decide how much to trust the measurement model.
pp_impressions = mmm.idata["posterior_predictive"]["funnel_sem_impressions_likelihood"]
observed_impressions = X["sem_impressions"].to_numpy()
impressions_coverage = {}
fig, ax = plt.subplots()
for prob, alpha in [(HDI_PROB, 0.2), (0.5, 0.4)]:
hdi = az.hdi(pp_impressions, prob=prob)
lower = hdi.sel(ci_bound="lower").to_numpy()
upper = hdi.sel(ci_bound="upper").to_numpy()
inside = (observed_impressions >= lower) & (observed_impressions <= upper)
impressions_coverage[prob] = float(inside.mean())
ax.fill_between(
dates, lower, upper, color="C2", alpha=alpha, label=f"{prob:.0%} HDI"
)
ax.plot(dates, observed_impressions, color="black", label="observed sem_impressions")
ax.legend(loc="upper right")
ax.set(
title="Posterior Predictive Check: SEM Impressions Including Measurement Noise",
xlabel="date",
ylabel="impressions",
)
n_weeks = len(observed_impressions)
print(
f"observed weeks inside the 94% HDI: {impressions_coverage[HDI_PROB]:.1%} "
f"({round(impressions_coverage[HDI_PROB] * n_weeks)} of {n_weeks}); "
f"inside the 50% HDI: {impressions_coverage[0.5]:.1%}"
)
observed weeks inside the 94% HDI: 97.6% (204 of 209); inside the 50% HDI: 43.5%
This second plot is the full posterior predictive of the measurement, with the model’s noise (funnel_sigma_s) added on top of the implied mean. The observed line should now sit inside each band about as often as the band level says.
At the wide level it does: \(204\) of the \(209\) weeks (\(97.6\%\)) fall inside the \(94\%\) HDI. That is slightly above the nominal level, because a single noise scale set largely by the four holiday peaks is generous for the quieter weeks.
The \(50\%\) HDI is a different story, covering only \(43.5\%\) of the weeks. The plot shows why. In the first year the observed line runs along the lower edge of the dark band, and in the last year along its upper edge. This is a slow drift in the level of search demand that a constant funnel_baseline cannot follow, and it is why we propose a time-varying demand baseline in the next steps.
So back to our audit question: does the model’s demand story match what the search team sees? The answer is yes on level, peaks, and tail coverage, with the slow drift as the documented gap.
Tip
These coverage numbers are also a statement the search team can hold the model to. If a future quarter falls outside the \(94\%\) HDI, we revisit the funnel story before we touch the budget.
Note
A skeptical reader may ask how this differs from adding sem_impressions as a channel, the mistake we warned about at the beginning. The difference is causal: the spend counterfactuals below move the latent demand, and with it the mediator, without re-conditioning on the observed impressions. A channel model would hold the impressions fixed and block the upper funnel path.
Media Deep Dive#
Now that we are happy with the model fit, we can look into the media and individual channel contributions.
Contribution Decomposition#
We start by plotting each component contribution over time.
mmm.plot.decomposition.contributions_over_time(
include=["channels", "controls", "baseline", "seasonality"],
hdi_prob=HDI_PROB,
figsize=(12, 7),
);
Next, we look deeper into the components with the waterfall chart.
mmm.plot.decomposition.waterfall(figsize=(15, 12));
The waterfall has a blind spot we need to handle, so we also compute the shares ourselves.
def contribution_shares(model: MMM, components: dict[str, str]) -> pd.Series:
"""Posterior-mean share of each original-scale component in observed total sales."""
posterior = model.idata["posterior"]
totals = {}
for label, name in components.items():
da = posterior[f"{name}_original_scale"]
total = da.sum([d for d in da.dims if d not in ("chain", "draw")]).mean()
# a constant intercept is stored once, not once per week
weeks = 1 if "date" in da.dims else len(y)
totals[label] = float(total) * weeks
return pd.Series(totals, name="share of observed sales") / float(y.sum())
funnel_components = {
"baseline (intercept)": "intercept_contribution",
"direct channels": "channel_contribution",
"mediated SEM (funnel)": "funnel_effect_contribution",
"controls": "control_contribution",
"seasonality": "yearly_seasonality_contribution",
}
funnel_shares = contribution_shares(mmm, funnel_components)
waterfall_denominator = funnel_shares.drop("mediated SEM (funnel)").sum()
print(
f"baseline share as the waterfall computes it (funnel term excluded): "
f"{funnel_shares['baseline (intercept)'] / waterfall_denominator:.1%}"
)
funnel_shares.round(3).to_frame()
baseline share as the waterfall computes it (funnel term excluded): 66.3%
| share of observed sales | |
|---|---|
| baseline (intercept) | 0.561 |
| direct channels | 0.289 |
| mediated SEM (funnel) | 0.147 |
| controls | -0.004 |
| seasonality | -0.000 |
Warning
The waterfall knows a fixed set of components and does not show the funnel contribution. The mediated SEM mass is missing from this chart, so the media share it displays understates the true media-driven share. The counterfactual decomposition below is the right tool for the funnel part.
The waterfall attributes \(66\%\) of sales to the baseline (printed above from the same original-scale contributions the chart uses), against roughly \(60\%\) in the first iteration’s prose. That figure and the \(56\%\) in the table below are different quantities: the waterfall divides the intercept by the components it knows about, which exclude the funnel term, while the table divides every component by observed sales. Careful though, the two charts are not comparable line by line! The mediated SEM contribution is missing from this one entirely, so part of what the first iteration would have counted as media mass is invisible here.
The table above puts every component on the same footing, observed sales. The direct channels carry \(28.9\%\) and the mediated SEM path another \(14.7\%\). The comparison section near the end of the notebook refits the first iteration’s specification so that we can read the two decompositions side by side.
Direct and Mediated Contributions#
Let’s get to the stakeholder question behind the funnel: how much of the TV, Radio, and Newspaper effect arrives indirectly through SEM?
The model answers it with a counterfactual. We set one upper channel’s spend to zero, push the data through the posterior, and read how much the mediated contribution falls. That difference is the channel’s indirect effect.
uplift_var_names = [
"y_original_scale",
"channel_contribution_original_scale",
"funnel_effect_contribution_original_scale",
]
def sample_plan(frame: pd.DataFrame) -> xr.Dataset:
"""Posterior predictive pass for a spend plan (DataFrame with funnel columns)."""
return mmm.sample_posterior_predictive(
frame,
var_names=uplift_var_names,
combined=False,
extend_idata=False,
progressbar=False,
random_seed=seed,
)
pp_base = sample_plan(X)
indirect_contributions = {}
direct_deltas = {}
for channel in upper_channels:
X_zero = X.copy()
X_zero[channel] = 0.0
pp_zero = sample_plan(X_zero)
indirect = (
pp_base["funnel_effect_contribution_original_scale"]
- pp_zero["funnel_effect_contribution_original_scale"]
)
if float(indirect.mean()) <= 0:
warnings.warn(
f"zeroing {channel} did not reduce the mediated contribution", stacklevel=1
)
indirect_contributions[channel] = indirect
direct_deltas[channel] = pp_base["channel_contribution_original_scale"].sel(
channel=channel
).sum("date") - pp_zero["channel_contribution_original_scale"].sel(
channel=channel
).sum("date")
indirect_da = xr.concat(
list(indirect_contributions.values()),
dim=pd.Index(upper_channels, name="channel"),
)
Sampling: [y]
Sampling: [y]
Sampling: [y]
Sampling: [y]
We stack the two parts for each upper funnel channel.
fig, axes = plt.subplots(
nrows=len(upper_channels),
ncols=1,
figsize=(12, 10),
sharex=True,
layout="constrained",
)
for ax, channel in zip(axes, upper_channels, strict=True):
direct_mean = (
pp_base["channel_contribution_original_scale"]
.sel(channel=channel)
.mean(("chain", "draw"))
)
indirect_mean = indirect_da.sel(channel=channel).mean(("chain", "draw"))
ax.stackplot(
dates,
direct_mean,
indirect_mean,
labels=["direct", "mediated through SEM"],
colors=["C0", "C1"],
alpha=0.8,
)
ax.legend(loc="upper right")
ax.set(title=channel, ylabel="contribution")
fig.suptitle(
"Direct and Mediated Contributions (posterior means)",
fontsize=18,
fontweight="bold",
);
For each upper funnel channel, the blue band is what the channel earns on its own through its direct path. The orange band on top is what it earns through search, the demand it created that SEM later harvested. Both bands are posterior means.
Here are two important remarks:
First, the orange band peaks where the blue band peaks (same spend, shifted by the mediator’s own adstock). The search team’s holiday quarters are partly bought by the brand team’s holiday campaigns!
Second, for Newspaper both bands vanish in the recent period. Recall its weekly spend fell to about \(33K\) over the last year and to effectively zero over the last quarter. A channel that once fed the demand pool has simply stopped doing so.
Business Value: this is the chart for the attribution dispute between the brand and the search teams. Search revenue is not all search’s doing, and the orange band is the size of the brand team’s claim on it.
We can turn that into a share per channel: how much of each channel’s total effect arrives through search?
direct_delta_da = xr.concat(
[direct_deltas[channel] for channel in upper_channels],
dim=pd.Index(upper_channels, name="channel"),
)
mediated_share = indirect_da.sum("date") / (direct_delta_da + indirect_da.sum("date"))
pc = azp.plot_forest(
xr.Dataset({"mediated_share": mediated_share}),
combined=True,
figure_kwargs={"figsize": (10, 4)},
)
fig = pc.viz["/"]["figure"].values.item()
value_ax = max(fig.axes, key=lambda a: len(a.lines) + len(a.collections))
value_ax.xaxis.set_major_formatter(mtick.PercentFormatter(xmax=1, decimals=0))
fig.suptitle(
"Share of Each Upper Funnel Channel's Effect That Arrives Through SEM",
fontsize=14,
fontweight="bold",
);
And the same numbers in a table:
mediated_share_summary = summary_table(mediated_share)
mediated_share_summary.round(2)
| mean | hdi_low | hdi_high | |
|---|---|---|---|
| channel | |||
| Newspaper | 0.06 | 0.0 | 0.21 |
| Radio | 0.25 | 0.0 | 0.77 |
| TV | 0.12 | 0.0 | 0.40 |
Before reading those shares, where does the mediator’s level come from? Decomposing it explains why the shares above are small.
posterior = mmm.idata["posterior"]
mediator_level = posterior["funnel_mediator"].mean("date")
baseline_share = float(
(posterior["funnel_baseline"] / mediator_level).mean(("chain", "draw"))
)
upper_share = float(
(
(posterior["funnel_demand"].mean("date") - posterior["funnel_baseline"])
/ mediator_level
).mean(("chain", "draw"))
)
sem_share = float(
(
(posterior["funnel_mediator"] - posterior["funnel_demand"]).mean("date")
/ mediator_level
).mean(("chain", "draw"))
)
print(
f"mediator level: baseline {baseline_share:.0%}, upper funnel media {upper_share:.0%}, SEM term {sem_share:.0%}"
)
mediator level: baseline 41%, upper funnel media 10%, SEM term 49%
The split is \(41\%\) baseline, \(10\%\) upper funnel media, and \(49\%\) the SEM term.
That last \(49\%\) is the lambda times spend term at historical SEM levels. In other words, about half of the pool SEM harvests is demand its own spend created. This is exactly what the DAG at the top draws.
So ambient demand and SEM activity carry most of the pool’s level, while the upper channels move it at the margin. This is the quantitative reason the mediated shares are small with wide intervals. There is real mediation here, but it rides on a pool that is mostly not built by the upper channels.
Now we can read the shares. The point estimates put the mediated share at \(6\%\) for Newspaper, \(12\%\) for TV, and \(25\%\) for Radio, with wide \(94\%\) HDIs whose lower ends reach zero (\(0\) to \(21\%\), \(0\) to \(40\%\), and \(0\) to \(77\%\)).
We should read this as a direction, not a precise measurement. The model sees a positive mediated path, but four years of observational data cannot pin down how much of each upper channel’s effect travels through it. This is exactly the quantity a future upper funnel lift test should target.
For the budget decision, the direction is what matters. Cutting Radio loses its direct sales plus about a quarter of its total effect that arrives through search, with an interval that reaches zero. The same logic applies with smaller shares to TV and Newspaper.
And for the attribution dispute, the mediator split gives us the size of the claim. Of the demand pool search harvests, about \(10\%\) is upper funnel media, \(41\%\) is ambient interest, and \(49\%\) is SEM’s own spend. The brand team’s claim on search revenue is real, but modest.
Return on Ad Spend (ROAS)#
To have a more quantitative metric of efficiency, we turn to ROAS. contribution_over_spend() runs the same counterfactual per channel, with two upgrades over our hand-rolled version. It widens the evaluation window so adstock carryover past the last observation is counted, and it does so for both paths because the funnel effect declared its reach through incrementality_spec.
For the four direct-only channels this is the same quantity the first iteration reported. For TV, Radio, and Newspaper it is the total ROAS: direct plus mediated.
roas = mmm.incrementality.contribution_over_spend(
frequency="all_time", num_samples=500, random_state=np.random.default_rng(seed)
).rename("roas")
pc = azp.plot_forest(
xr.Dataset({"roas": roas}), combined=True, figure_kwargs={"figsize": (10, 6)}
)
fig = pc.viz["/"]["figure"].values.item()
fig.suptitle("Return on Ad Spend (funnel-aware)", fontsize=18, fontweight="bold");
The same estimates as a table:
roas_summary = summary_table(roas)
roas_summary.round(2)
| mean | hdi_low | hdi_high | |
|---|---|---|---|
| channel | |||
| Direct Mail | 4.69 | 2.04 | 7.27 |
| Insert | 57.62 | 0.43 | 141.59 |
| Newspaper | 14.86 | 0.37 | 31.97 |
| Online Display | 27.94 | 21.06 | 35.06 |
| Radio | 35.74 | 1.83 | 85.17 |
| Social Media | 26.79 | 0.08 | 67.36 |
| TV | 45.74 | 3.55 | 92.81 |
Is the built-in ROAS doing what our explicit counterfactual did? We can check the two against each other.
comparison_rows = []
for channel in upper_channels:
hand_rolled = float(
(
(direct_deltas[channel] + indirect_da.sel(channel=channel).sum("date"))
/ X[channel].sum()
).mean(("chain", "draw"))
)
comparison_rows.append(
{
"channel": channel,
"incrementality module": float(roas.sel(channel=channel).mean()),
"hand-rolled counterfactual": hand_rolled,
}
)
roas_check = pd.DataFrame(comparison_rows).set_index("channel")
roas_check["relative gap"] = (
roas_check["incrementality module"] / roas_check["hand-rolled counterfactual"] - 1
)
print(
f"largest relative gap between the two ROAS computations: {roas_check['relative gap'].abs().max():.1%}"
)
roas_check.round(2)
largest relative gap between the two ROAS computations: 2.4%
| incrementality module | hand-rolled counterfactual | relative gap | |
|---|---|---|---|
| channel | |||
| Newspaper | 14.86 | 14.81 | 0.00 |
| Radio | 35.74 | 36.64 | -0.02 |
| TV | 45.74 | 46.47 | -0.02 |
The two computations agree to within \(2.4\%\) (the largest relative gap, printed above). The small differences come from the module’s wider evaluation window, which counts the carryover tail past the last observation, and from posterior subsampling. Good, the built-in ROAS is doing exactly what our explicit counterfactual does.
Now the interesting part. How do the two calibrated channels compare against the experiments we fed the model?
pc = azp.plot_dist(
xr.Dataset({"roas": roas}),
coords={"channel": list(roas_experiments)},
ci_prob=HDI_PROB,
ci_kind="hdi",
figure_kwargs={"figsize": (14, 5), "layout": "constrained"},
)
for channel, experiment in roas_experiments.items():
ax = pc.viz["plot"]["roas"].sel(channel=channel).item()
lower = experiment["roas"] - 1.96 * experiment["sigma"]
upper = experiment["roas"] + 1.96 * experiment["sigma"]
ax.axvspan(lower, upper, color="C1", alpha=0.2, label="experiment 95% range")
ax.axvline(experiment["roas"], color="C1", linestyle="--", label="experiment value")
ax.legend(loc="upper right")
ax.set(title=channel, xlabel="ROAS")
fig = pc.viz["/"]["figure"].values.item()
fig.suptitle(
"Calibrated ROAS Posteriors vs the Experiments", fontsize=16, fontweight="bold"
);
Each panel shows the posterior density of the channel’s ROAS, the bar under it (its \(94\%\) HDI, with the posterior mean labelled), and the shaded band (the experiment’s reported \(95\%\) range, its center dashed).
What should you look for here? A posterior inside the experiment range is the expected outcome when the calibration is doing its job, so that is a good start. The information is in the position within the range. A posterior at the center means the observational data and the experiment agree. A posterior pulled against an edge means the two were competing, and that is what failure would have looked like.
Both posteriors sit at their range centers: Online Display at a mean of \(27.9\) (\(94\%\) HDI \(21.1\) to \(35.1\)) against a midpoint of \(28.0\), and Direct Mail at \(4.7\) (\(2.0\) to \(7.3\)) against \(5.0\). The observational data and the experiments agree.
Now contrast these widths with the uncalibrated channels in the table above, where Insert spans \(0.4\) to \(141.6\) and Radio \(1.8\) to \(85.2\). Calibration turns the two decision-critical ROAS numbers from wide bands into precise estimates.
Business Value: uncertainty reduction is what converts a recommendation into an approved budget line. The first iteration could only say “Online Display looks great, but the band is wide”. This iteration can commit to a number.
Budget Optimization#
In this model iteration we work in a more realistic scenario: we have no test set any more, so let’s do what a marketing team actually does, namely plan the next quarter. Our window is the \(13\) weeks after the last observation, 2018-08-05 to 2018-10-28.
We make three planning choices here, and you should feel free to challenge each of them:
Budget anchor. We anchor the weekly pot and the starting plan on the last \(52\) weeks of spend. The last quarter alone would freeze
Newspaperat a near-zero weekly level, which is not a real decision space. The full four years would resurrect a spend regime that no longer exists.One pot, SEM included. SEM is money, so it competes with the seven channels for the same total. The optimizer handles this through
spend_vars, where the SEM budget joins the budget sum constraint and gets optimized as a single weekly level.The objective is a choice. The default objective prices the direct path only. Our model routes part of the media response through the funnel, so the objective we want is
total_response_original_scale. We run both on the same posterior to show that the choice matters.
Note
Requirements for optimizing a funnel model, from Advanced Funnel-Aware MMM: Geo-Level Mediation Through Lower-Funnel Spend: the effect must read the model’s own channel data (ours does), the objective must be named explicitly (total_response_original_scale), and the optimization window arithmetic must respect the effect’s extra carryover (handled automatically by budget_optimizer() via effective_carryover_lags). We compile with the C linker because the default backend cannot compile the gradient of the chained adstock.
First, the historical weekly spend that anchors the plan.
pot_columns = [*channel_columns, "sem_spend"]
weekly_hist = X[pot_columns].tail(52).mean()
media_hist = weekly_hist[channel_columns]
sem_hist = float(weekly_hist["sem_spend"])
total_weekly_budget = float(weekly_hist.sum())
media_pot = float(media_hist.sum())
print(
f"weekly budget pot: {total_weekly_budget:,.0f} (media {media_pot:,.0f} + SEM {sem_hist:,.0f})"
)
weekly_hist.round(0).to_frame(name="weekly spend, last 52 weeks")
weekly budget pot: 2,086,890 (media 1,364,471 + SEM 722,419)
| weekly spend, last 52 weeks | |
|---|---|
| Direct Mail | 566469.0 |
| Insert | 49379.0 |
| Newspaper | 32624.0 |
| Online Display | 281534.0 |
| Radio | 121969.0 |
| Social Media | 175997.0 |
| TV | 136500.0 |
| sem_spend | 722419.0 |
Next we set the guardrails. As in the first iteration, we allow a \(50\%\) change around the historical weekly spend of each channel, plus a band of the same width on the SEM level.
start_date, end_date = "2018-08-05", "2018-10-28" # 13 weeks after the last observation
n_plan_weeks = 13
percentage_change = 0.5
def media_box_bounds(pct: float) -> xr.DataArray:
"""Per-channel weekly bounds: a plus or minus pct box around the historical plan."""
return xr.DataArray(
np.stack(
[(1 - pct) * media_hist.to_numpy(), (1 + pct) * media_hist.to_numpy()],
axis=1,
),
dims=["channel", "bound"],
coords={"channel": channel_columns, "bound": ["lower", "upper"]},
)
# The box bounds reach the media cells only. The SEM budget needs its own band,
# expressed as two scalar inequality constraints; passing any custom constraint
# replaces the default budget-sum constraint, so it is re-added every time.
def sem_band_constraints(pct: float) -> list:
"""Build the budget-sum constraint plus a plus or minus pct band on the SEM weekly level."""
sem_lo = (1 - pct) * sem_hist
sem_hi = (1 + pct) * sem_hist
return [
build_default_sum_constraint("default"),
Constraint(
key="sem_lower",
constraint_type="ineq",
constraint_fun=lambda budgets, total, opt, lo=sem_lo: (
opt.optimization_variables.variable_slice("sem_spend").sum() - lo
),
),
Constraint(
key="sem_upper",
constraint_type="ineq",
constraint_fun=lambda budgets, total, opt, hi=sem_hi: (
hi - opt.optimization_variables.variable_slice("sem_spend").sum()
),
),
]
budget_bounds = media_box_bounds(percentage_change)
Now we are ready to pass all these business requirements to the optimizer.
CVM = {"mode": Mode(linker="cvm")}
optimizer = mmm.budget_optimizer(
start_date=start_date,
end_date=end_date,
response_variable="total_response_original_scale",
spend_vars=["sem_spend"],
constraints=sem_band_constraints(percentage_change),
compile_kwargs=CVM,
)
# The decision vector must reach the mediated path, or the funnel objective
# would silently price the direct path alone.
mediated_response = optimizer.extract_response_distribution(
"funnel_effect_contribution"
)
if optimizer.optimization_variables.flat not in set(ancestors([mediated_response])):
raise RuntimeError("the decision vector does not reach the mediated path")
x0 = {
"channel_data": xr.DataArray(
media_hist.to_numpy(), dims=["channel"], coords={"channel": channel_columns}
),
"sem_spend": xr.DataArray(sem_hist),
}
objective_at_x0 = optimizer.evaluate_plan(x0).objective
# SLSQP's ftol is an absolute precision target. Anchor it to the objective scale,
# and anchor it tightly: the response total is dominated by the baseline, so even
# large reallocation gains are a small relative change of the objective.
ftol = 1e-9 * abs(objective_at_x0)
minimize_kwargs = {"method": "SLSQP", "options": {"ftol": ftol, "maxiter": 2_000}}
result = optimizer.allocate_budget(
total_budget=total_weekly_budget,
budget_bounds=budget_bounds,
x0=x0,
minimize_kwargs=minimize_kwargs,
)
print(
f"solver: {result.scipy_result.message} after {result.scipy_result.nit} iterations"
)
sem_allocation = float(result.spend_var_allocations["sem_spend"].sum())
np.testing.assert_allclose(
float(result.budgets.sum()) + sem_allocation, total_weekly_budget, rtol=1e-6
)
print(
f"funnel-aware optimum: media {float(result.budgets.sum()):,.0f} + SEM {sem_allocation:,.0f}"
)
solver: Optimization terminated successfully after 29 iterations
funnel-aware optimum: media 1,480,238 + SEM 606,652
Direct-Only Comparison#
What would we have planned without the funnel? Let’s run the same posterior with the default objective.
This optimizer cannot see the mediated path, so we must not give it the SEM lever (the library refuses a spend variable the objective cannot see). It re-allocates the media pot while SEM stays at its historical level, so both plans spend the same total.
direct_optimizer = mmm.budget_optimizer(
start_date=start_date,
end_date=end_date,
response_variable="total_media_contribution_original_scale",
compile_kwargs=CVM,
)
x0_media = x0["channel_data"]
direct_objective_at_x0 = direct_optimizer.evaluate_plan(x0_media).objective
result_direct = direct_optimizer.allocate_budget(
total_budget=media_pot,
budget_bounds=budget_bounds,
x0=x0_media,
minimize_kwargs={
"method": "SLSQP",
"options": {"ftol": 1e-9 * abs(direct_objective_at_x0), "maxiter": 2_000},
},
)
print(
f"solver: {result_direct.scipy_result.message} after {result_direct.scipy_result.nit} iterations"
)
np.testing.assert_allclose(float(result_direct.budgets.sum()), media_pot, rtol=1e-6)
solver: Optimization terminated successfully after 31 iterations
Objective vs SEM Lever#
There is a subtlety here. The two runs above differ in two ways at once: the objective, and whether the SEM budget may move at all.
So we add a third run to separate the two effects. We use the funnel objective with SEM pinned at its historical level, which forces exactly the direct run’s media pot. Comparing the pinned run with the direct run isolates the objective, and comparing the full funnel plan with the pinned run isolates the SEM lever.
Tip
Swapping the constraint set with set_constraints reuses the already compiled optimizer, so we do not pay the compilation cost again.
optimizer.set_constraints(
[
build_default_sum_constraint("default"),
Constraint(
key="sem_pinned",
constraint_type="eq",
constraint_fun=lambda budgets, total, opt: (
opt.optimization_variables.variable_slice("sem_spend").sum() - sem_hist
),
),
]
)
result_pinned = optimizer.allocate_budget(
total_budget=total_weekly_budget,
budget_bounds=budget_bounds,
x0=x0,
minimize_kwargs=minimize_kwargs,
)
np.testing.assert_allclose(
float(result_pinned.budgets.sum()) + sem_hist, total_weekly_budget, rtol=1e-6
)
np.testing.assert_allclose(float(result_pinned.budgets.sum()), media_pot, rtol=1e-6)
print(
f"pinned-SEM optimum: media {float(result_pinned.budgets.sum()):,.0f}, the direct run's pot"
)
pinned-SEM optimum: media 1,364,471, the direct run's pot
Before comparing the plans, one structural check. At the historical plan, the difference in marginal returns between the two objectives is exactly the marginal mediated response: positive for the three upper funnel channels, zero for the four direct-only channels.
marginal_total = optimizer.evaluate_plan(x0).utility_gradient["channel_data"]
marginal_direct = direct_optimizer.evaluate_plan(x0_media).utility_gradient[
"channel_data"
]
mediated_marginal = marginal_total - marginal_direct
# Both gradients are scattered over the full channel axis, and a cell outside
# an optimizer's mask comes back as zero. Subtracting across two optimizers is
# the mediated marginal only where both optimize the same cells.
xr.testing.assert_equal(
optimizer.budgets_to_optimize, direct_optimizer.budgets_to_optimize
)
direct_only_channels = [c for c in channel_columns if c not in upper_channels]
if not (mediated_marginal.sel(channel=upper_channels) > 0).all():
raise RuntimeError(
"expected a positive mediated marginal on every upper funnel cell"
)
np.testing.assert_allclose(
mediated_marginal.sel(channel=direct_only_channels).to_numpy(),
0.0,
atol=float(abs(marginal_total).max()) * 1e-9,
)
print(
"gradient decomposition holds: the mediated marginal sits on the upper funnel channels only"
)
gradient decomposition holds: the mediated marginal sits on the upper funnel channels only
That same gradient gives every spending line its marginal return at the historical plan, which is the number the reallocation logic actually runs on. It is also the only place SEM gets a return figure of its own, since the ROAS table earlier covers the seven direct channels while SEM’s value runs through the mediator.
Note
A note on units. The decision variable is a weekly level held for the \(13\) planning weeks, so the raw gradient is response over the planning window (including its carry-over tail) per sustained weekly dollar. Divided by \(13\) it becomes an average per-dollar marginal, directly comparable with the ROAS table above. Every channel lands within about \(15\%\) of its average ROAS near the operating point (the ratio column of the table), with Online Display the widest gap.
marginal = optimizer.evaluate_plan(x0).utility_gradient
per_sustained_weekly_dollar = pd.concat(
[
marginal["channel_data"].to_series(),
pd.Series({"sem_spend": float(marginal["sem_spend"])}),
]
).loc[pot_columns]
marginal_response = pd.DataFrame(
{
"per sustained weekly dollar": per_sustained_weekly_dollar,
"per dollar": per_sustained_weekly_dollar / n_plan_weeks,
}
)
marginal_response["per dollar over average ROAS"] = marginal_response[
"per dollar"
] / roas_summary["mean"].reindex(marginal_response.index)
gap_to_roas = (marginal_response["per dollar over average ROAS"] - 1).abs()
print(
f"largest gap between the per-dollar marginal and the average ROAS: "
f"{gap_to_roas.max():.1%} ({gap_to_roas.idxmax()})"
)
marginal_response.round(2)
largest gap between the per-dollar marginal and the average ROAS: 14.9% (Online Display)
| per sustained weekly dollar | per dollar | per dollar over average ROAS | |
|---|---|---|---|
| Direct Mail | 60.70 | 4.67 | 1.00 |
| Insert | 765.20 | 58.86 | 1.02 |
| Newspaper | 207.11 | 15.93 | 1.07 |
| Online Display | 309.18 | 23.78 | 0.85 |
| Radio | 440.66 | 33.90 | 0.95 |
| Social Media | 307.70 | 23.67 | 0.88 |
| TV | 627.29 | 48.25 | 1.05 |
| sem_spend | 158.54 | 12.20 | NaN |
In per-dollar terms, a marginal SEM dollar returns about \(12\). That is above Direct Mail’s \(4.7\) and below every other media line (Online Display \(24\), TV \(48\), Insert \(59\)).
Keep this ranking in mind. The optimizer moves money from low-marginal lines toward high-marginal ones until a guardrail binds, so the SEM trim we are about to see is exactly this comparison at work. SEM outearns Direct Mail at the margin, but once Direct Mail is pinned at its floor, SEM becomes the weakest line still free to give.
We can now visualize the three plans side by side.
allocation_df = pd.DataFrame(
{
"historical (52w avg)": weekly_hist,
"direct-only objective": pd.concat(
[result_direct.budgets.to_series(), pd.Series({"sem_spend": sem_hist})]
),
"funnel objective, SEM pinned": pd.concat(
[result_pinned.budgets.to_series(), pd.Series({"sem_spend": sem_hist})]
),
"funnel-aware objective": pd.concat(
[result.budgets.to_series(), pd.Series({"sem_spend": sem_allocation})]
),
}
).loc[[*channel_columns, "sem_spend"]]
fig, ax = plt.subplots(figsize=(12, 8))
allocation_df.plot.barh(color=["gray", "C0", "C2", "C1"], ax=ax)
ax.xaxis.set_major_formatter(mtick.FuncFormatter(lambda v, _: f"{v / 1e3:,.0f}K"))
ax.set(title="Weekly Budget Allocation for the Next Quarter", xlabel="weekly spend");
And the same allocations as a table, with the change against the historical plan:
pinned_vs_direct = (
allocation_df["funnel objective, SEM pinned"]
- allocation_df["direct-only objective"]
).abs()
sem_trim = sem_allocation / sem_hist - 1
sem_freed = sem_hist - sem_allocation
print(
f"largest cell difference, pinned vs direct-only: {pinned_vs_direct.max():,.2f} "
f"({pinned_vs_direct.idxmax()})"
)
print(
f"SEM trim: {sem_trim:.1%} ({sem_freed:,.0f} per week freed for the media ceilings)"
)
allocation_df.assign(
**{
"funnel-aware vs historical": allocation_df["funnel-aware objective"]
/ allocation_df["historical (52w avg)"]
- 1
}
).round(2)
largest cell difference, pinned vs direct-only: 0.05 (Online Display)
SEM trim: -16.0% (115,767 per week freed for the media ceilings)
| historical (52w avg) | direct-only objective | funnel objective, SEM pinned | funnel-aware objective | funnel-aware vs historical | |
|---|---|---|---|---|---|
| Direct Mail | 566468.87 | 283234.44 | 283234.44 | 283234.45 | -0.50 |
| Insert | 49378.59 | 74067.88 | 74067.88 | 74067.88 | 0.50 |
| Newspaper | 32623.57 | 16311.79 | 16311.79 | 48935.36 | 0.50 |
| Online Display | 281533.89 | 371293.24 | 371293.20 | 422300.83 | 0.50 |
| Radio | 121968.93 | 182953.39 | 182953.39 | 182953.39 | 0.50 |
| Social Media | 175997.19 | 231860.32 | 231860.36 | 263995.78 | 0.50 |
| TV | 136499.97 | 204749.96 | 204749.96 | 204749.95 | 0.50 |
| sem_spend | 722418.76 | 722418.76 | 722418.76 | 606652.13 | -0.16 |
It also helps to see where each cell of the plan sits inside its own guardrail band.
plan = allocation_df["funnel-aware objective"]
band = np.linspace(1 - percentage_change, 1 + percentage_change, 101)
def spend_dataset(values: np.ndarray) -> xr.Dataset:
"""Wrap weekly spend levels of shape (sample, channel) for plot_forest."""
return xr.Dataset(
{
"weekly spend": xr.DataArray(
values,
dims=["sample", "channel"],
coords={"channel": [*channel_columns, "SEM"]},
)
}
)
# A uniform grid over each guardrail band draws the band as the full-range
# interval; because the band is centered on the historical plan, the mean of
# that grid marks the historical spend itself.
pc = azp.plot_forest(
{
f"allowed range (+/-{percentage_change:.0%}), centered on historical spend": spend_dataset(
np.outer(band, weekly_hist.to_numpy())
),
"funnel-aware optimum": spend_dataset(np.tile(plan.to_numpy(), (len(band), 1))),
},
sample_dims=["sample"],
ci_kind="eti",
ci_probs=(1.0, 1.0),
point_estimate="mean",
labels=["channel"],
visuals={"point_estimate": {"size": 60}},
figure_kwargs={"figsize": (12, 7), "width_ratios": [1, 5]},
)
pc.add_legend("model", title="", loc="outside lower center", ncols=2)
fig = pc.viz["/"]["figure"].values.item()
ax = max(fig.axes, key=lambda a: len(a.lines) + len(a.collections))
ax.set_xscale("log")
ax.set_xticks([25e3, 50e3, 100e3, 200e3, 400e3, 800e3])
ax.xaxis.set_major_formatter(mtick.FuncFormatter(lambda v, _: f"{v / 1e3:,.0f}K"))
ax.xaxis.set_minor_formatter(mtick.NullFormatter())
ax.set_xlabel("weekly spend (log scale)")
fig.suptitle(
"Where the Funnel-Aware Plan Sits Inside Its Guardrails",
fontsize=14,
fontweight="bold",
);
The range plot makes the shape of the solution “obvious”: every media cell in the funnel-aware plan sits on a guardrail. Its orange marker lands on an end of the band, with Direct Mail at its floor and the other six at their ceilings. The media allocation is a corner of the plus or minus \(50\%\) box.
The only genuinely interior number in the plan is SEM, whose marker sits inside its band. It is the residual left once the ceilings are funded. We measure how much of the plan is really the guardrail in the sensitivity section below.
Here is an interesting point: the pinned column is numerically identical to the direct-only column! The two agree to within five cents in every cell, because the objective is flat between two channels it values identically.
So at the same media pot, the funnel-aware and direct-only objectives choose the same corner. The mediated marginals raise the upper channels’ returns without reordering the channels.
Every visible difference between the direct-only and funnel-aware plans is therefore the SEM lever. Trimming SEM by \(16\%\) frees about \(116K\) per week, which funds the remaining ceilings. Newspaper then flips from its floor to its ceiling simply because it is the weakest line still below its cap when the freed money arrives (its marginal of \(207\) per sustained weekly dollar, about \(16\) per dollar, is the lowest media line after Direct Mail).
In this model, the useful role of the funnel is that it tells us how much to spend on SEM, not that it changes the order of the other media channels.
Expected Response Under the Plan#
The optimizer’s own graph gives the posterior response distribution under any plan, including the carryover into and out of the window, through evaluate_response_distribution(). We use it to answer the planning question directly: how much more response should the funnel-aware plan produce than the direct-only plan, when both are scored by the full (funnel-aware) objective?
def report_distribution(
da: xr.DataArray,
label: str,
fmt: str = ",.0f",
dim: list[str] | None = None,
) -> None:
"""Print the mean, the 94% HDI and P(> 0) of a scalar-per-draw quantity.
`dim` names the sample dimensions. It defaults to ArviZ's `chain` and
`draw`; pass `["sample"]` for a quantity that came out of the optimizer's
own graph, where the two are stacked.
"""
hdi = az.hdi(da, prob=HDI_PROB, **({"dim": dim} if dim else {}))
low, high = float(hdi.sel(ci_bound="lower")), float(hdi.sel(ci_bound="upper"))
print(
f"{label}: {float(da.mean()):{fmt}} "
f"(94% HDI [{low:{fmt}}, {high:{fmt}}]), "
f"P(> 0) = {float((da > 0).mean()):.1%}"
)
plans = {
"funnel-aware": {
"channel_data": result.budgets,
"sem_spend": result.spend_var_allocations["sem_spend"],
},
"direct-only": {
"channel_data": result_direct.budgets,
"sem_spend": xr.DataArray(sem_hist),
},
"SEM pinned": {
"channel_data": result_pinned.budgets,
"sem_spend": xr.DataArray(sem_hist),
},
}
response = {
label: optimizer.evaluate_response_distribution(
plan, "total_response_original_scale"
)
for label, plan in plans.items()
}
response_gain = response["funnel-aware"] - response["direct-only"]
gain_from_objective = response["SEM pinned"] - response["direct-only"]
gain_from_sem_lever = response["funnel-aware"] - response["SEM pinned"]
pc = azp.plot_dist(
xr.Dataset({"x": response_gain}),
sample_dims=["sample"],
figure_kwargs={"figsize": (10, 6)},
)
fig = pc.viz["/"]["figure"].values.item()
ax = fig.axes[0]
ax.axvline(x=0, color="C1", linestyle="--", linewidth=2)
ax.set_title(
"Funnel-Aware Plan Minus Direct-Only Plan: Total Response over the Quarter",
fontsize=14,
)
report_distribution(response_gain, "expected gain", dim=["sample"])
print(
f"gain from the objective alone (pinned vs direct): "
f"{float(gain_from_objective.mean()):,.2f} "
f"({float(gain_from_objective.mean() / response_gain.mean()):.2%} of the total)"
)
print(
f"gain from the SEM lever (funnel-aware vs pinned): "
f"{float(gain_from_sem_lever.mean()):,.2f} "
f"({float(gain_from_sem_lever.mean() / response_gain.mean()):.2%} of the total)"
)
expected gain: 9,997,983 (94% HDI [-19,546,052, 38,518,106]), P(> 0) = 75.0%
gain from the objective alone (pinned vs direct): -0.52 (-0.00% of the total)
gain from the SEM lever (funnel-aware vs pinned): 9,997,983.93 (100.00% of the total)
The expected gain of the funnel-aware plan over the direct-only plan is about \(10.0\) million in total response over the quarter (\(94\%\) HDI minus \(19.5\) to plus \(38.5\) million), with a \(75\%\) probability of being the better plan.
The decomposition is one-sided. The objective alone contributes essentially nothing, since the pinned and direct plans coincide, and the SEM lever contributes the entire \(10.0\) million.
Warning
Two caveats keep this number honest:
The funnel-aware plan is by construction the argmax of the metric scoring it, so a positive expected gain is guaranteed. The win probability only reflects posterior variation around that optimum, and is not an external validation.
The gain comes from money the direct framing could not move at all. It measures the value of having the SEM lever, and not the value of a smarter media split.
Guardrail Sensitivity#
Every media cell of the funnel-aware plan sits on a bound. This means the quotable numbers, including the SEM trim, are functions of the \(50\%\) guardrail, a number we picked for operational comfort rather than derived from the model.
So which conclusions are properties of the response surface, and which are properties of the box? We re-solve the same problem under tighter and looser boxes to find out. We reuse the compiled optimizer and only change the constraint set and the bounds per box.
sensitivity_rows = []
for pct in [0.2, 0.5, 1.0]:
if pct == percentage_change:
box_result = result
else:
optimizer.set_constraints(sem_band_constraints(pct))
box_result = optimizer.allocate_budget(
total_budget=total_weekly_budget,
budget_bounds=media_box_bounds(pct),
x0=x0,
minimize_kwargs=minimize_kwargs,
return_if_fail=True,
)
box = media_box_bounds(pct)
at_bound = (
np.isclose(box_result.budgets, box.sel(bound="lower"), rtol=1e-4)
| np.isclose(box_result.budgets, box.sel(bound="upper"), rtol=1e-4)
).sum()
box_sem = float(box_result.spend_var_allocations["sem_spend"].sum())
sensitivity_rows.append(
{
"guardrail": f"+/-{pct:.0%}",
"Direct Mail": float(box_result.budgets.sel(channel="Direct Mail")),
"SEM": box_sem,
"SEM vs historical": box_sem / sem_hist - 1,
"media cells at a bound": int(at_bound),
"solver success": bool(box_result.scipy_result.success),
}
)
pd.DataFrame(sensitivity_rows).round(2)
| guardrail | Direct Mail | SEM | SEM vs historical | media cells at a bound | solver success | |
|---|---|---|---|---|---|---|
| 0 | +/-20% | 453175.10 | 676112.11 | -0.06 | 7 | True |
| 1 | +/-50% | 283234.45 | 606652.13 | -0.16 | 7 | True |
| 2 | +/-100% | 0.00 | 490885.49 | -0.32 | 7 | True |
Let’s read this table carefully: in every box the optimizer picks the same corner, with Direct Mail at its floor and every other media cell at its ceiling. SEM is then the arithmetic residual, the pot minus the media corner, and the table matches that identity to the cent. So the three rows are really one observation.
What survives is the ranking: SEM prices below every media line still free to give (every line except Direct Mail, which is already at its floor), over the whole range we tested.
The minus \(6\%\), minus \(16\%\), and minus \(32\%\) magnitudes, on the other hand, are properties of the box. On a real engagement you should pick the box from operational constraints and report the plan conditional on it.
There is one probe that does move independently of that arithmetic. The SEM band machinery is inactive at every optimum above, as SEM lands interior in all three boxes. But it binds once the band is tighter than about plus or minus \(16\%\), because the residual at the \(50\%\) media box (\(606,652\)) sits \(16\%\) below the historical level.
So we tighten the band to plus or minus \(10\%\). This forces about \(43.5K\) per week back out of media, and which ceiling unwinds is genuine information: the marginal table is evaluated at the historical plan, not at the optimum, so we cannot read the unwind order off it.
optimizer.set_constraints(sem_band_constraints(0.10))
result_tight_band = optimizer.allocate_budget(
total_budget=total_weekly_budget,
budget_bounds=media_box_bounds(percentage_change),
x0=x0,
minimize_kwargs=minimize_kwargs,
return_if_fail=True,
)
tight_sem = float(result_tight_band.spend_var_allocations["sem_spend"].sum())
print(
f"success={result_tight_band.scipy_result.success}, "
f"SEM={tight_sem:,.0f} (band floor {0.9 * sem_hist:,.0f}), "
f"moved back out of media: {tight_sem - sem_allocation:,.0f} per week"
)
unwound = (result_tight_band.budgets - result.budgets).to_series().round(0)
unwound[unwound != 0].to_frame(name="change vs the 50% box optimum")
success=True, SEM=650,177 (band floor 650,177), moved back out of media: 43,525 per week
| change vs the 50% box optimum | |
|---|---|
| channel | |
| Newspaper | -32624.0 |
| Online Display | -6566.0 |
| Social Media | -4335.0 |
The unwind is ordered. Newspaper gives back its entire ceiling premium (\(32,624\), landing exactly on its floor), while Online Display (\(6,566\)) and Social Media (\(4,335\)) supply the remaining \(10,901\), both settling at interior levels below their ceilings.
The order matches the per-dollar ranking, as Newspaper is the weakest funded media line. But the solve adds what that ranking could not tell us: Newspaper does not just shrink, it exits to its floor before the other two concede anything material, because at these levels its saturation curve offers the least response per dollar reclaimed.
In-Sample Uplift#
The first iteration closed by re-weighting the training spend with the optimized allocation and reading the counterfactual uplift. Let’s repeat that exercise, with one improvement.
The first iteration scaled each channel by an optimal-to-historical multiplier, which entangles the answer with how spend drifted over the years. Here we re-weight in share space instead: every week keeps its observed total media pot, and we redistribute that pot across the eight lines (seven channels plus SEM) according to the optimal shares.
This makes the question clean. If we keep the same money, same weeks, but allocate it differently: what changes?
optimal_allocation = pd.concat(
[result.budgets.to_series(), pd.Series({"sem_spend": sem_allocation})]
).loc[pot_columns]
optimal_shares = optimal_allocation / optimal_allocation.sum()
weekly_pot = X[pot_columns].sum(axis=1)
X_weighted = X.copy()
for column in pot_columns:
X_weighted[column] = weekly_pot * optimal_shares[column]
np.testing.assert_allclose(
X_weighted[pot_columns].sum(axis=1).to_numpy(), weekly_pot.to_numpy()
)
np.testing.assert_allclose(
X_weighted[pot_columns].to_numpy().sum(), X[pot_columns].to_numpy().sum()
)
optimal_shares.round(3).to_frame(name="optimal share of the weekly pot")
| optimal share of the weekly pot | |
|---|---|
| Direct Mail | 0.136 |
| Insert | 0.035 |
| Newspaper | 0.023 |
| Online Display | 0.202 |
| Radio | 0.088 |
| Social Media | 0.127 |
| TV | 0.098 |
| sem_spend | 0.291 |
Now we push the re-weighted history through the posterior and read the uplift in media-driven response.
pp_weighted = sample_plan(X_weighted)
def media_driven_response(pp: xr.Dataset) -> xr.DataArray:
"""Total media-driven response: direct channel contributions plus the mediated path."""
return pp["channel_contribution_original_scale"].sum(("date", "channel")) + pp[
"funnel_effect_contribution_original_scale"
].sum("date")
media_uplift = media_driven_response(pp_weighted) / media_driven_response(pp_base) - 1
pc = azp.plot_dist(xr.Dataset({"x": media_uplift}), figure_kwargs={"figsize": (10, 6)})
fig = pc.viz["/"]["figure"].values.item()
ax = fig.axes[0]
ax.axvline(x=0, color="C1", linestyle="--", linewidth=2)
ax.xaxis.set_major_formatter(mtick.PercentFormatter(xmax=1, decimals=0))
ax.set_title("Uplift in Media-Driven Response Under the Optimal Mix", fontsize=16)
report_distribution(media_uplift, "mean uplift", fmt=".1%")
Sampling: [y]
mean uplift: 23.0% (94% HDI [6.5%, 41.4%]), P(> 0) = 99.7%
The optimal mix would have produced a \(23\%\) higher media-driven response on the same historical spend (\(94\%\) HDI \(6\%\) to \(41\%\)), with a \(99.7\%\) probability of a positive uplift. That is a remarkable figure in view of the typical numbers in marketing spend! (of course, this is a synthetic data set 😅).
Note
Note the numerator here: direct channel contributions plus the mediated contribution. The first iteration summed the channel contributions only, which undercounts in a funnel model, because part of what the upper channels earn arrives through SEM.
Finally, let’s look at the uplift in sales.
sales_uplift = (
pp_weighted["y_original_scale"].sum("date")
/ pp_base["y_original_scale"].sum("date")
- 1
)
pc = azp.plot_dist(xr.Dataset({"x": sales_uplift}), figure_kwargs={"figsize": (10, 6)})
fig = pc.viz["/"]["figure"].values.item()
ax = fig.axes[0]
ax.axvline(x=0, color="C1", linestyle="--", linewidth=2)
ax.xaxis.set_major_formatter(mtick.PercentFormatter(xmax=1, decimals=1))
ax.set_title("Uplift in Total Sales Under the Optimal Mix", fontsize=16)
report_distribution(sales_uplift, "mean sales uplift", fmt=".1%")
mean sales uplift: 9.4% (94% HDI [2.4%, 16.9%]), P(> 0) = 99.6%
At the sales line the expected uplift is \(9.4\%\) (\(94\%\) HDI \(2\%\) to \(17\%\)). It is diluted by everything media does not control, as baseline, seasonality, and observation noise all sit in the denominator. The media-driven number above is the decision metric, while this one shows what would actually be visible in the top line.
How does this compare to the first iteration? It reported \(26\%\) and \(10\%\), but on a different channel set and with a different re-weighting recipe. Treat the agreement in magnitude as reassuring, and not as a like-for-like comparison. We make the proper like-for-like comparison in the next section, on the quantity both models define the same way: the ROAS.
Note
We score both plans with the same predictive-noise draws (the same seed), so observation noise largely cancels in this ratio and the spread mostly reflects posterior parameter uncertainty. Scored with independent noise draws, the distribution would be wider.
Warning
Critical Assessment: this is an in-sample counterfactual and not a forecast. The optimal shares push some channels beyond spend levels the model has seen, where the saturation curves are extrapolations. The SEM impressions series stays at its observed values because it is an outcome of the plan, not an input. And the model does not net out the extra SEM spend the new demand would induce. Pricing that feedback is a candidate for the next iteration.
Comparison with the First Iteration#
The first iteration reported its results as figures and did not save its posterior. So we refit its exact specification here (the same seven channels, priors, adstock, saturation, controls, and seasonality, with no SEM, no funnel, and no calibration) on the same \(209\) weeks, and put the two models side by side.
We compare them on two things: what each model can and cannot say, and which numbers actually changed. The ROAS is the right quantity for this because both models define it the same way, as the incremental contribution of a channel divided by its spend.
Warning
Critical Assessment: this is the first iteration’s specification refit on all \(209\) weeks, and not the first iteration’s posterior, which was fit on \(183\) training weeks. Also, a narrower interval is not automatically a more correct one. The calibrated channels are narrower because two extra likelihood terms constrain them, and that is only as good as the experiments behind them. The first iteration’s results are about \(60\%\) baseline share, \(26\%\) and \(10\%\) in-sample uplift.
mmm_baseline = MMM(
model_config=make_model_config(X),
target_column=target_column,
date_column=date_column,
adstock=GeometricAdstock(l_max=l_max),
saturation=LogisticSaturation(),
channel_columns=channel_columns,
control_columns=control_columns,
yearly_seasonality=5,
)
mmm_baseline.build_model(X=X, y=y)
mmm_baseline.add_original_scale_contribution_variable(
var=["channel_contribution", "intercept_contribution", "y"]
)
mmm_baseline.fit(
X=X,
y=y,
target_accept=0.9,
chains=4,
cores=4,
draws=1_000,
tune=1_500,
nuts_sampler="nutpie",
random_seed=seed,
)
n_divergences_baseline = int(
mmm_baseline.idata["sample_stats"]["diverging"].sum().item()
)
print(f"baseline model divergences: {n_divergences_baseline}")
if n_divergences_baseline != 0:
warnings.warn(
f"the baseline sampler reported {n_divergences_baseline} divergences",
stacklevel=1,
)
NUTS[nutpie]: [y_sigma, gamma_fourier, gamma_control, adstock_alpha, saturation_lam, saturation_beta, intercept_contribution]
baseline model divergences: 0
With both models fitted, we can put their ROAS estimates on the same plot.
roas_baseline = mmm_baseline.incrementality.contribution_over_spend(
frequency="all_time", num_samples=500, random_state=np.random.default_rng(seed)
).rename("roas")
model_labels = ["Part I specification", "Part II funnel-aware, calibrated"]
roas_by_model = xr.concat(
[roas_baseline, roas], dim=pd.Index(model_labels, name="model")
)
pc = azp.plot_forest(
xr.Dataset({"roas": roas_by_model}),
combined=True,
aes={"color": ["model"]},
figure_kwargs={"figsize": (11, 8)},
)
pc.add_legend("model")
fig = pc.viz["/"]["figure"].values.item()
fig.suptitle(
"ROAS: Part I Specification vs Part II",
fontsize=16,
fontweight="bold",
);
The table makes the shifts and the width changes explicit.
roas_part1 = summary_table(roas_baseline)
roas_comparison = pd.concat(
{"Part I specification": roas_part1, "Part II": roas_summary}, axis=1
)
roas_comparison[("comparison", "mean shift (II minus I)")] = (
roas_summary["mean"] - roas_part1["mean"]
)
roas_comparison[("comparison", "HDI width ratio (II over I)")] = (
roas_summary["hdi_high"] - roas_summary["hdi_low"]
) / (roas_part1["hdi_high"] - roas_part1["hdi_low"])
roas_comparison.round(2)
| Part I specification | Part II | comparison | ||||||
|---|---|---|---|---|---|---|---|---|
| mean | hdi_low | hdi_high | mean | hdi_low | hdi_high | mean shift (II minus I) | HDI width ratio (II over I) | |
| channel | ||||||||
| Direct Mail | 7.26 | 0.00 | 19.50 | 4.69 | 2.04 | 7.27 | -2.57 | 0.27 |
| Insert | 55.89 | 0.13 | 125.93 | 57.62 | 0.43 | 141.59 | 1.73 | 1.12 |
| Newspaper | 14.43 | 0.19 | 30.07 | 14.86 | 0.37 | 31.97 | 0.43 | 1.06 |
| Online Display | 56.41 | 0.75 | 118.76 | 27.94 | 21.06 | 35.06 | -28.47 | 0.12 |
| Radio | 32.07 | 0.04 | 84.66 | 35.74 | 1.83 | 85.17 | 3.68 | 0.98 |
| Social Media | 30.65 | 0.19 | 73.95 | 26.79 | 0.08 | 67.36 | -3.86 | 0.91 |
| TV | 43.58 | 0.69 | 86.02 | 45.74 | 3.55 | 92.81 | 2.16 | 1.05 |
And we can compare how the two models split total sales.
baseline_components = {
"baseline (intercept)": "intercept_contribution",
"direct channels": "channel_contribution",
}
share_comparison = pd.DataFrame(
{
"Part I specification": contribution_shares(mmm_baseline, baseline_components),
"Part II": funnel_shares.loc[
["baseline (intercept)", "direct channels", "mediated SEM (funnel)"]
],
}
)
share_comparison.loc["media-driven (direct plus mediated)"] = share_comparison.drop(
"baseline (intercept)"
).sum()
share_comparison.round(3)
| Part I specification | Part II | |
|---|---|---|
| baseline (intercept) | 0.625 | 0.561 |
| direct channels | 0.369 | 0.289 |
| mediated SEM (funnel) | NaN | 0.147 |
| media-driven (direct plus mediated) | 0.369 | 0.437 |
Before reading the numbers, it is worth asking what each model is even able to say.
The first iteration’s model answers one question per channel: what did a dollar of spend return through its own path? It has no place for SEM, no experiment, and one validation split.
The second iteration’s model answers three questions: the same direct one, how much of an upper funnel channel’s return arrives through search, and what a marginal SEM dollar is worth. It does so with two experiments as extra data and four validation folds. The forest and the two tables above measure what that buys us.
Precision: the Calibrated Channels#
The two calibrated lines are where the gain is largest, and where it matters most for the budget.
Online Display moves from a mean of \(56.4\) with a \(94\%\) HDI of \(0.8\) to \(118.8\) under the first specification, to \(27.9\) with an HDI of \(21.1\) to \(35.1\). The interval is \(12\%\) of its former width. Direct Mail moves from \(7.3\) (\(0.0\) to \(19.5\)) to \(4.7\) (\(2.0\) to \(7.3\)), \(27\%\) of the former width.
Both first-iteration intervals reach zero, so on that model alone you could not defend either channel in a budget meeting. The calibrated intervals exclude zero by a wide margin.
The narrowing is by construction, since we added two extra likelihood terms, and the calibrated means moved toward the experiment centers. That is exactly the point: an experiment turns a band into a number a finance team can approve.
Note also the direction of the correction for Online Display. The first specification put its return at \(56\), about twice the experiment’s center. So the wide band was not just imprecise, its center was misleading.
Interpretability: the Upper Funnel Channels#
TV, Radio, and Newspaper have similar widths in both models (ratios \(1.05\), \(0.98\), and \(1.06\)), because neither model has an experiment on them. What changes is the meaning of the number.
The first specification’s figure is a direct-only return that cannot say where the money came back. The second is direct plus mediated, and the mediated share section above tells us how the two parts split.
The mean shifts follow those mediated shares. Radio, with the largest mediated share (\(25\%\)), gains the most, moving from \(32.1\) to \(35.7\), a shift of \(3.7\) ROAS points. TV (\(12\%\) mediated) and Newspaper (\(6\%\)) shift by only \(2.2\) and \(0.4\) points.
The shifts are smaller than the shares alone would predict, because the funnel model also reallocates part of the pool to ambient demand and to SEM’s own spend (the \(41\%\) and \(49\%\) in the mediator split). This is the honest reading: the brand team’s claim on search revenue is real, quantified, and modest.
The Control Group#
Insert and Social Media are direct-only in both models and calibrated in neither, so they should not move. And they do not! Their means shift by \(+1.7\) (Insert) and \(-3.9\) (Social Media) ROAS points, which is small against posterior widths of about \(140\) and \(67\) points, and their HDI width ratios are \(1.12\) and \(0.91\).
This is our check that the funnel and the calibration changed what they were meant to change, and nothing else.
What Only the Second Iteration Has#
SEM does not appear in the first specification at all. That means \(26.9\%\) of the media budget had no return figure and no seat in the budget optimization.
The second iteration prices a marginal SEM dollar at about \(12\), puts SEM in the pot, and finds that the funnel’s planning value in this dataset is exactly that pricing.
The decomposition table tells the same story from the sales side. The first specification credits the direct channels with \(36.9\%\) of sales and the baseline with \(62.5\%\). The funnel model credits the channels with \(28.9\%\) directly plus \(14.7\%\) through search, a media-driven total of \(43.7\%\), and the baseline with \(56.1\%\). So the flat model both undercounts media and misattributes the part that travels through search.
Why the Funnel Results Are More Actionable#
A recommendation is actionable when it names a number, states the condition it depends on, and survives a validation the stakeholder accepts.
The second iteration’s plan (cut Direct Mail to its guardrail floor, fund the other ceilings, trim SEM by \(16\%\)) rests on two ROAS values that are commitments rather than bands, on a funnel checked against the search team’s own impressions series (\(204\) of \(209\) weeks inside the \(94\%\) HDI), and on a specification whose parameters held steady across four quarterly folds.
The first iteration’s recommendation (increase Online Display and Insert, reduce Direct Mail) pointed the same way. But it rested on intervals that included zero and on one validation split, and it could not price the second largest line in the budget at all.
The gain in precision is what lets the plan be approved. The gain in interpretability is what lets the brand and search teams agree on it.
Conclusion#
This second iteration answered the two questions the stakeholders brought back. SEM entered the model as a demand harvester rather than a flat channel, two experiments calibrated the ROAS of Online Display and Direct Mail, four quarterly cross-validation folds replaced the single split, the final model used all \(209\) weeks, and we planned the next quarter’s budget with a funnel-aware objective that co-optimizes the SEM budget.
Key Findings#
The funnel is supported, but read its evidence precisely.
funnel_lambda(the spend-to-impressions link) is tight at \(0.67\) (\(94\%\) HDI \(0.59\) to \(0.75\)) and stable across all four validation folds. The conversion parametersat_lf_betais positive-leaning but wide. The mediator’s level is mostly ambient demand and SEM activity (\(41\%\) baseline, \(10\%\) upper funnel media, \(49\%\) SEM term), so the mediated shares come out positive but imprecise (\(6\%\) to \(25\%\), with \(94\%\) HDIs reaching zero).The measurement model holds up. It covers \(204\) of \(209\) observed impression weeks inside its \(94\%\) HDI, with a slow level drift as its documented gap.
Calibration works. It pinned the two decision-critical ROAS values at the centers of their experiment ranges,
Online Displayat \(27.9\) (\(21.1\) to \(35.1\)) andDirect Mailat \(4.7\) (\(2.0\) to \(7.3\)), while uncalibrated channels still span order-of-magnitude ranges. Against the first iteration’s specification refit on the same data, those two intervals shrank to \(12\%\) and \(27\%\) of their former widths, the upper funnel means moved up in step with their mediated shares (Radioby \(3.7\) ROAS points), and the direct-only channels did not move.The quarter plan is a corner solution. Cut
Direct Mailto its guardrail floor, fund every other channel’s ceiling, and trim SEM by \(16\%\). The whole expected gain over a direct-only plan (about \(10\) million, \(75\%\) win probability) comes from the SEM lever, not from the objective reordering media. The SEM trim follows from SEM pricing below every free media line at every guardrail we tested, with its size set by the guardrail.The optimal mix is worth real money. Re-weighting the whole history to it would have lifted the media-driven response by \(23\%\) (\(99.7\%\) probability positive) and sales by \(9.4\%\).
Model Limitations#
One shared demand pool serves all three upper channels.
A single measurement series (SEM impressions) anchors the mediator.
The calibration constrains direct paths only.
The in-sample uplift extrapolates saturation curves beyond observed spend.
The uplift figures are not like-for-like with the first iteration because the decomposition changed. The ROAS comparison is like-for-like, because both models define ROAS the same way.
Recommendations#
Fund the funnel-aware plan for one quarter with the \(50\%\) guardrails in place.
Log realized SEM impressions against the model’s \(94\%\) HDI predictive band, as a live check the search team can own. Today \(204\) of \(209\) weeks sit inside it.
Hold out a region for the
Direct Mailreduction if operations allow it.
Next Steps#
Where would a third iteration go?
Price the induced SEM spend. Demand created by upper funnel media pulls SEM spend along with it. We should report ROAS per incremental dollar including that feedback, as Advanced Funnel-Aware MMM: Geo-Level Mediation Through Lower-Funnel Spend does.
Use real experiments. Replace the imagined ranges with the real studies they stand in for, and extend the calibration to an upper funnel channel. There the experiment measures the total effect (direct plus mediated), so the calibration target must span both paths. That deserves a notebook of its own, and Lift Test Calibration covers the lift-test mechanics it will build on.
Score the plan. After the quarter, score the realized \(13\) weeks against the plan’s response distribution (CRPS) and roll the cross-validation forward one fold.
Structural extensions. A time-varying demand baseline, separate branded and non-branded demand pools, and geo-level disaggregation if regional data become available.
Operationalize. Track iterations with MLflow. The appendix below shows the YAML registry version of this model, following Building funnel-aware MMMs from YAML, which is the first step of that handover.
Appendix: Building This Model from a YAML Registry#
A model that lives only as notebook code is hard to hand over. build_mmm_from_yaml() builds a model from a configuration file, and Building funnel-aware MMMs from YAML explains the file format.
Here we do it for this exact model and check that the result is the same graph: the same named variables, and the same log-probability at the same point. That check matters, because a YAML file with a subtly different prior or a missing effect would build and fit without complaint.
Two conventions of the builder decide the shape of the file:
Every object is a
classname plus itskwargs, nested as deep as the constructor arguments go. A class is found either by a short name registered on the builder’sREGISTRYor by its dotted import path.Every prior is written in the dictionary form
Prior.to_dict()produces (dist,kwargs,dims), so we can serialize the priors the notebook already holds instead of retyping them.
We register the two classes defined in this notebook (FunnelEffect and FunnelMMM) by short name. We compute the data-derived parts (the spend-share prior widths and the two scaling constants) in Python: the scales go into the file as numbers, and model_config is passed through model_kwargs, which override the file. Keeping data-derived numbers out of the file is a choice here, not a requirement.
Important
Registering notebook-local classes on REGISTRY is a demonstration pattern. A production registry puts FunnelEffect and FunnelMMM in an importable module and references them by dotted path, so that the YAML file, the saved model, and the prediction service all resolve the same class. The prediction-time override in FunnelMMM is still needed until #2945 lands.
REGISTRY["FunnelEffect"] = FunnelEffect
REGISTRY["FunnelMMM"] = FunnelMMM
def class_path(obj) -> str:
"""Dotted import path of an object's class."""
return f"{type(obj).__module__}.{type(obj).__name__}"
def priors_spec(component) -> dict:
"""Serialize a component's priors in the dictionary form the builder reads."""
return {name: prior.to_dict() for name, prior in component.function_priors.items()}
def transformation_spec(transformation: MediaTransformation) -> dict:
"""Build spec for a MediaTransformation, nested class by class."""
adstock, saturation = transformation.adstock, transformation.saturation
return {
"class": class_path(transformation),
"kwargs": {
"adstock": {
"class": class_path(adstock),
"kwargs": {
"l_max": adstock.l_max,
"prefix": adstock.prefix,
"priors": priors_spec(adstock),
},
},
"saturation": {
"class": class_path(saturation),
"kwargs": {
"prefix": saturation.prefix,
"priors": priors_spec(saturation),
},
},
"adstock_first": transformation.adstock_first,
"dims": list(transformation.dims),
},
}
def to_plain(value):
"""Turn tuples into lists recursively so the file is plain YAML."""
if isinstance(value, dict):
return {key: to_plain(item) for key, item in value.items()}
if isinstance(value, list | tuple):
return [to_plain(item) for item in value]
return value
registry_config = to_plain(
{
"model": {
"class": "FunnelMMM",
"kwargs": {
"date_column": date_column,
"target_column": target_column,
"channel_columns": channel_columns,
"control_columns": control_columns,
"yearly_seasonality": 5,
"adstock": {
"class": class_path(mmm.adstock),
"kwargs": {"l_max": l_max},
},
"saturation": {"class": class_path(mmm.saturation), "kwargs": {}},
},
},
"extra_vars": funnel_vars,
"effects": [
{
"class": "FunnelEffect",
"kwargs": {
"data_vars": funnel_vars,
"prefix": "funnel",
"upper_channels": upper_channels,
"sem_scale": funnel_effect.sem_scale,
"impressions_scale": funnel_effect.impressions_scale,
"upper_transform": transformation_spec(
funnel_effect.upper_transform
),
"demand_transform": transformation_spec(
funnel_effect.demand_transform
),
},
}
],
"original_scale_vars": [
"channel_contribution",
"control_contribution",
"intercept_contribution",
"yearly_seasonality_contribution",
"funnel_effect_contribution",
"y",
],
}
)
config_path = Path(tempfile.mkdtemp()) / "mmm_case_study2.yml"
config_path.write_text(yaml.safe_dump(registry_config, sort_keys=False))
print(config_path.read_text())
model:
class: FunnelMMM
kwargs:
date_column: date
target_column: sales
channel_columns:
- Direct Mail
- Insert
- Newspaper
- Online Display
- Radio
- Social Media
- TV
control_columns:
- hldy_Black Friday
- hldy_Christmas Day
- hldy_Christmas Eve
- hldy_Columbus Day
- hldy_Cyber Monday
- hldy_Day after Christmas
- hldy_Easter
- hldy_Father's Day
- hldy_Green Monday
- hldy_July 4th
- hldy_Labor Day
- hldy_MLK
- hldy_Memorial Day
- hldy_Mother's Day
- hldy_NYE
- hldy_New Year's Day
- hldy_Pre Thanksgiving
- hldy_Presidents Day
- hldy_Prime Day
- hldy_Thanksgiving
- hldy_Valentine's Day
- hldy_Veterans Day
yearly_seasonality: 5
adstock:
class: pymc_marketing.mmm.components.adstock.GeometricAdstock
kwargs:
l_max: 6
saturation:
class: pymc_marketing.mmm.components.saturation.LogisticSaturation
kwargs: {}
extra_vars:
- sem_spend
- sem_impressions
effects:
- class: FunnelEffect
kwargs:
data_vars:
- sem_spend
- sem_impressions
prefix: funnel
upper_channels:
- Newspaper
- Radio
- TV
sem_scale: 3134565.13
impressions_scale: 790102.0
upper_transform:
class: pymc_marketing.mmm.media_transformation.MediaTransformation
kwargs:
adstock:
class: pymc_marketing.mmm.components.adstock.GeometricAdstock
kwargs:
l_max: 6
prefix: adstock_uf
priors:
alpha:
dist: Beta
kwargs:
alpha: 2
beta: 3
dims:
- channel_upper
saturation:
class: pymc_marketing.mmm.components.saturation.LogisticSaturation
kwargs:
prefix: sat_uf
priors:
lam:
dist: Gamma
kwargs:
mu: 2.0
sigma: 1.0
dims:
- channel_upper
beta:
dist: HalfNormal
kwargs:
sigma: 1.0
dims:
- channel_upper
adstock_first: true
dims:
- channel_upper
demand_transform:
class: pymc_marketing.mmm.media_transformation.MediaTransformation
kwargs:
adstock:
class: pymc_marketing.mmm.components.adstock.GeometricAdstock
kwargs:
l_max: 6
prefix: adstock_lf
priors:
alpha:
dist: Beta
kwargs:
alpha: 2
beta: 3
saturation:
class: pymc_marketing.mmm.components.saturation.LogisticSaturation
kwargs:
prefix: sat_lf
priors:
lam:
dist: Gamma
kwargs:
mu: 1.5
sigma: 0.5
beta:
dist: HalfNormal
kwargs:
sigma: 1.0
adstock_first: true
dims: []
original_scale_vars:
- channel_contribution
- control_contribution
- intercept_contribution
- yearly_seasonality_contribution
- funnel_effect_contribution
- y
Now we build the model from that file and compare it against the one we fitted above.
mmm_yaml = build_mmm_from_yaml(
config_path,
X=make_funnel_dataset(X, y),
model_kwargs={"model_config": make_model_config(X)},
)
mmm_yaml.add_cost_per_target_calibration(
data=X,
calibration_data=roas_calibration_df,
name_prefix="roas_calibration",
target_column="roas",
target_per_cost=True,
)
python_vars = set(mmm.model.named_vars)
yaml_vars = set(mmm_yaml.model.named_vars)
if python_vars != yaml_vars:
raise RuntimeError(f"variable sets differ: {sorted(python_vars ^ yaml_vars)}")
point = mmm.model.initial_point()
logp_python = float(mmm.model.compile_logp()(point))
logp_yaml = float(mmm_yaml.model.compile_logp()(point))
np.testing.assert_allclose(logp_yaml, logp_python)
print(
f"named variables in both models: {len(python_vars)}, differences: {len(python_vars ^ yaml_vars)}"
)
print(
f"log-probability at the initial point: python {logp_python:,.3f}, "
f"yaml {logp_yaml:,.3f}, absolute difference {abs(logp_yaml - logp_python):.2e}"
)
named variables in both models: 43, differences: 0
log-probability at the initial point: python -9,551.650, yaml -9,551.650, absolute difference 0.00e+00
The two builds agree! There are \(43\) named variables in both graphs with no differences, and the log-probability at the same initial point matches to the last digit (absolute difference printed above). The YAML file describes this model, and not a near relative of it.
%load_ext watermark
%watermark -n -u -v -iv -w -p pymc_marketing,pytensor,nutpie
Last updated: Tue, 22 Sep 2026
Python implementation: CPython
Python version : 3.14.2
IPython version : 9.15.0
pymc_marketing: 1.1.0
pytensor : 3.3.2
nutpie : 0.16.11
arviz : 1.2.0
arviz_plots : 1.2.0
graphviz : 0.21
matplotlib : 3.11.1
numpy : 2.4.6
pandas : 3.0.5
pydantic : 2.13.4
pymc : 6.3.2
pymc_extras : 0.15.1
pymc_marketing: 1.1.0
pytensor : 3.3.2
seaborn : 0.13.2
xarray : 2026.4.0
yaml : 6.0.3
Watermark: 2.6.0