Advanced Funnel-Aware MMM: Geo-Level Mediation Through Lower-Funnel Spend#
Upper-funnel media create demand, and part of this demand converts later through lower-funnel spend such as paid search. A media mix model can treat this lower-funnel spend as one more channel, leave it out, or control for a symptom of demand instead. Each choice changes the ROAS of the upper-funnel channels, and with it the budget plan.
The Funnel-Aware MMM: Upper to Lower-Funnel Mediation via a custom MuEffect notebook introduced the core idea: a custom MuEffect with its own likelihood encodes the mediation in the model. That example kept the funnel minimal, with one upper-funnel channel, one mediator and one national time series. In this notebook, we scale the same idea up to a setup that is closer to production work.
Here we are interested in the following question: how do the ROAS estimates and the budget plan change when a model mishandles the funnel? We answer it on a synthetic weekly geo panel with a known ground truth, so we can score every model against the truth.
Five features make this example advanced:
Several upper-funnel channels:
tv_spendandsocial_spendfeed the same latent demand pool, each with its own carryover and saturation.A geo panel: we fit the model across three regions at once, with some parameters pooled across geos and others estimated per geo.
A mediator that is paid media: lower-funnel spend (
lower_spend) is paid search, and it converts. Part of it responds to the demand that upper-funnel activity creates, and an exogenous budget of its own drives the rest.A demand indicator: branded search volume (
search_volume) measures the same latent demand but causes nothing. Most of the modeling problem is to keep a cause and an indicator apart.Confounded media plans: category demand and the media plans share an annual cycle, a growth trend and a small amount of correlated noise. Every model absorbs the cycle, only some models adjust for the trend, and no model can observe the noise.
We compare six models on the same data: five naive specifications, which mishandle the funnel in different ways, and a funnel model that encodes the mediation explicitly. Then we put the fitted funnel model through the budget optimizer with two objectives, and we check both plans against the truth.
The notebook has the following sections:
Funnel Structure: the causal graph, its equations and the identification argument.
Candidate Models: the six model structures that we compare.
Prepare Notebook: libraries and configuration.
Model Specification: the custom
BudgetDualIndicatorEffect.Data Generating Process: synthetic data with a known causal decomposition.
Exploratory Data Analysis: the observed series and the latent layer.
D-Separation Analysis: which effect each model identifies, computed on the graph.
Model Fitting: the six fits.
Model Diagnostics: divergences and \(\hat{R}\).
Posterior Predictive Checks: the three likelihoods of the funnel model.
Parameter Recovery: the funnel model against the true parameters.
Media Deep Dive: contributions, ROAS, bias and interval coverage of every model against the truth.
Predictive Fit: why in-sample \(R^2\) cannot rank the models.
Built-In Decomposition: what the library decomposition shows for a custom effect.
Funnel-Aware Budget Optimization: the same fitted model under two objectives, and a joint optimization of the lower funnel.
Conclusion: key findings, model limitations and recommendations.
Next Steps: extensions of this example.
References: the case-study series and the causal inference texts.
Note
This notebook builds on other notebooks. See Funnel-Aware MMM: Upper to Lower-Funnel Mediation via a custom MuEffect for the basic funnel effect, Generating Synthetic MMM Data for the do-operator workflow that generates the data, and Introduction to Understanding Causal Relationships in Media Mix Modeling for adjustment sets in PyMC-Marketing. The References section lists the Nürnberger Versicherung case-study series, which works through the same funnel problem in production.
Funnel Structure#
We model a weekly geo panel. Each model in this notebook is a claim about which causal arrows connect media, demand and sales, so we start with the causal graph (DAG) that generates the data. We refer back to it throughout the notebook.
Tip
How to read the figure: filled nodes are observed, so they appear in a dataset. A dashed outline marks the three quantities that no dataset contains (\(D\), \(M^{*}\) and \(N\)).
Dotted arrows are deterministic calendar components. Dashed arrows are measurement or indicator links: a latent quantity causes its noisy record or its symptom, and d-separation treats them as ordinary edges. Solid arrows are the other structural edges.
Within every geo \(g\) and week \(t\), the graph has the following nodes, from media to target:
tv_spend,social_spend(\(U_{t,g,c}\)): the two upper-funnel channels.category_demand(\(C_{t,g}\)): an exogenous, seasonal index of how much the whole category is in-market, for example a category-level search-trend index. It drives lower-funnel demand but does not act on the target directly. This missing arrow is an exclusion restriction on the mediator equation.lf_budget(\(B_{t,g}\)): an exogenous lower-funnel budget of always-on spend plus promotional pushes. The search team sets it on its own calendar, not in response to demand.\(D_{t,g}\): latent demand, which we never observe.
\(M^{*}_{t,g}\): structural lower-funnel spend, the quantity that buys conversions. The column
lower_spendis this quantity observed with reporting noise.search_volume(\(S_{t,g}\)): branded search volume, an indicator of \(D\).y(\(Y_{t,g}\)): the target.
The arrows in the figure correspond to four equations. First, upper-funnel activity and the category cycle create demand:
Second, the search team’s spend responds to this demand and to its own budget. Both parts are paid media, and both convert:
Third, branded search buys nothing. It is a pure indicator, because people search for the brand when they are in-market:
Finally, the target responds to lower-funnel spend, to the direct effect of each upper-funnel channel and to its own seasonality:
In the code, \(\gamma\) is budget_dual_parent_coef, \(\lambda\) is budget_dual_budget_coef and \(\kappa\) is budget_dual_indicator_coef.
Therefore, every upper-funnel channel reaches the target along two paths: directly, and indirectly along \(U \to D \to M^{*} \to Y\).
Here are two important remarks about the figure.
First, \(M^{*}\) and lower_spend are separate nodes on purpose. A dataset contains only the filled node at the end of the dashed measurement edge. This gap between the structural quantity and its noisy observation is why a model with lower_spend as a channel estimates the direct effect only approximately.
Second, the market driver appears as three nodes, and we observe only two of them:
season: the annual cycle. It is a deterministic function of the calendar, and every model below includes a Fourier basis for it.t: the common growth trend of the media plans and the category index. It is also deterministic and observable, but an annual Fourier basis cannot represent a trend, so the seasonal term does not adjust for it.N: the shared noise of the drivers, the only strand that nobody observes. It is small but not zero, and we measure it after we draw the covariates.
The three strands play different roles in the models below. Every model adjusts for the season, and only a model with a category-demand control or a trend control adjusts for the trend. No model can adjust for \(N\).
If we merged the three strands into one seasonality node, the graph would hide the confounding of a model without the mediator. It would also mark a trend control as a valid adjustment set, although a trend control is valid only if the noise strand is negligible.
Note that a confounder can bias a model even when we observe it, if the model leaves it out.
Note
Why \(N\) is on the graph: the covariate simulation below gives the two media plans and the category index a shared error term (pm.LKJCholeskyCov). For d-separation, a shared error term is equivalent to an unobserved common cause of all three drivers.
A single latent node is a rank-one stand-in for a general \(3 \times 3\) correlation, and splitting \(N\) into three pairwise causes changes none of the verdicts below. The correlation is small, but any nonzero shared noise adds \(N\) to the graph, and \(N\) changes the d-separation verdicts. The size of the correlation sets how much bias the open path can carry.
Mediator and Indicator#
Lower-funnel spend and branded search both sit downstream of demand, but only lower-funnel spend has an arrow into the target.
Lower-funnel spend belongs in the target equation because it is paid media: money goes in and conversions come out. A model that leaves it out misses a driver of sales.
Lower-funnel spend is a mediator and not an ordinary channel, because its level is not chosen freely. Part of it responds to the demand that upper-funnel activity created, and the credit for that part belongs to the upper-funnel channels.
Branded search volume does not belong in the target equation. Nobody buys search volume, because it is a symptom of being in-market. An arrow \(S \to Y\) would invent a causal channel that does not exist.
A model that uses branded search as a channel blocks part of the path that we want to measure, as we show below.
Important
The practical test: ask “if we doubled this variable, would sales move?” Both variables predict sales well, so predictive power cannot separate a spend line from a tracking metric. The doubling question separates them, and model fit cannot answer it.
Note
Latent scale anchor: \(D\) is latent, so the model needs a fixed unit for it. We fix the coefficient on \(D\) in the \(M^{*}\) equation at one, so we measure latent demand in units of the lower-funnel spend that it induces.
Then \(\kappa\) is free, and \(\lambda\) converts the budget into the same units. If that coefficient were free too, we could rescale \(D\) and compensate in \(\kappa\) without changing the likelihood.
Confounding and Identification#
Two more features make this example harder than the single-channel case.
First, the media plans are confounded. The media plans and category demand follow the same market rhythm, so upper-funnel spend and category demand move together through the three strands above (season, trend and noise).
Before we fit anything, we measure how much confounding each strand carries. Adjustment can remove the trend-driven part of category demand. It cannot remove the demand that each channel creates downstream, and the likelihood decides which response curve absorbs that demand.
Second, both channels push on one shared, saturating demand pool, so their indirect effects do not add up. Switching off tv_spend alone and switching off social_spend alone remove less demand in total than switching off both at once.
Each channel is partly shielded by the contribution of the other channel further along the concave curve. We quantify this below.
Identification rests on two exclusion restrictions, that is, two arrows that the graph leaves out on purpose. They do not carry the same weight:
Budget restriction: the lower-funnel budget \(B\) enters the spend equation but not the search equation. This restriction identifies \(\lambda\), because it lets the model separate \(\lambda B\) from \(D\) instead of merging the two into one level. Nothing else in the model can replace it.
Category restriction: category demand \(C\) enters the demand equation but not the target equation. This restriction is a modeling convenience, because \(C\) is observed. We could absorb a direct path \(C \to Y\) with a control column for \(C\) without any loss of identification.
Parameter Pooling#
We replicate the whole graph across geos. Which parameters we share across geos and which we estimate per geo is a modeling choice. With about \(130\) weeks per geo, this choice has a large effect on precision:
Parameter |
Dimensions |
Why |
|---|---|---|
|
|
Carryover and curvature are properties of the medium, not of the region. Pooling them improves precision. |
|
|
The response amplitude differs by region (market size, competition). |
|
|
Regional levels and volatility. How strongly the budget converts into spend is a regional operating choice. |
|
|
The conversion path from lower-funnel spend to the target. How fast paid-search clicks convert and how hard they saturate depends on the regional market. |
|
scalar |
There is only one mediator, so there is no channel axis. With one series per geo, we pool the curvature and let the per-geo amplitude carry the regional differences. |
|
scalar |
The search index is the same instrument everywhere, so we share its loading. |
Here pymc.dims xtensors help. They align operands by dimension name, so a ("channel",) decay multiplies a ("date", "geo", "channel") media array without any manual broadcasting.
Candidate Models#
We compare six model structures, so let’s define them before any code runs. Every model is a geo-level MMM with the same adstock and saturation families, the same priors, the same Fourier seasonal basis and the same sampler settings.
Every model includes the two upper-funnel channels. The models differ only in what else they condition on, and how.
Model |
Beyond the two upper channels |
How it treats the funnel |
Question it can answer |
|---|---|---|---|
|
adds |
conditions on the mediator, which blocks the indirect path |
approximately the direct effect of the upper channels |
|
nothing |
ignores the mediator |
the total effect, but a backdoor through the shared market driver stays open |
|
adds |
ignores the mediator, but blocks that backdoor through the observed demand index |
the total effect, a valid adjustment set on the graph |
|
adds the growth trend |
blocks only the trend strand of the backdoor, so the shared noise of the drivers stays open |
the total effect, valid only if the shared noise is negligible |
|
adds |
conditions on a symptom of demand, which partially blocks the indirect path |
neither effect cleanly |
|
the |
models the mediation explicitly |
direct, indirect and total, separately |
The D-Separation Analysis checks the last column on the graph, and the later sections test it against the known ground truth.
Only the Funnel model goes on to the budget optimization at the end. There we compare what the same fitted model recommends under two objectives: one that ignores the mediated path and one that includes it.
Prepare Notebook#
We load the necessary libraries and set the notebook’s configuration.
import itertools
import warnings
import arviz as az
import arviz_plots as azp
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import pandas as pd
import pymc as pm
import pymc.dims as pmd
import pytensor
import pytensor.tensor as pt
import xarray as xr
from pydantic import InstanceOf
from pymc import do
from pymc_extras.prior import Prior
from pytensor.graph import rewrite_graph
from pymc_marketing.mmm import GeometricAdstock, LogisticSaturation
from pymc_marketing.mmm.additive_effect import DataVarMuEffect, IncrementalitySpec
from pymc_marketing.mmm.budget_optimizer import BudgetOptimizer
from pymc_marketing.mmm.media_transformation import MediaTransformation
from pymc_marketing.mmm.mmm import MMM
# A harmless Numba object-mode notice for the truncated-normal likelihood.
warnings.filterwarnings("ignore", message="Numba will use object mode")
az.style.use("arviz-darkgrid")
plt.rcParams["figure.figsize"] = [10, 6]
plt.rcParams["figure.dpi"] = 100
plt.rcParams["figure.facecolor"] = "white"
%config InlineBackend.figure_format = "retina"
The imports cell also silences a harmless warning that Numba runs the truncated-normal likelihood of lower_spend in object mode, which affects speed only.
Next, we fix the random seed so that the simulation is reproducible.
seed: int = sum(map(ord, "geo funnel model"))
rng: np.random.Generator = np.random.default_rng(seed=seed)
We simulate \(130\) weeks of data for three geos (north, south and west) and two upper-funnel channels. Every adstock uses a maximum lag of \(8\) weeks.
n_dates = 130
l_max = 8
geos = ["north", "south", "west"]
channels = ["tv_spend", "social_spend"]
date_range = pd.date_range(start="2021-01-04", freq="W-MON", periods=n_dates)
coords = {"date": date_range, "geo": geos}
Model Specification#
Funnel Effect#
We encode the funnel in a custom effect class, BudgetDualIndicatorEffect. The class name describes the subgraph topology, and “funnel” names the use case. Other topologies need their own effect class or composable building blocks.
This class is a notebook example, not a public library API.
A MuEffect adds a term to the target mean through three methods. create_data registers the data that the effect needs, create_effect builds and returns the term, and set_data updates the data at prediction time.
Because create_effect runs inside the model context, it can declare its own observed likelihoods. This is how one effect fits extra structural equations jointly with the target.
DataVarMuEffect removes most of this boilerplate, and both this notebook and Funnel-Aware MMM: Upper to Lower-Funnel Mediation via a custom MuEffect build on it. It takes a list of data_vars (variables in the training dataset) and implements create_data and set_data for us. Each variable keeps the dimensions that it has in the dataset, so the effect also works out of sample.
This notebook adds one more helper. MediaTransformation does for an adstock and saturation pair what DataVarMuEffect does for the data. The effect in the basic notebook holds four separate components, while this one holds two transformation objects.
A single adstock_first flag sets the order within each pair, and the priors of every component carry their own dimensions. We need this because every prior in the geo panel has dimensions.
The method incrementality_spec matters for the ROAS section. The effect reads mmm.channel_data_scaled, so a spend counterfactual reaches the target through the effect as well as through channel_contribution.
Returning an IncrementalitySpec declares this to the incrementality module. Without the declaration, the module raises an error instead of reporting the direct path as if it were the total.
The mediator chains a second adstock after the adstock of the base model, so a change in spend can move the target for extra weeks. The declaration records this extra reach, and the yearly evaluation windows below use it.
The method create_effect contains the structural equations of the funnel:
class BudgetDualIndicatorEffect(DataVarMuEffect):
"""Geo-level mediation through budget-linked lower-funnel spend.
Encodes ``upper channels -> demand -> lower-funnel spend -> target`` as an
additive effect on the target mean, while jointly fitting observed spend
and a second indicator of the unobserved demand hub.
"""
hub_input_transform: InstanceOf[MediaTransformation]
contribution_transform: InstanceOf[MediaTransformation]
model_config = {"arbitrary_types_allowed": True}
def to_dict(self) -> dict:
"""Serialize the effect (transformations delegate to their own dicts)."""
return {
"data_vars": self.data_vars,
"prefix": self.prefix,
"hub_input_transform": self.hub_input_transform.to_dict(),
"contribution_transform": self.contribution_transform.to_dict(),
}
def incrementality_spec(self) -> IncrementalitySpec:
"""Opt in to incrementality and declare the reach of the second adstock.
A spend counterfactual moves this contribution because the effect reads
``mmm.channel_data_scaled``. Opting in tells the incrementality module to
include that movement in the increment. The mediator chains a second
adstock after the adstock of the base model, so a change in spend can move
this contribution up to ``contribution_transform.adstock.l_max`` weeks
after the base ``l_max``. The declaration sizes the evaluation windows of
the yearly breakdown below. An empty ``IncrementalitySpec()`` also works:
the module then measures the reach on the fitted graph, as the basic
funnel notebook does.
"""
return IncrementalitySpec(
additional_carryover_lags=self.contribution_transform.adstock.l_max
)
def create_effect(self, mmm):
"""Build the budget-linked dual-indicator subgraph from the geo funnel DAG."""
model = mmm.model
# ("date", "geo", "channel") -> ("date", "geo"): all upper
# channels push on one demand pool, so sum away channel.
parents_into_hub = self.hub_input_transform(
mmm.channel_data_scaled, dim="date"
).sum(dim="channel")
# Structural equation 1: upper spend and category demand create D.
baseline = pmd.HalfNormal(f"{self.prefix}_baseline", sigma=1.0, dims=("geo",))
parent_coef = pmd.HalfNormal(
f"{self.prefix}_parent_coef", sigma=1.0, dims=("geo",)
)
hub = pmd.Deterministic(
f"{self.prefix}_hub",
baseline + parent_coef * model["category_demand"] + parents_into_hub,
)
# Structural equation 2: spend responds to D with unit coefficient
# (anchoring D's scale) and to its own exogenous budget.
budget_coef = pmd.HalfNormal(
f"{self.prefix}_budget_coef", sigma=0.5, dims=("geo",)
)
spend = pmd.Deterministic(
f"{self.prefix}_spend",
hub + budget_coef * model["lf_budget"],
)
# Measurement: reported lower-funnel spend, up to reporting noise.
pmd.TruncatedNormal(
f"{self.prefix}_spend_measurement",
mu=spend,
sigma=pmd.HalfNormal(
f"{self.prefix}_spend_measurement_sigma", sigma=1.0, dims=("geo",)
),
lower=0.0,
observed=model["lower_spend"],
)
# Indicator: branded search loads on demand, not on spend.
indicator_coef = pmd.HalfNormal(f"{self.prefix}_indicator_coef", sigma=1.0)
pmd.Normal(
f"{self.prefix}_hub_indicator",
mu=indicator_coef * hub,
sigma=pmd.HalfNormal(
f"{self.prefix}_hub_indicator_sigma", sigma=1.0, dims=("geo",)
),
observed=model["search_volume"],
)
# Structural equation 3: lower-funnel spend enters the target mean.
return pmd.Deterministic(
f"{self.prefix}_effect_contribution",
self.contribution_transform(spend, dim="date"),
)
Next, a factory function builds the two transformations and sets the dimensions of each prior from the pooling table above. GeometricAdstock and LogisticSaturation accept a priors dictionary.
In the demand equation, alpha and lam get dims="channel" (pooled across geos), while beta gets dims=("geo", "channel"). In the conversion equation, alpha and beta vary by geo and lam is a scalar.
Unique prefix values keep these parameters separate from the direct-path transformations of the base MMM.
def make_budget_dual_indicator_effect() -> BudgetDualIndicatorEffect:
"""Construct the funnel effect with uniquely-prefixed transformations."""
return BudgetDualIndicatorEffect(
data_vars=["lower_spend", "search_volume", "category_demand", "lf_budget"],
prefix="budget_dual",
hub_input_transform=MediaTransformation(
adstock=GeometricAdstock(
l_max=l_max,
prefix="adstock_hub_input",
priors={"alpha": Prior("Beta", alpha=2, beta=3, dims="channel")},
),
saturation=LogisticSaturation(
prefix="sat_hub_input",
priors={
"lam": Prior("Gamma", mu=2.0, sigma=1.0, dims="channel"),
"beta": Prior("HalfNormal", sigma=1.0, dims=("geo", "channel")),
},
),
adstock_first=True,
dims=("geo", "channel"),
),
contribution_transform=MediaTransformation(
adstock=GeometricAdstock(
l_max=l_max,
prefix="adstock_contribution",
priors={"alpha": Prior("Beta", alpha=2, beta=3, dims="geo")},
),
saturation=LogisticSaturation(
prefix="sat_contribution",
priors={
"lam": Prior("Gamma", mu=1.5, sigma=0.5),
"beta": Prior("HalfNormal", sigma=1.0, dims="geo"),
},
),
adstock_first=True,
dims=("geo",),
),
)
Data Plumbing#
DataVarMuEffect looks up its variables in the dataset that the MMM stores when it builds the model.
Warning
DataFrames drop the funnel columns: when X is a pandas.DataFrame, the conversion keeps only the channel and control columns, and every other column disappears without an error. To make extra variables visible to the effect, build the model from an xarray.Dataset, which keeps every data variable.
Also, MMM.fit needs a flat (date, *dims) frame for its bookkeeping, and a Dataset with a channel dimension cannot become one.
We build the model with the Dataset and fit it with the equivalent long DataFrame. The Dataset alone determines the model.
As a consequence, the stored fit_data does not contain the funnel columns, so a save and build_from_idata round trip would not rebuild this model.
Data Generating Process#
As in Generating Synthetic MMM Data, we build a model, clamp every parameter to a known value with the do() operator, and sample forward. The generative model is the funnel MMM, so one forward pass produces the target and both demand proxies at once.
Exogenous Covariates#
Three drivers (the two media plans and the category demand index) share an annual seasonal component, and two of them share a strong growth trend. The amplitudes and trends in the next cell define the confounding:
tv_spend: it follows the category cycle (amplitude \(0.80\), against \(0.70\) for category demand) and the same growth trend (loading \(1.20\), against \(1.00\) for the category). This setup mimics a growing category whose main channel grows with it.social_spend: it has a much milder cycle (\(0.40\)) and an almost flat trend (\(0.15\)), so it moves with category demand through the seasonal strand only.
This way, the two channels are confounded with demand through different strands of the shared driver, and this holds by construction for any random draw.
The noise around these shared means comes from an LKJ prior with eta=50 and scales close to \(0.3\). This concentration keeps the noise close to independent across drivers, so the raw correlations between plans and demand follow the deterministic means.
The residual correlations depend much more on the draw. With \(\eta = 50\) and \(d = 3\), the marginal correlation sd of the LKJ prior is \(1/\sqrt{2\eta + 2} \approx 0.10\), so a residual correlation of this size can change from seed to seed.
The ranking of the residuals does not depend on the draw, because the trend loading of tv_spend is \(1.20\) and that of social_spend is \(0.15\). We set these loadings, we do not draw them, and the causal argument below only needs them.
We still measure the realized noise, because part of the argument depends on its size. A softplus keeps every series non-negative, and a size multiplier scales each region.
t = np.arange(n_dates) / n_dates
week_of_year = date_range.isocalendar().week.to_numpy() / 52.0
season = np.sin(2 * np.pi * week_of_year)
drivers = ["tv_spend", "social_spend", "category_demand"]
geo_size = np.array([1.0, 0.7, 0.45])
level = np.array([0.10, 0.05, 0.30])
trend = np.array([1.20, 0.15, 1.00])
seasonal_amplitude = np.array([0.80, 0.40, 0.70])
cov_coords = {"date": date_range, "geo": geos, "driver": drivers}
with pm.Model(coords=cov_coords) as covariates_model:
L, _, _ = pm.LKJCholeskyCov(
"L", n=3, eta=50, sd_dist=pm.Gamma.dist(mu=0.30, sigma=0.05)
)
mu_cov = (
level + trend * t[:, None, None] + seasonal_amplitude * season[:, None, None]
)
x_raw = pm.MvNormal("x_raw", mu=mu_cov, chol=L, dims=("date", "geo", "driver"))
x = pm.Deterministic("x", pt.softplus(x_raw), dims=("date", "geo", "driver"))
Next, we draw one realization of the three drivers and scale each geo by its size. We also rescale the category index so that its maximum over all weeks and geos is one.
x_data = pm.draw(covariates_model.x, draws=1, random_seed=rng) * geo_size[None, :, None]
media_raw = x_data[..., :2]
category_demand = x_data[..., 2] / x_data[..., 2].max()
We generate the lower-funnel budget differently on purpose. It has a low always-on level with occasional promotional pushes on the calendar of the search team. It is not seasonal, and it does not follow the media plans.
The pushes are exogenous shocks to lower-funnel spend, so they let the model separate the budget-driven part of this spend from the demand-driven part. The table shows the mean and the standard deviation of the budget per geo.
promo = rng.binomial(n=1, p=0.35, size=(n_dates, len(geos)))
magnitude = rng.gamma(shape=2.0, scale=0.30, size=(n_dates, len(geos)))
lf_budget = geo_size[None, :] * (0.12 + promo * magnitude)
budget_stats = pd.DataFrame(
{
"mean": lf_budget.mean(axis=0),
"sd": lf_budget.std(axis=0),
"sd/mean": lf_budget.std(axis=0) / lf_budget.mean(axis=0),
},
index=pd.Index(geos, name="geo"),
)
budget_stats.style.format({"mean": "{:.3f}", "sd": "{:.3f}", "sd/mean": "{:.2f}"})
| mean | sd | sd/mean | |
|---|---|---|---|
| geo | |||
| north | 0.329 | 0.350 | 1.07 |
| south | 0.249 | 0.283 | 1.14 |
| west | 0.133 | 0.138 | 1.04 |
The standard deviation is \(4\%\) to \(14\%\) larger than the mean in every geo (sd/mean from \(1.04\) to \(1.14\)), so the budget varies a lot. This variation is necessary, because a budget that stayed at its always-on level would carry no information, and the demand baseline would absorb \(\lambda\).
Next, we check that the budget is exogenous. A promotional calendar that followed the category cycle would be a second confounder instead of a source of clean variation. The table shows the raw correlations between the four drivers in each geo:
driver_names = ["tv_spend", "social_spend", "category_demand", "lf_budget"]
corr_frames = {}
for gi, geo in enumerate(geos):
stacked = np.column_stack(
[
media_raw[:, gi, 0],
media_raw[:, gi, 1],
category_demand[:, gi],
lf_budget[:, gi],
]
)
corr_frames[geo] = pd.DataFrame(
np.corrcoef(stacked, rowvar=False), index=driver_names, columns=driver_names
)
corr_table = pd.concat(corr_frames, names=["geo", "driver"])
max_budget_corr = (
corr_table["lf_budget"].drop(index="lf_budget", level="driver").abs().max()
)
corr_table.style.background_gradient(cmap="RdBu_r", vmin=-1, vmax=1).format(
"{:+.2f}"
).set_caption(
f"Raw driver correlations; max |corr| involving lf_budget: {max_budget_corr:.3f}"
)
| tv_spend | social_spend | category_demand | lf_budget | ||
|---|---|---|---|---|---|
| geo | driver | ||||
| north | tv_spend | +1.00 | +0.46 | +0.84 | -0.05 |
| social_spend | +0.46 | +1.00 | +0.41 | -0.13 | |
| category_demand | +0.84 | +0.41 | +1.00 | -0.03 | |
| lf_budget | -0.05 | -0.13 | -0.03 | +1.00 | |
| south | tv_spend | +1.00 | +0.53 | +0.83 | +0.04 |
| social_spend | +0.53 | +1.00 | +0.53 | +0.11 | |
| category_demand | +0.83 | +0.53 | +1.00 | +0.07 | |
| lf_budget | +0.04 | +0.11 | +0.07 | +1.00 | |
| west | tv_spend | +1.00 | +0.53 | +0.86 | -0.06 |
| social_spend | +0.53 | +1.00 | +0.47 | -0.06 | |
| category_demand | +0.86 | +0.47 | +1.00 | -0.08 | |
| lf_budget | -0.06 | -0.06 | -0.08 | +1.00 |
The budget passes the check: its largest absolute correlation with any other driver is \(0.127\). In contrast, tv_spend tracks category demand much more closely than social_spend does (\(0.83\) to \(0.86\) against \(0.41\) to \(0.53\)), as we configured.
A natural guess is that the naive models will be biased mostly on tv_spend, the channel that moves with demand. We test this guess in the Media Deep Dive.
Also, a raw correlation mixes the strands of the shared driver. Every model below adjusts for the seasonal strand, and only two of them adjust for the trend strand. The next table separates the two strands.
Residual Co-Movement#
We quantify the split with partial correlations. The residual correlation after we project out the Fourier basis is the co-movement that the seasonal adjustment cannot absorb, so it is the part that can bias the models below.
The residual after we also project out \(t\) is what remains for a model with a trend control. The table caption also reports the \(R^2\) of the seasonal cycle of the drivers on the Fourier basis of the models.
doy_frac = date_range.dayofyear.to_numpy() / 365.25
fourier_basis = np.column_stack(
[f(2 * np.pi * k * doy_frac) for k in (1, 2) for f in (np.sin, np.cos)]
)
def partial_out(v: np.ndarray, controls: np.ndarray) -> np.ndarray:
"""Residual of ``v`` after an OLS projection on ``[1, controls]``."""
design = np.column_stack([np.ones(len(v)), controls])
return v - design @ np.linalg.lstsq(design, v, rcond=None)[0]
# The drivers' cycle and the models' basis are different functions of the
# calendar (isoweek/52 against dayofyear/365.25), so their overlap is a
# measurement, not an identity.
season_residual = partial_out(season, fourier_basis)
season_r2 = 1 - season_residual.var() / season.var()
residual_rows = []
for gi, geo in enumerate(geos):
demand_g = category_demand[:, gi]
for ci, ch in enumerate(channels):
spend_g = media_raw[:, gi, ci]
basis_t = np.column_stack([fourier_basis, t])
residual_rows.append(
{
"geo": geo,
"channel": ch,
"raw": np.corrcoef(spend_g, demand_g)[0, 1],
"after Fourier basis": np.corrcoef(
partial_out(spend_g, fourier_basis),
partial_out(demand_g, fourier_basis),
)[0, 1],
"after Fourier basis + t": np.corrcoef(
partial_out(spend_g, basis_t), partial_out(demand_g, basis_t)
)[0, 1],
}
)
residual_corr = pd.DataFrame(residual_rows).set_index(["geo", "channel"])
residual_corr.style.format("{:+.3f}").set_caption(
"Plan/demand correlations after partialing. "
f"R² of the drivers' seasonal cycle on the model's Fourier basis: {season_r2:.4f}"
)
| raw | after Fourier basis | after Fourier basis + t | ||
|---|---|---|---|---|
| geo | channel | |||
| north | tv_spend | +0.842 | +0.627 | +0.050 |
| social_spend | +0.406 | -0.053 | -0.206 | |
| south | tv_spend | +0.834 | +0.628 | +0.011 |
| social_spend | +0.527 | +0.119 | -0.030 | |
| west | tv_spend | +0.861 | +0.673 | +0.048 |
| social_spend | +0.473 | +0.045 | -0.122 |
The two channels behave very differently. The co-movement of tv_spend with category demand survives the seasonal projection (\(+0.63\) to \(+0.67\)), because the tv_spend plan follows the same growth trend as the category. No annual basis absorbs a trend.
The residual co-movement of social_spend is already at noise level (\(|r| \le 0.12\)), because its raw correlation with demand was almost entirely seasonal. Adding \(t\) to the projection reduces the residual of tv_spend to \(|r| \le 0.05\). Hence, the remaining confounding has the shape of a trend, and in this linear diagnostic a trend control removes it.
The social_spend column also moves, from \(-0.05\), \(+0.12\) and \(+0.05\) to \(-0.21\), \(-0.03\) and \(-0.12\). These negative residuals do not mean that the trend control over-corrects. They are already in the draw, as the next table shows.
Keep both results in mind for the model comparison. Both channels have a backdoor through the shared trend, and the d-separation analysis reports it for both. However, only the backdoor of tv_spend is large in numbers: \(+0.63\) to \(+0.67\) against \(|r| \le 0.12\).
Therefore, if a naive model gets social_spend wrong below, confounding is unlikely to be the main reason. D-separation gives a yes-or-no verdict, while this measurement gives a size, and we need both.
Note
Two seasonal functions: the season node in the DAG stands for two nearly collinear functions. The drivers follow \(\sin(2\pi \cdot \text{isoweek}/52)\), while every model uses a two-mode Fourier basis in dayofyear/365.25. The \(R^2\) of \(0.9998\) in the table caption confirms that the Fourier basis absorbs the seasonal cycle of the drivers.
Ground-Truth Parameters#
We give every structural parameter an explicit true value with the shape that its dimensions imply. Instead of relying on memory for which parameters are pooled, we build the model first and read the dimensions off it.
true = {
# target equation: direct paths and seasonality
"intercept_contribution": np.array([0.35, 0.28, 0.22]),
"adstock_alpha": np.array([0.55, 0.30]),
"saturation_lam": np.array([3.0, 4.0]),
"saturation_beta": np.array([[0.55, 0.35], [0.45, 0.30], [0.35, 0.25]]),
"gamma_fourier": np.array(
[
[0.05, 0.03, 0.02, -0.02],
[0.04, 0.02, 0.02, -0.01],
[0.03, 0.02, 0.01, -0.01],
]
),
"y_sigma": np.array([0.04, 0.04, 0.04]),
# mediator equation: upper spend + category demand -> latent demand
"budget_dual_baseline": np.array([0.20, 0.16, 0.12]),
"budget_dual_parent_coef": np.array([0.35, 0.30, 0.25]),
"adstock_hub_input_alpha": np.array([0.60, 0.35]),
"sat_hub_input_lam": np.array([2.5, 3.0]),
"sat_hub_input_beta": np.array([[0.85, 0.55], [0.70, 0.45], [0.55, 0.35]]),
# lower-funnel spend equation and its observation
"budget_dual_budget_coef": np.array([0.50, 0.45, 0.40]),
"budget_dual_spend_measurement_sigma": np.array([0.05, 0.05, 0.05]),
# branded-search indicator
"budget_dual_indicator_coef": np.array(0.8),
"budget_dual_hub_indicator_sigma": np.array([0.05, 0.05, 0.05]),
# conversion path: lower-funnel spend -> target
"adstock_contribution_alpha": np.array([0.35, 0.30, 0.25]),
"sat_contribution_lam": np.array(1.5),
"sat_contribution_beta": np.array([1.05, 0.90, 0.72]),
}
Model Configuration#
The generative model and the fitted models use the same priors, so we define them once. Note the scaling block: dims=() means “take the maximum over the date dimension only”, which gives one scale per geo and per channel. This is the usual choice for a geo panel where regions differ in size.
The helper make_mmm builds the base model that every specification shares, and make_dataset assembles the xarray.Dataset that the funnel model needs.
SCALING = {
"channel": {"method": "max", "dims": ()},
"target": {"method": "max", "dims": ()},
}
MODEL_CONFIG = {
"intercept": Prior("Normal", mu=0.2, sigma=0.2, dims="geo"),
"gamma_fourier": Prior("Normal", mu=0, sigma=0.1, dims=("geo", "fourier_mode")),
"likelihood": Prior(
"Normal",
sigma=Prior("HalfNormal", sigma=0.5, dims="geo"),
dims=("date", "geo"),
),
}
def make_mmm(channel_columns: list[str], **kwargs) -> MMM:
"""Build the base MMM shared by the generative and fitted models."""
return MMM(
date_column="date",
target_column="y",
channel_columns=channel_columns,
dims=("geo",),
scaling=SCALING,
adstock=GeometricAdstock(
l_max=l_max,
priors={"alpha": Prior("Beta", alpha=2, beta=3, dims="channel")},
),
saturation=LogisticSaturation(
priors={
"lam": Prior("Gamma", mu=3.0, sigma=1.0, dims="channel"),
"beta": Prior("HalfNormal", sigma=1.0, dims=("geo", "channel")),
}
),
yearly_seasonality=2,
model_config=MODEL_CONFIG,
**kwargs,
)
def make_dataset(media, lower_spend, search_volume, cat_demand, budget) -> xr.Dataset:
"""Assemble the xarray Dataset that the funnel model is built from.
The media variable must carry a ``channel`` dimension and must be named
``media``, not ``channel``. With the explicit ``channel`` coordinate below,
xarray rejects a data variable of the same name with a ``ValueError``
(found in both data_vars and coords).
Every variable in the ``data_vars`` of the effect must appear here.
Otherwise ``build_model`` raises the ``KeyError`` of xarray ("No variable
named ..."), because the MMM does not validate these variables itself.
"""
return xr.Dataset(
{
"media": xr.DataArray(media, dims=("date", "geo", "channel")),
"lower_spend": xr.DataArray(lower_spend, dims=("date", "geo")),
"search_volume": xr.DataArray(search_volume, dims=("date", "geo")),
"category_demand": xr.DataArray(cat_demand, dims=("date", "geo")),
"lf_budget": xr.DataArray(budget, dims=("date", "geo")),
},
coords={**coords, "channel": channels},
)
Forward Simulation#
We build the generative model on placeholder observations and a placeholder target of ones, so the target scale is one and the true parameters are on the natural scale. The table lists the dimensions of every model variable next to the shape of its true value.
Warning
Placeholders and real inputs: lower_spend and search_volume are placeholders here, because they are outputs of the forward pass. lf_budget is not a placeholder. It is an input to the mediator equation and must be the real series.
With ones in its place, \(\lambda\) would multiply a constant. Sampling would still succeed, and the posterior of \(\lambda\) would return its prior without raising an error.
ones = np.ones((n_dates, len(geos)))
ds_gen = make_dataset(media_raw, ones, ones, category_demand, lf_budget)
y_gen = xr.DataArray(ones, dims=("date", "geo"), coords=coords)
gen = make_mmm(channels)
gen.add_mu_effect(make_budget_dual_indicator_effect())
gen.build_model(X=ds_gen, y=y_gen)
gen.add_original_scale_contribution_variable(
var=["channel_contribution", "budget_dual_effect_contribution", "y"]
)
pd.DataFrame(
[(v, gen.model.named_vars_to_dims.get(v), true[v].shape) for v in true],
columns=["parameter", "dims", "true value shape"],
)
| parameter | dims | true value shape | |
|---|---|---|---|
| 0 | intercept_contribution | (geo,) | (3,) |
| 1 | adstock_alpha | (channel,) | (2,) |
| 2 | saturation_lam | (channel,) | (2,) |
| 3 | saturation_beta | (geo, channel) | (3, 2) |
| 4 | gamma_fourier | (geo, fourier_mode) | (3, 4) |
| 5 | y_sigma | (geo,) | (3,) |
| 6 | budget_dual_baseline | (geo,) | (3,) |
| 7 | budget_dual_parent_coef | (geo,) | (3,) |
| 8 | adstock_hub_input_alpha | (channel,) | (2,) |
| 9 | sat_hub_input_lam | (channel,) | (2,) |
| 10 | sat_hub_input_beta | (geo, channel) | (3, 2) |
| 11 | budget_dual_budget_coef | (geo,) | (3,) |
| 12 | budget_dual_spend_measurement_sigma | (geo,) | (3,) |
| 13 | budget_dual_indicator_coef | () | () |
| 14 | budget_dual_hub_indicator_sigma | (geo,) | (3,) |
| 15 | adstock_contribution_alpha | (geo,) | (3,) |
| 16 | sat_contribution_lam | () | () |
| 17 | sat_contribution_beta | (geo,) | (3,) |
Every true value has the shape of its model variable. For example, saturation_beta has dimensions (geo, channel) and a true value of shape (3, 2), while budget_dual_indicator_coef and sat_contribution_lam are scalars.
The method table() gives a summary of the whole model, with every variable, its dimensions and its prior in one place:
gen.table()
Variable Expression Dimensions ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────── channel_scale = Data geo[3] × channel[2] target_scale = Data geo[3] channel_data = Data date[130] × geo[3] × channel[2] target_data = Data date[130] × geo[3] lower_spend = Data date[130] × geo[3] search_volume = Data date[130] × geo[3] category_demand = Data date[130] × geo[3] lf_budget = Data date[130] × geo[3] dayofyear = Data date[130] intercept_contribution ~ Normal(0.2, 0.2) geo[3] adstock_alpha ~ Beta(2, 3) channel[2] saturation_lam ~ Gamma(<constant>, <constant>) channel[2] saturation_beta ~ HalfNormal(0, 1) geo[3] × channel[2] gamma_fourier ~ Normal(0, 0.1) geo[3] × fourier_mode[4] adstock_hub_input_alpha ~ Beta(2, 3) channel[2] sat_hub_input_lam ~ Gamma(<constant>, <constant>) channel[2] sat_hub_input_beta ~ HalfNormal(0, 1) geo[3] × channel[2] budget_dual_baseline ~ HalfNormal(0, 1) geo[3] budget_dual_parent_coef ~ HalfNormal(0, 1) geo[3] budget_dual_budget_coef ~ HalfNormal(0, 0.5) geo[3] budget_dual_spend_measurement_sigma HalfNormal(0, 1) geo[3] ~ budget_dual_indicator_coef ~ HalfNormal(0, 1) budget_dual_hub_indicator_sigma ~ HalfNormal(0, 1) geo[3] adstock_contribution_alpha ~ Beta(2, 3) geo[3] sat_contribution_lam ~ Gamma(<constant>, <constant>) sat_contribution_beta ~ HalfNormal(0, 1) geo[3] y_sigma ~ HalfNormal(0, 0.5) geo[3] Parameter count = 61 channel_contribution = f(saturation_beta, saturation_lam, date[130] × geo[3] × channel[2] adstock_alpha) fourier_contribution = f(gamma_fourier) date[130] × geo[3] × fourier_mode[4] yearly_seasonality_contribution = f(gamma_fourier) date[130] × geo[3] budget_dual_hub = f(budget_dual_baseline, geo[3] × date[130] sat_hub_input_beta, budget_dual_parent_coef, sat_hub_input_lam, adstock_hub_input_alpha) budget_dual_spend = f(budget_dual_budget_coef, geo[3] × date[130] budget_dual_baseline, sat_hub_input_beta, budget_dual_parent_coef, sat_hub_input_lam, adstock_hub_input_alpha) budget_dual_effect_contribution = f(sat_contribution_beta, geo[3] × date[130] sat_contribution_lam, budget_dual_budget_coef, adstock_contribution_alpha, budget_dual_baseline, sat_hub_input_beta, budget_dual_parent_coef, sat_hub_input_lam, adstock_hub_input_alpha) mu = f(sat_contribution_beta, geo[3] × date[130] intercept_contribution, gamma_fourier, saturation_beta, sat_contribution_lam, budget_dual_budget_coef, saturation_lam, adstock_contribution_alpha, budget_dual_baseline, sat_hub_input_beta, budget_dual_parent_coef, adstock_alpha, sat_hub_input_lam, adstock_hub_input_alpha) total_media_contribution_original_s… f(saturation_beta, saturation_lam, = adstock_alpha) total_response_original_scale = f(sat_contribution_beta, intercept_contribution, gamma_fourier, saturation_beta, sat_contribution_lam, budget_dual_budget_coef, saturation_lam, adstock_contribution_alpha, budget_dual_baseline, sat_hub_input_beta, budget_dual_parent_coef, adstock_alpha, sat_hub_input_lam, adstock_hub_input_alpha) channel_contribution_original_scale f(saturation_beta, saturation_lam, date[130] × geo[3] × channel[2] = adstock_alpha) budget_dual_effect_contribution_ori… f(sat_contribution_beta, date[130] × geo[3] = sat_contribution_lam, budget_dual_budget_coef, adstock_contribution_alpha, budget_dual_baseline, sat_hub_input_beta, budget_dual_parent_coef, sat_hub_input_lam, adstock_hub_input_alpha) y_original_scale = f(<normal>) date[130] × geo[3] budget_dual_spend_measurement ~ TruncatedNormal(f(budget_dual_budget… geo[3] × date[130] budget_dual_baseline, sat_hub_input_beta, budget_dual_parent_coef, sat_hub_input_lam, adstock_hub_input_alpha), budget_dual_spend_measurement_sigma, 0, inf) budget_dual_hub_indicator ~ Normal(f(budget_dual_indicator_coef, geo[3] × date[130] budget_dual_baseline, sat_hub_input_beta, budget_dual_parent_coef, sat_hub_input_lam, adstock_hub_input_alpha), budget_dual_hub_indicator_sigma) y ~ Unknown(Normal(f(sat_contribution_be… date[130] × geo[3] intercept_contribution, gamma_fourier, saturation_beta, sat_contribution_lam, budget_dual_budget_coef, saturation_lam, adstock_contribution_alpha, budget_dual_baseline, sat_hub_input_beta, budget_dual_parent_coef, adstock_alpha, sat_hub_input_lam, adstock_hub_input_alpha), y_sigma))
The summary lists \(61\) parameters. The data block contains the four funnel series (lower_spend, search_volume, category_demand and lf_budget) next to the channel and target data.
We read the dimensions off the model before we write the intervention, because a clamped value of the wrong length can pass do without raising an error. The mistake only appears at the first predictive draw, two steps away from its cause. With the pooling choices above, the shapes are not all obvious.
The model graph shows the three likelihoods. The target and the two demand proxies all depend on the same set of latent parameters.
pm.model_to_graphviz(gen.model)
The three observed nodes are y, budget_dual_spend_measurement and budget_dual_hub_indicator. Next, we clamp every parameter to its true value with do and draw one forward sample of the target, the two proxies and the true contributions.
gen.model = do(gen.model, true)
var_names = [
"y_original_scale",
"channel_contribution_original_scale",
"budget_dual_effect_contribution_original_scale",
"budget_dual_hub",
"budget_dual_spend",
]
pp_names = ["budget_dual_spend_measurement", "budget_dual_hub_indicator"]
with gen.model:
idata_gen = pm.sample_prior_predictive(
draws=1, var_names=var_names + pp_names, random_seed=rng
)
prior = idata_gen["prior"].sel(chain=0, draw=0)
prior_pp = idata_gen["prior_predictive"].sel(chain=0, draw=0)
y_obs = prior["y_original_scale"].transpose("date", "geo").to_numpy()
demand_true = prior["budget_dual_hub"].transpose("date", "geo").to_numpy()
lf_spend_true = prior["budget_dual_spend"].transpose("date", "geo").to_numpy()
lower_spend_obs = (
prior_pp["budget_dual_spend_measurement"].transpose("date", "geo").to_numpy()
)
search_obs = prior_pp["budget_dual_hub_indicator"].transpose("date", "geo").to_numpy()
direct_true = (
prior["channel_contribution_original_scale"]
.transpose("date", "geo", "channel")
.to_numpy()
)
indirect_raw = (
prior["budget_dual_effect_contribution_original_scale"]
.transpose("date", "geo")
.to_numpy()
)
True Causal Decomposition#
The raw budget_dual_effect_contribution includes the demand from the baseline and the category index, which upper-funnel spend never created. To isolate the causal indirect contribution of each channel, we run the clamped generative model again with the spend of that channel set to zero, and we take the difference.
The same counterfactual with both channels off measures the non-additivity.
This is the standard mixed-regime construction of mediation analysis: an expectation in which the treatment takes one value in the outcome equation and another value in the mediator equation. Mediation Analysis and (In)Direct Effects with PyMC writes these expectations as \(E_{t, t', t''}\) and evaluates them by intervening on the generative model.
The next cell does the same. The do operator above clamped the parameters to their true values, and pm.set_data swaps the inputs. The two regimes are the observed spend of a channel and zero spend.
def forward_with_zeroed(channel_mask) -> dict[str, np.ndarray]:
"""Mediated contribution and structural LF spend with masked channels zeroed."""
media_cf = media_raw.copy()
media_cf[..., channel_mask] = 0.0
# The generative model is restored whatever happens in between, so a raise
# here cannot leave every later truth number on the perturbed media.
with gen.model:
try:
pm.set_data({"channel_data": media_cf})
idata_cf = pm.sample_prior_predictive(
draws=1,
var_names=[
"budget_dual_effect_contribution_original_scale",
"budget_dual_spend",
],
random_seed=rng,
)
finally:
pm.set_data({"channel_data": media_raw})
draw = idata_cf["prior"].sel(chain=0, draw=0)
return {
var: draw[var].transpose("date", "geo").to_numpy()
for var in [
"budget_dual_effect_contribution_original_scale",
"budget_dual_spend",
]
}
cf_gen = {
"tv_spend": forward_with_zeroed(np.array([True, False])),
"social_spend": forward_with_zeroed(np.array([False, True])),
"both": forward_with_zeroed(np.array([True, True])),
}
indirect_true = np.stack(
[
indirect_raw - cf_gen[ch]["budget_dual_effect_contribution_original_scale"]
for ch in channels
],
axis=-1,
)
indirect_both = (
indirect_raw - cf_gen["both"]["budget_dual_effect_contribution_original_scale"]
)
# The same counterfactual also prices each channel's *induced* lower-funnel spend:
# the structural spend M* with the channel on, minus M* with it off. Because the
# coefficient on D in the M* equation is exactly one, a unit of induced demand is a
# unit of incremental lower-funnel spend; we will need this when defining ROAS.
induced_spend_true = np.stack(
[lf_spend_true - cf_gen[ch]["budget_dual_spend"] for ch in channels], axis=-1
)
# Each counterfactual pass restores the generative model's media data itself.
np.testing.assert_allclose(gen.model["channel_data"].get_value(), media_raw)
total_true = direct_true + indirect_true
Note
Full date range: unlike the single-channel example, we keep all weeks and do not drop the first l_max weeks. The generative model and the fitted models apply the same zero-padded adstock, so the early weeks are not an artifact.
Dropping them would introduce one, because the fitted model would have to explain weeks whose carryover history we had deleted.
Next, we collect the data and the true decompositions in one xarray.Dataset. Then we compute the share of the total effect that travels through the mediator, for every geo and channel:
data = xr.Dataset(
{
"tv_spend": (("date", "geo"), media_raw[..., 0]),
"social_spend": (("date", "geo"), media_raw[..., 1]),
"category_demand": (("date", "geo"), category_demand),
"lf_budget": (("date", "geo"), lf_budget),
"lower_spend": (("date", "geo"), lower_spend_obs),
"search_volume": (("date", "geo"), search_obs),
"y": (("date", "geo"), y_obs),
"demand_true": (("date", "geo"), demand_true),
"lf_spend_true": (("date", "geo"), lf_spend_true),
"direct_true": (("date", "geo", "channel"), direct_true),
"indirect_true": (("date", "geo", "channel"), indirect_true),
"total_true": (("date", "geo", "channel"), total_true),
"induced_spend_true": (("date", "geo", "channel"), induced_spend_true),
},
coords={**coords, "channel": channels},
)
shares = (
data["indirect_true"].mean("date") / data["total_true"].mean("date")
).to_dataframe(name="indirect share")
shares["indirect share"] = shares["indirect share"].map("{:.1%}".format)
shares
| indirect share | ||
|---|---|---|
| geo | channel | |
| north | tv_spend | 34.9% |
| social_spend | 35.4% | |
| south | tv_spend | 37.7% |
| social_spend | 36.1% | |
| west | tv_spend | 37.2% |
| social_spend | 33.7% |
Between \(33.7\%\) and \(37.7\%\) of the effect of every upper-funnel channel travels through the mediator, a bit over a third. Next, we measure the non-additivity:
sum_of_singles = float(data["indirect_true"].sum("channel").mean())
joint = float(indirect_both.mean())
pd.DataFrame(
{
"sum of single-channel indirect effects": [sum_of_singles],
"indirect effect of switching both off": [joint],
"shortfall vs additive": [1 - sum_of_singles / joint],
},
index=["mean per week"],
).style.format(
{
"sum of single-channel indirect effects": "{:.3f}",
"indirect effect of switching both off": "{:.3f}",
"shortfall vs additive": "{:.1%}",
}
)
| sum of single-channel indirect effects | indirect effect of switching both off | shortfall vs additive | |
|---|---|---|---|
| mean per week | 0.282 | 0.326 | 13.3% |
The two channels share one saturating demand pool. Therefore, the indirect effects computed one channel at a time understate their joint effect by \(13.3\%\) (a mean of \(0.282\) per week against \(0.326\)).
This is expected. With a concave response, there is no unique way to split a joint effect into per-channel parts, and leave-one-out attribution is one convention among several. State the convention whenever these numbers enter a budget discussion.
The mediation literature gives this remainder a name. Mediation Analysis and (In)Direct Effects with PyMC decomposes a total effect with two mediators as \(\text{TE} = \text{DE} + \text{IIE}_1 + \text{IIE}_2 + \text{INT} + \text{DEP}\). The interaction and dependence terms are the part that the individual indirect effects do not add up to.
Our setting differs, because two treatments push on one shared mediator, while that decomposition has one treatment that reaches the outcome through two mediators. The lesson still applies: the gap is a named term, and you should tell your reader which convention you report.
Exploratory Data Analysis#
Let’s look at the synthetic data. The figure shows the target, the two upper-funnel channels, the two exogenous inputs, the mediator and the indicator, with one line per geo.
All series except lf_budget share the annual cycle, and tv_spend, category demand and the target also grow over the two and a half years. In contrast, lf_budget shows short promotional spikes without a seasonal pattern. north is the largest market and west the smallest in most series.
Business Context: upper-funnel plans, lower-funnel spend and branded search rise and fall together in this data. A direct comparison of these series cannot tell which of them drives sales.
Next, we zoom in on the latent layer of each geo: latent demand \(D\), structural lower-funnel spend \(M^{*}\), its observed version lower_spend and branded search \(S\).
In every geo, structural lower-funnel spend \(M^{*}\) sits above latent demand \(D\) by the budget wedge \(\lambda_g B_{t,g}\), with spikes in the promotional weeks. The observed lower_spend follows \(M^{*}\) up to reporting noise. Branded search follows the shape of \(D\) at a lower level, because \(\kappa = 0.8\).
D-Separation Analysis#
Before we fit anything, we can query the DAG itself. In causal terms, each model below is a conditioning set: the variables that it regresses the target on. Pearl’s backdoor criterion says that a set \(Z\) identifies the total effect of a treatment on \(Y\) when two conditions hold:
No-descendant condition: \(Z\) contains no descendant of the treatment.
Separation condition: \(Z\) d-separates the treatment from the target in the graph without the outgoing edges of the treatment.
Both conditions are mechanical graph checks, and networkx implements them. This means that we can compute which conditioning sets are valid for the total effect instead of arguing about it.
See also
If backdoor paths and d-separation are new to you, Pearl, Glymour and Jewell, Causal Inference in Statistics: A Primer (Wiley, 2016) is the standard introduction to both. The Book of Why (Pearl and Mackenzie, 2018) covers the ideas without the formalism.
For how these checks drive model specification in PyMC-Marketing, see Introduction to Understanding Causal Relationships in Media Mix Modeling and Introducing causal discovery to PyMC-Marketing.
Note that the graph is a single time slice. Both adstocks create lagged edges, \(U_{t-k} \to D_t\) and \(M^{*}_{t-k} \to Y_t\), and season and t change with the date. As a result, each node stands for a whole series, not for the value of one week.
This reading fits an MMM, whose control columns condition on the whole series. However, the figure does not show it, and a confounder that acts only with a lag would need a graph with explicit time.
Backdoor Verdicts#
We build the graph from the same edge list that drew the figure. We check every conditioning set with tv_spend (\(U_1\)) as the treatment, and the same loop confirms that each verdict also holds with social_spend (\(U_2\)) as the treatment:
G = nx.DiGraph([(parent, child) for parent, child, _ in DAG_EDGES])
treatment, outcome = "U1", "Y"
descendants = nx.descendants(G, treatment)
# Backdoor surgery: remove the treatment's outgoing edges.
G_do = G.copy()
G_do.remove_edges_from(list(G.out_edges(treatment)))
# What each model conditions on when estimating tv_spend's effect. Every model
# includes the other channel and the Fourier basis (= the `season` node); they
# differ only in what else they add.
conditioning_sets = {
"Naive A (mediator as channel)": {"U2", "season", "M"},
"Naive B (mediator omitted)": {"U2", "season"},
"Naive B+ (demand as control)": {"U2", "season", "C"},
"Naive B++ (trend as control)": {"U2", "season", "t"},
"Naive C (indicator as channel)": {"U2", "season", "S"},
"(hypothetical: condition on M*)": {"U2", "season", "MS"},
}
rows = {}
for label, Z in conditioning_sets.items():
forbidden = sorted(Z & descendants)
blocks = nx.is_d_separator(G_do, {treatment}, {outcome}, Z)
rows[label] = {
"conditions on": ", ".join(sorted(Z)),
"descendant of U₁ in Z": ", ".join(forbidden) or "none",
"blocks all backdoors": blocks,
"valid for total effect": blocks and not forbidden,
}
# The verdicts are symmetric in the treatment: re-run every check with U2 as
# treatment (its conditioning sets swap in the *other* channel, U1).
G_do_u2 = G.copy()
G_do_u2.remove_edges_from(list(G.out_edges("U2")))
descendants_u2 = nx.descendants(G, "U2")
for label, Z in conditioning_sets.items():
Z_u2 = {("U1" if z == "U2" else z) for z in Z}
blocks_u2 = nx.is_d_separator(G_do_u2, {"U2"}, {outcome}, Z_u2)
valid_u2 = blocks_u2 and not (Z_u2 & descendants_u2)
if (
blocks_u2 != rows[label]["blocks all backdoors"]
or valid_u2 != rows[label]["valid for total effect"]
):
raise ValueError(f"verdict differs with U2 as treatment: {label}")
pd.DataFrame(rows).T.style.set_caption(
f"Backdoor verdicts for treatment {treatment}. "
f"descendants({treatment}) = {{{', '.join(sorted(descendants))}}}. "
"Every verdict also holds with U2 as the treatment (checked)"
)
| conditions on | descendant of U₁ in Z | blocks all backdoors | valid for total effect | |
|---|---|---|---|---|
| Naive A (mediator as channel) | M, U2, season | M | False | False |
| Naive B (mediator omitted) | U2, season | none | False | False |
| Naive B+ (demand as control) | C, U2, season | none | True | True |
| Naive B++ (trend as control) | U2, season, t | none | False | False |
| Naive C (indicator as channel) | S, U2, season | S | False | False |
| (hypothetical: condition on M*) | MS, U2, season | MS | True | False |
Only Naive B+ blocks every backdoor path and is valid for the total effect. Naive A and Naive C condition on a descendant of the treatment (M and S), and Naive B and Naive B++ leave backdoor paths open.
Open Backdoor Paths#
A verdict is a yes-or-no answer, so it cannot say how a set fails or which paths it leaves open. A count of the open backdoor paths shows this. The helper below counts them, and we check every count against nx.is_d_separator, so that the helper and the library agree.
backdoor_paths = list(nx.all_simple_paths(G_do.to_undirected(), treatment, outcome))
def open_paths(Z):
"""Backdoor paths that conditioning set ``Z`` leaves open."""
left_open = []
for path in backdoor_paths:
blocked = False
for before, node, after in zip(path, path[1:], path[2:], strict=False):
if G_do.has_edge(before, node) and G_do.has_edge(after, node):
# A collider blocks unless it, or a descendant, is in Z.
blocked = not ({node} | nx.descendants(G_do, node)) & Z
else:
blocked = node in Z
if blocked:
break
if not blocked:
left_open.append(path)
return left_open
def render(path):
"""``U1 <- N -> C -> D -> MS -> Y``: arrows point the way the graph does."""
out = [path[0]]
for before, node in itertools.pairwise(path):
out += ["->" if G_do.has_edge(before, node) else "<-", node]
return " ".join(out)
adjustment_sets = {
"the seasonal basis alone": {"season"},
"Naive B": {"U2", "season"},
"Naive B+": {"U2", "season", "C"},
"Naive B++": {"U2", "season", "t"},
"Naive B+ without the Fourier basis": {"U2", "C"},
}
# Check the helper against the library: a set leaves zero paths open exactly
# when `is_d_separator` says it blocks them all.
for Z in adjustment_sets.values():
if (len(open_paths(Z)) == 0) != nx.is_d_separator(G_do, {treatment}, {outcome}, Z):
raise ValueError(f"open_paths disagrees with is_d_separator for {sorted(Z)}")
pd.DataFrame(
{
"conditioning set": {
label: "{" + ", ".join(sorted(Z)) + "}"
for label, Z in adjustment_sets.items()
},
"open backdoor paths": {
label: len(open_paths(Z)) for label, Z in adjustment_sets.items()
},
}
).style.set_caption(
f"Of {len(backdoor_paths)} backdoor paths from {treatment} to {outcome}"
)
| conditioning set | open backdoor paths | |
|---|---|---|
| the seasonal basis alone | {season} | 6 |
| Naive B | {U2, season} | 4 |
| Naive B+ | {C, U2, season} | 0 |
| Naive B++ | {U2, season, t} | 1 |
| Naive B+ without the Fourier basis | {C, U2} | 9 |
Of the \(43\) backdoor paths from \(U_1\) to \(Y\), the seasonal basis alone leaves \(6\) open, Naive B leaves \(4\), Naive B++ leaves \(1\) and Naive B+ leaves none. Without the Fourier basis, the conditioning set of Naive B+ leaves \(9\) paths open.
Tip
Colliders: conditioning on a collider opens paths that something else must close again. The last row (Naive B+ without its Fourier basis) shows this effect at C, where the seasonal basis closes the paths.
Next, we list the open paths of every set that leaves between one and four paths open:
pd.DataFrame(
[
{"conditioning set": label, "open path": render(path)}
for label, Z in adjustment_sets.items()
if 1 <= len(open_paths(Z)) <= 4
for path in open_paths(Z)
]
).set_index("conditioning set")
| open path | |
|---|---|
| conditioning set | |
| Naive B | U1 <- t -> U2 <- N -> C -> D -> MS -> Y |
| Naive B | U1 <- t -> C -> D -> MS -> Y |
| Naive B | U1 <- N -> U2 <- t -> C -> D -> MS -> Y |
| Naive B | U1 <- N -> C -> D -> MS -> Y |
| Naive B++ | U1 <- N -> C -> D -> MS -> Y |
The listing shows the four open paths of Naive B and the single open path of Naive B++, the noise fork \(U_1 \leftarrow N \to C \to D \to M^{*} \to Y\). Let’s read the verdicts and the paths together.
Naive B fails the separation condition. With only the other channel and the seasonal basis in the conditioning set, the backdoor \(U_1 \leftarrow t \to C \to D \to M^{*} \to Y\) stays open, and so does its twin through the shared noise.
Omitting the mediator is the correct choice for a total effect, so on the graph the only problem of Naive B is the open backdoor. Its trend strand is the co-movement of tv_spend that survived the seasonal adjustment above.
Naive B+ is a valid adjustment set, and Naive B++ is not. Both fork nodes, \(t\) and \(N\), point at the treatment and at \(C\), and \(C\) lies on both paths after them. Therefore, conditioning on the observed category index blocks the trend strand and the noise strand at once.
Conditioning on \(t\) blocks only the trend strand. The noise path stays open, and \(N\) is unobservable, so Naive B++ cannot close it with another column. Both sets meet the no-descendant condition, so the difference is which paths each one blocks.
Of the four open paths of Naive B, two are the forks themselves. The other two enter on one fork and leave on the other, joined at the other channel. season, \(t\) and \(N\) all point into \(U_2\), so the other channel is a collider as well as a regressor.
Conditioning on \(U_2\) closes the paths that run through it and opens the paths that meet at it. The count table shows this: the seasonal basis alone leaves six paths open, and adding \(U_2\) leaves four. All four leave through \(C \to D \to M^{*} \to Y\), so conditioning on \(C\) closes all of them.
The validity of Naive B+ depends on the Fourier basis. \(C\) is also a collider, because season, \(t\) and \(N\) point into it as they do into \(U_2\). Hence, conditioning on \(C\) opens every path that meets there, and the seasonal basis, which every model includes, closes them again.
Without the basis, the same conditioning set leaves nine paths open, as the last row of the count table shows. With yearly_seasonality switched off, Naive B+ is not a valid adjustment set.
Naive A and Naive C fail the no-descendant condition. Each conditions on a descendant of the treatment (M through \(U_1 \to D \to M^{*} \to M\), and S through \(U_1 \to D \to S\)), and conditioning on a proxy of a node on the causal path partly blocks that path. Extra controls cannot repair this.
Descendants of \(D\) also open collider paths at \(D\). On this graph, \(U_2\) blocks every such path again, and both conditioning sets contain \(U_2\). This protection comes from this particular graph, not from the conditioning choice.
Finally, conditioning on the structural \(M^{*}\) blocks every backdoor, but it is still invalid for the total effect, because \(M^{*}\) is a descendant of the treatment. Conditioning on it defines the direct effect instead.
Direct Effect Check#
We can make this last point precise. If we remove the direct edge \(U_1 \to Y\) and the conditioning set still separates the treatment from the target, then the set isolates the direct path:
G_direct = G.copy()
G_direct.remove_edge(treatment, "Y")
pd.DataFrame(
{
"isolates the direct effect": {
"M* (structural)": nx.is_d_separator(
G_direct, {treatment}, {outcome}, {"U2", "season", "MS"}
),
"M (noisy observation, what Naive A has)": nx.is_d_separator(
G_direct, {treatment}, {outcome}, {"U2", "season", "M"}
),
}
}
)
| isolates the direct effect | |
|---|---|
| M* (structural) | True |
| M (noisy observation, what Naive A has) | False |
The structural \(M^{*}\) isolates the direct effect, but the noisy observation M that Naive A uses does not. Hence, Naive A estimates the direct effect approximately, and the approximation is good only as far as the measurement noise \(\sigma^M\) is small.
Keep this in mind for the comparisons below, where we score Naive A against the direct effect as well as against the total effect.
Minimal Adjustment Sets#
Next, we ask the graph which adjustment set it recommends. find_minimal_d_separator returns one minimal separator, first without restrictions and then restricted to observed variables that are not descendants of the treatment:
# Sorted: a set's repr order is not stable across processes, and these outputs
# are meant to reproduce exactly.
unrestricted = sorted(nx.find_minimal_d_separator(G_do, {treatment}, {outcome}))
# `N` is deliberately absent from `observable`: no dataset carries the drivers'
# shared noise. A practitioner can only act on the restricted answer.
observable = {"U2", "season", "C", "B", "M", "S", "t"} - descendants
restricted = sorted(
nx.find_minimal_d_separator(
G_do,
{treatment},
{outcome},
included={"U2", "season"},
restricted=observable | {"U2", "season"},
)
)
pd.DataFrame(
{
"a minimal separator": {
"unrestricted": str(unrestricted),
"observables only": str(restricted),
}
}
)
| a minimal separator | |
|---|---|
| unrestricted | ['N', 'season', 't'] |
| observables only | ['C', 'U2', 'season'] |
Without restrictions, the function returns \(\{N, \text{season}, t\}\), which contains the unobservable \(N\). Restricted to observed variables, it returns \(\{C, U_2, \text{season}\}\).
The function returns a minimal separator, not all of them. The graph is small enough (\(2^{10}\) subsets) to enumerate every minimal separator:
# `find_minimal_d_separator` returns *a* minimal separator, not *the* minimal
# one. On a graph this size the whole lattice is enumerable (2**10 subsets),
# so count them all with `nx.is_minimal_d_separator` instead of over-reading
# one answer.
def minimal_separators(candidates):
"""Every minimal d-separator of (treatment, outcome) drawn from ``candidates``."""
return sorted(
[
sorted(Z)
for size in range(len(candidates) + 1)
for Z in itertools.combinations(sorted(candidates), size)
if nx.is_minimal_d_separator(G_do, {treatment}, {outcome}, set(Z))
],
key=lambda Z: (len(Z), Z),
)
pd.DataFrame(
[
{"restriction": "unrestricted", "minimal separator": str(Z)}
for Z in minimal_separators(set(G) - {treatment, outcome})
]
+ [
{"restriction": "observables only", "minimal separator": str(Z)}
for Z in minimal_separators(observable | {"U2", "season"})
]
).set_index("restriction")
| minimal separator | |
|---|---|
| restriction | |
| unrestricted | ['C', 'U2', 'season'] |
| unrestricted | ['D', 'U2', 'season'] |
| unrestricted | ['MS', 'U2', 'season'] |
| unrestricted | ['N', 'season', 't'] |
| observables only | ['C', 'U2', 'season'] |
Without restrictions, the enumeration finds four minimal separators of the same size. The first answer above, \(\{N, \text{season}, t\}\), contains the unobservable \(N\). However, \(\{C, U_2, \text{season}\}\) is just as small, so the latent variable in the first answer only reflects the order of the search.
The other two, \(\{D, U_2, \text{season}\}\) and \(\{M^{*}, U_2, \text{season}\}\), fail the no-descendant condition, because \(D\) and \(M^{*}\) are descendants of the treatment (see the caption of the verdict table). Hence, every minimal separator is unobservable, contains a descendant, or is \(\{C, U_2, \text{season}\}\).
Restricted to observed variables that are not descendants, \(\{C, U_2, \text{season}\}\) is the only minimal separator, and it needs \(U_2\). This is the same \(U_2\) that opens two paths as a collider in the path count. Both statements hold, because a separation test and a path count answer different questions.
In short, the graph recommends the conditioning set of Naive B+. It does not need \(t\), because conditioning on the observed category index blocks the trend strand and the noise strand at once.
Note
Library versions: networkx is not pinned. A new version could change the order of the open-path listing and which separator find_minimal_d_separator returns first. The counts and the verdicts do not depend on the version.
Naive B++ is the cheaper specification, because its control column comes from the calendar and needs no extra data series. However, it closes the trend fork and leaves the noise fork open. Therefore, using it as a total-effect estimator assumes that the noise correlation of the drivers is negligible.
We measured that correlation above: \(|r| \le 0.18\) per geo against a floor of \(0.13\), and \(|r| \le 0.13\) pooled against \(0.11\). A trend control rests on this assumption, so we state it as an assumption and not as a verdict.
We fit both Naive B+ and Naive B++ below, so we can check the ranking of the graph against the truth.
Note that the verdict loop checked all six verdicts with \(U_2\) as the treatment too, while we run the direct-effect and separator checks for \(U_1\) only.
Note
The funnel model has no adjustment set: the funnel model specifies the full structural system instead. This is correct only if the demand equation is correctly specified (including the role of \(C\)) and nothing unobserved confounds \(M^{*} \to Y\).
The backdoor checks cannot verify these assumptions. Here they hold by construction. On real data, posterior predictive checks of the proxies can reveal some violations, but not an unobserved confounder of \(M^{*} \to Y\).
Testable Implication#
The graph also gives us a check on the data. The exclusion restriction “the budget enters lower-funnel spend but not branded search” implies that \(B\) and \(S\) are independent marginally. Their only connecting path runs through the collider at \(M^{*}\), which stays blocked as long as we do not condition on anything downstream of it:
pd.DataFrame(
{
"B ⊥ S given Z": {
"Z = {}": nx.is_d_separator(G, {"B"}, {"S"}, set()),
"Z = {Y}": nx.is_d_separator(G, {"B"}, {"S"}, {"Y"}),
"Z = {M}": nx.is_d_separator(G, {"B"}, {"S"}, {"M"}),
}
}
)
| B ⊥ S given Z | |
|---|---|
| Z = {} | True |
| Z = {Y} | False |
| Z = {M} | False |
The graph confirms that \(B\) and \(S\) are independent marginally, but not given \(Y\) or given M, because conditioning on a descendant of \(M^{*}\) opens the collider. Next, we check the implication on the raw data, per geo:
pd.DataFrame(
{
"corr(lf_budget, search_volume)": {
geo: np.corrcoef(lf_budget[:, gi], search_obs[:, gi])[0, 1]
for gi, geo in enumerate(geos)
}
}
).style.format("{:+.3f}")
| corr(lf_budget, search_volume) | |
|---|---|
| north | -0.070 |
| south | +0.046 |
| west | -0.101 |
All three correlations are small: \(-0.070\) in north, \(+0.046\) in south and \(-0.101\) in west. Even the largest is only a little larger than the sampling standard error of about \(0.09\) on \(130\) weeks.
The restriction holds by construction in the DGP, so these departures from zero come from the draw, not from a failed exclusion restriction. On real data, the check matters, because a systematic failure of the restriction would show up here.
Here are two important remarks. The check works per geo only, because pooled regions make market size a common cause of both series. It also works on raw series only, because conditioning on a descendant of \(M^{*}\), including the target, opens the collider.
Note
Automated adjustment sets: for production use, CausalGraphModel selects adjustment sets from a DAG string automatically. Its recommendation is coarser than the checks above. We ran the checks with networkx to keep the reasoning visible.
Model Fitting#
We now fit the six candidate models on the same synthetic data. The five naive models differ only in which observed series they condition on, so the differences among them come from structure.
The funnel model observes more than any of them: two extra likelihood series (lower_spend and search_volume) and two extra inputs (category_demand and lf_budget). Where it does better, the gain comes from structure and information together. Keep this in mind when you read the comparisons.
First, we build the long data frame of each specification and set the sampler: \(4\) chains with \(1{,}000\) tuning steps and \(1{,}000\) draws each, target_accept=0.95 and the nutpie sampler.
frame = (
data[
[
"tv_spend",
"social_spend",
"lower_spend",
"search_volume",
"category_demand",
"y",
]
]
.to_dataframe()
.reset_index()
)
X_upper = frame[["date", "geo", *channels]]
X_with_mediator = frame[["date", "geo", *channels, "lower_spend"]]
X_with_search = frame[["date", "geo", *channels, "search_volume"]]
X_with_control = frame[["date", "geo", *channels, "category_demand"]]
# The growth trend as a control regressor: observable from the calendar alone.
frame["t"] = frame["date"].map(dict(zip(date_range, t, strict=True)))
X_with_trend = frame[["date", "geo", *channels, "t"]]
y = frame["y"]
sample_kwargs = dict(chains=4, tune=1_000, draws=1_000, target_accept=0.95)
sampler_config = {"nuts_sampler": "nutpie"}
Naive A#
Naive A (mediator as channel) lists upper- and lower-funnel spend side by side as channels, a common choice when both appear in the same table. In this DAG, lower-funnel spend drives the target, so the coefficient that the model estimates for it has a causal meaning.
However, as the d-separation analysis showed, conditioning on the mediator blocks the indirect path. The upper-funnel channels get credit only for their direct path, and the lower-funnel channel keeps the credit for the demand that they created.
In other words, Naive A answers a different question, the direct effect. We score it against the direct effect as well as against the total effect.
naive_a = make_mmm([*channels, "lower_spend"], sampler_config=sampler_config)
naive_a.build_model(X_with_mediator, y)
naive_a.add_original_scale_contribution_variable(var=["channel_contribution", "y"])
naive_a.fit(X_with_mediator, y, random_seed=rng, **sample_kwargs)
naive_a.sample_posterior_predictive(X_with_mediator, random_seed=rng);
Naive B#
In Naive B (mediator omitted), the two upper-funnel channels are the only media. This conditioning set targets the total effect, because a total effect does not condition on the mediator.
Naive B has two problems, and they affect different channels. The first is the open backdoor through the shared growth trend, which the residual-correlation table showed is large for tv_spend and at noise level for social_spend.
The graph says nothing about the second problem. Without the mediator, each response curve must also absorb the demand that its channel creates downstream, through a saturating pool that the model does not represent.
naive_b = make_mmm(channels, sampler_config=sampler_config)
naive_b.build_model(X_upper, y)
naive_b.add_original_scale_contribution_variable(var=["channel_contribution", "y"])
naive_b.fit(X_upper, y, random_seed=rng, **sample_kwargs)
naive_b.sample_posterior_predictive(X_upper, random_seed=rng);
Naive B+#
Naive B+ (demand as control) is the adjustment that the d-separation analysis recommended. We block the remaining backdoors with the observed category index, which lies downstream of both forks of the shared driver.
MMM takes the index as an ordinary control column, a linear term with a per-geo coefficient and no adstock or saturation:
naive_bplus = make_mmm(
channels, control_columns=["category_demand"], sampler_config=sampler_config
)
naive_bplus.build_model(X_with_control, y)
naive_bplus.add_original_scale_contribution_variable(var=["channel_contribution", "y"])
naive_bplus.fit(X_with_control, y, random_seed=rng, **sample_kwargs)
naive_bplus.sample_posterior_predictive(X_with_control, random_seed=rng);
Here are two important remarks before we see the results. First, on the graph, Naive B+ is a valid adjustment set for the total effect, at the cost of one extra column. If you only need the total ROAS of each upper-funnel channel, a specification like this is the baseline that the funnel model has to beat.
Second, its control is a linear stand-in for a nonlinear path: \(C\) enters through the demand equation and then through a saturating conversion. Therefore, some bias can remain even with a valid graph verdict, because d-separation says which variables suffice and says nothing about functional form.
Naive B++#
Naive B++ (trend as control) is the cheapest specification here, because its regressor is a function of the calendar and needs no extra data series. On the full graph, it is valid only if the shared noise of the drivers is negligible.
Trend controls usually rest on this kind of assumption, and we measured its size above. The functional-form remark about Naive B+ also applies here: the trend reaches the target through softplus-rectified plans and a saturating conversion, and a linear term stands in for all of it.
naive_bpp = make_mmm(channels, control_columns=["t"], sampler_config=sampler_config)
naive_bpp.build_model(X_with_trend, y)
naive_bpp.add_original_scale_contribution_variable(var=["channel_contribution", "y"])
naive_bpp.fit(X_with_trend, y, random_seed=rng, **sample_kwargs)
naive_bpp.sample_posterior_predictive(X_with_trend, random_seed=rng);
Naive C#
Naive C (indicator as channel) adds branded search volume to the channel list. This variable is tempting, because it correlates strongly with sales, it costs nothing and it arrives every week.
However, it causes nothing. It is a descendant of latent demand, so conditioning on it partly blocks \(U \to D \to M^{*} \to Y\), as conditioning on a noisy record of the mediator does. This is post-treatment conditioning on a variable that has no causal role at all.
naive_c = make_mmm([*channels, "search_volume"], sampler_config=sampler_config)
naive_c.build_model(X_with_search, y)
naive_c.add_original_scale_contribution_variable(var=["channel_contribution", "y"])
naive_c.fit(X_with_search, y, random_seed=rng, **sample_kwargs)
naive_c.sample_posterior_predictive(X_with_search, random_seed=rng);
Funnel Model#
The funnel model has the two upper-funnel channels (their direct effects) plus the BudgetDualIndicatorEffect, which reads the mediator data from the Dataset and fits the target and both proxies jointly. As described in the Data Plumbing section, we build the model from the Dataset and fit it with the long frame.
ds_fit = make_dataset(
media_raw, lower_spend_obs, search_obs, category_demand, lf_budget
)
y_fit = xr.DataArray(y_obs, dims=("date", "geo"), coords=coords)
funnel = make_mmm(channels, sampler_config=sampler_config)
funnel.add_mu_effect(make_budget_dual_indicator_effect())
funnel.build_model(X=ds_fit, y=y_fit)
funnel.add_original_scale_contribution_variable(
var=["channel_contribution", "budget_dual_effect_contribution", "y"]
)
funnel.fit(X=X_upper, y=y, random_seed=rng, **sample_kwargs)
funnel.sample_posterior_predictive(X_upper, random_seed=rng);
Model Diagnostics#
Before looking into the results, let’s check some diagnostics. We count the divergent transitions and compute the largest \(\hat{R}\) over the scalar parameters of each model:
models = {
"Naive A": naive_a,
"Naive B": naive_b,
"Naive B+": naive_bplus,
"Naive B++": naive_bpp,
"Naive C": naive_c,
"Funnel": funnel,
}
# Scalar parameters only: summarizing the date x geo deterministics would
# dominate this cell's runtime.
scalar_params = [
"intercept_contribution",
"adstock_alpha",
"saturation_lam",
"saturation_beta",
"gamma_fourier",
"gamma_control",
"y_sigma",
"budget_dual_baseline",
"budget_dual_parent_coef",
"budget_dual_budget_coef",
"budget_dual_indicator_coef",
"budget_dual_spend_measurement_sigma",
"budget_dual_hub_indicator_sigma",
"adstock_hub_input_alpha",
"sat_hub_input_lam",
"sat_hub_input_beta",
"adstock_contribution_alpha",
"sat_contribution_lam",
"sat_contribution_beta",
]
sampler_diagnostics = {}
for label, model in models.items():
present = [v for v in scalar_params if v in model.idata.posterior]
sampler_diagnostics[label] = {
"divergences": int(model.idata["sample_stats"]["diverging"].sum()),
"max r-hat": float(az.summary(model.idata, var_names=present)["r_hat"].max()),
}
pd.DataFrame(sampler_diagnostics).T.style.format({"max r-hat": "{:.3f}"})
| divergences | max r-hat | |
|---|---|---|
| Naive A | 0.000000 | 1.000 |
| Naive B | 0.000000 | 1.000 |
| Naive B+ | 0.000000 | 1.000 |
| Naive B++ | 0.000000 | 1.010 |
| Naive C | 0.000000 | 1.000 |
| Funnel | 0.000000 | 1.010 |
All six fits have \(0\) divergent transitions, and the largest \(\hat{R}\) is at most \(1.01\) for every model, which is good!
Posterior Predictive Checks#
As our model and posterior samples are looking good, we can now look into the posterior predictive checks. The funnel model has three observed quantities, and it should recover all of them. First, the target:
funnel.plot_suite = "new"
pc = funnel.plot.diagnostics.posterior_predictive(
hdi_prob=0.94, return_as_pc=True, figsize=(13, 8)
)
fig = pc.viz["figure"].item()
for i, geo in enumerate(geos):
ax = pc.viz["plot"].sel(geo=geo).item()
ax.set_xlabel("date" if i == len(geos) - 1 else "")
ax.set_ylabel("")
ax.set_title(geo)
fig.supylabel("target")
fig.suptitle(
"Funnel model: target posterior predictive", fontsize=14, fontweight="bold"
)
fig.autofmt_xdate()
fig.set_layout_engine("tight", rect=(0.01, 0, 1, 0.97))
The posterior predictive mean follows the observed target closely in all three geos, and the \(94\%\) HDI band covers almost every week. The posterior predictive target looks quite reasonable.
The other two observed quantities are the likelihoods of the effect. They have no built-in plot, so we draw them by hand. Most figures from here on show a posterior HDI, so we first define a small helper for the interval bounds:
def hdi_bounds(da, prob: float = 0.94):
"""Return (lower, upper) HDI DataArrays for a posterior DataArray."""
h = az.hdi(da, prob=prob)
if isinstance(h, xr.Dataset):
h = h[next(iter(h.data_vars))]
return h.sel(ci_bound="lower"), h.sel(ci_bound="upper")
Next, we plot the posterior predictive of the two demand proxies, lower_spend and search_volume, against the observed series:
The posterior mean of both proxies follows the observed series in every geo, including the promotional spikes of lower_spend. The observed series stay inside the \(94\%\) HDI bands in almost every week.
Parameter Recovery#
Because we generated the data, we can check the funnel model against the truth. The MMM fits on an internally scaled target, so we must rescale some parameters before we compare them with the truth.
We multiply every parameter that acts on the scale of the target mean by target_scale: the intercept, the two saturation amplitudes and the Fourier coefficients. With dims=() scaling, this factor is per geo, and xarray applies it by label.
Parameters that act on the inputs (adstock decays and saturation curvatures) or on the mediator scale (everything in the demand equation) are scale-free, so we compare them directly. The first figure shows the demand and conversion equations, with the truth as a dashed line:
The posteriors of the two loadings of the new structure, budget_dual_budget_coef and budget_dual_indicator_coef, are narrow and close to the truth, which is good! These loadings decide how much of lower-funnel spend is demand-driven and how strongly branded search tracks demand.
Not every parameter of the new block is this sharp. For example, the posterior mean of adstock_contribution_alpha in west is \(0.159\) against a true \(0.250\), with a wide interval that still covers the truth.
Next, we check that the data identify \(\lambda\), because an unidentified parameter can also have a posterior that contains the truth. We look at how much the data narrow \(\lambda\) relative to its prior, and at whether \(\lambda\) trades off against the demand baseline, which shares its level:
prior_lambda_sd = 0.5 * np.sqrt(1 - 2 / np.pi) # sd of HalfNormal(sigma=0.5)
lambda_draws = funnel.idata.posterior["budget_dual_budget_coef"]
identification = pd.DataFrame(
{
"posterior mean": lambda_draws.mean(dim=["chain", "draw"]).to_series(),
"true": pd.Series(true["budget_dual_budget_coef"], index=geos),
"posterior sd / prior sd": (
lambda_draws.std(dim=["chain", "draw"]) / prior_lambda_sd
).to_series(),
"corr(lambda, baseline)": pd.Series(
{
geo: float(
xr.corr(
lambda_draws.sel(geo=geo),
funnel.idata.posterior["budget_dual_baseline"].sel(geo=geo),
)
)
for geo in geos
}
),
}
)
identification.style.format("{:.3f}")
| posterior mean | true | posterior sd / prior sd | corr(lambda, baseline) | |
|---|---|---|---|---|
| north | 0.515 | 0.500 | 0.042 | -0.072 |
| south | 0.447 | 0.450 | 0.045 | 0.042 |
| west | 0.371 | 0.400 | 0.095 | -0.067 |
The posterior sd of \(\lambda\) is \(4\%\) to \(10\%\) of its prior sd (\(0.042\) to \(0.095\)). Its correlation with the demand baseline is close to zero (\(-0.072\) to \(+0.042\)).
The budget exclusion restriction explains this. The budget enters lower-funnel spend but not branded search, so the search equation fixes the baseline while the promotional spikes of the budget fix \(\lambda\). With a flat, always-on budget, the two parameters would not be separately identified, and more sampling would not help.
Next, let’s look at the media transformations:
recovery_plot(
[
"adstock_hub_input_alpha",
"sat_hub_input_lam",
"sat_hub_input_beta",
"adstock_alpha",
"saturation_lam",
"saturation_beta",
],
"Parameter recovery: media transformations (dashed = truth)",
(14, 14),
);
Most true values of the media parameters fall inside the bulk of their posteriors. The demand-path amplitudes sat_hub_input_beta of social_spend have long right tails, and the direct amplitude saturation_beta of social_spend in south sits clearly above its true value.
The two figures have \(37\) panels, so we also show the same information as numbers: the posterior mean, the \(94\%\) HDI, the true value and whether the interval covers it.
| mean | hdi 3% | hdi 97% | true | covered | ||
|---|---|---|---|---|---|---|
| parameter | cell | |||||
| budget_dual_baseline | north | 0.217 | 0.104 | 0.330 | 0.200 | True |
| south | 0.209 | 0.139 | 0.270 | 0.160 | True | |
| west | 0.147 | 0.080 | 0.212 | 0.120 | True | |
| budget_dual_parent_coef | north | 0.297 | 0.225 | 0.379 | 0.350 | True |
| south | 0.306 | 0.202 | 0.413 | 0.300 | True | |
| west | 0.237 | 0.059 | 0.421 | 0.250 | True | |
| budget_dual_budget_coef | north | 0.515 | 0.493 | 0.540 | 0.500 | True |
| south | 0.447 | 0.420 | 0.472 | 0.450 | True | |
| west | 0.371 | 0.316 | 0.424 | 0.400 | True | |
| budget_dual_indicator_coef | none | 0.798 | 0.791 | 0.805 | 0.800 | True |
| adstock_contribution_alpha | north | 0.326 | 0.187 | 0.459 | 0.350 | True |
| south | 0.265 | 0.126 | 0.393 | 0.300 | True | |
| west | 0.159 | 0.016 | 0.304 | 0.250 | True | |
| sat_contribution_lam | none | 1.123 | 0.604 | 1.653 | 1.500 | True |
| sat_contribution_beta | north | 0.974 | 0.715 | 1.252 | 1.050 | True |
| south | 1.108 | 0.748 | 1.456 | 0.900 | True | |
| west | 0.880 | 0.460 | 1.372 | 0.720 | True | |
| adstock_hub_input_alpha | tv_spend | 0.610 | 0.569 | 0.656 | 0.600 | True |
| social_spend | 0.337 | 0.241 | 0.429 | 0.350 | True | |
| sat_hub_input_lam | tv_spend | 2.693 | 2.301 | 3.133 | 2.500 | True |
| social_spend | 2.347 | 1.026 | 3.605 | 3.000 | True | |
| sat_hub_input_beta | north, tv_spend | 0.880 | 0.771 | 0.990 | 0.850 | True |
| north, social_spend | 0.632 | 0.435 | 0.831 | 0.550 | True | |
| south, tv_spend | 0.746 | 0.644 | 0.847 | 0.700 | True | |
| south, social_spend | 0.394 | 0.252 | 0.540 | 0.450 | True | |
| west, tv_spend | 0.578 | 0.478 | 0.670 | 0.550 | True | |
| west, social_spend | 0.344 | 0.216 | 0.482 | 0.350 | True | |
| adstock_alpha | tv_spend | 0.575 | 0.457 | 0.694 | 0.550 | True |
| social_spend | 0.344 | 0.197 | 0.489 | 0.300 | True | |
| saturation_lam | tv_spend | 3.163 | 2.290 | 4.033 | 3.000 | True |
| social_spend | 4.110 | 2.570 | 5.516 | 4.000 | True | |
| saturation_beta | north, tv_spend | 0.600 | 0.465 | 0.734 | 0.550 | True |
| north, social_spend | 0.353 | 0.183 | 0.534 | 0.350 | True | |
| south, tv_spend | 0.403 | 0.245 | 0.547 | 0.450 | True | |
| south, social_spend | 0.470 | 0.261 | 0.696 | 0.300 | True | |
| west, tv_spend | 0.366 | 0.238 | 0.503 | 0.350 | True | |
| west, social_spend | 0.285 | 0.121 | 0.458 | 0.250 | True |
All \(37\) intervals contain their true values. Most direct-path parameters are less precise than their funnel-side counterparts (compare the interval widths), because both paths start at the same upper-funnel spend. The posterior can be much more certain about the total response of a channel than about how this response splits between the two paths.
The path amplitudes can trade off, but they do not trade off in every cell. In south, social_spend, the demand-path amplitude (sat_hub_input_beta) is about \(12\%\) below the truth, while the direct-path amplitude (saturation_beta) is about \(57\%\) above it.
The amplitude errors also have opposite signs in south, tv_spend and, barely, in west, social_spend. In the other three cells, the two errors have the same sign. For this reason, we check the total effects directly below, instead of inferring them from one loading.
The smallest market has the weakest signal by design, but it is not the worst recovered market in this run. The observation noise y_sigma is \(0.04\) in every geo, while the media inputs shrink with market size (\(1.0\), \(0.7\) and \(0.45\) for north, south and west).
Yet the mean relative parameter errors in the caption are \(7.7\%\), \(17.1\%\) and \(13.2\%\), so south has the largest error, not west. This is consistent with sampling variation on top of the signal-to-noise design.
Coverage does not mean that every posterior mean is exact. The largest total-effect error below is in south, social_spend, where the direct loading is high. For a budget decision, the relevant quantity is the total effect, so we check it next.
Media Deep Dive#
Now that we are happy with the model fit, we can look into the media and individual channel contributions.
Channel Contributions#
First, we reconstruct the estimated decomposition of each channel. The direct path is the base channel contribution. The indirect path is the mediated contribution minus the mediated contribution with the spend of that channel set to zero.
This is the same counterfactual that we ran on the generative model, applied to the posterior. The incrementality module automates it in the ROAS section below, and we compute it by hand once to show what the module computes.
The manual version also gives the induced lower-funnel spend \(\Delta M^{*}\). The module does not report it, because its increment is in target units while the induced spend is a spend series.
def counterfactual_indirect(channel: str) -> dict[str, xr.DataArray]:
"""Posterior mediated contribution and induced LF spend attributable to ``channel``."""
X_zero = X_upper.copy()
X_zero[channel] = 0.0
cf = funnel.sample_posterior_predictive(
X_zero,
extend_idata=False,
combined=False,
var_names=[
"budget_dual_effect_contribution_original_scale",
"budget_dual_spend",
],
random_seed=rng,
)
factual = funnel.idata.posterior["budget_dual_effect_contribution_original_scale"]
counterfactual = cf["budget_dual_effect_contribution_original_scale"]
# Guard the plumbing: when X is a DataFrame the conversion keeps only the channel
# columns, so the effect's own data variables (including lf_budget) correctly
# retain their training values. If that ever stopped holding, the mediated
# contribution would not move and the indirect effect would silently read zero.
if float(counterfactual.mean()) >= float(factual.mean()):
raise RuntimeError(
f"zeroing {channel} did not reduce the mediated contribution: "
"the effect's data variables were probably overwritten"
)
return {
"indirect": factual - counterfactual,
# M* with the channel on minus M* with it off: the lower-funnel spend the
# channel *induces*. Same units as observed lower_spend (the likelihood
# ties them), which is what lets it enter a ROAS denominator later.
"induced_spend": (
funnel.idata.posterior["budget_dual_spend"] - cf["budget_dual_spend"]
),
}
counterfactuals = {channel: counterfactual_indirect(channel) for channel in channels}
direct_post = funnel.idata.posterior["channel_contribution_original_scale"]
indirect_post = xr.concat(
[counterfactuals[channel]["indirect"] for channel in channels],
dim=pd.Index(channels, name="channel"),
)
induced_spend_post = xr.concat(
[counterfactuals[channel]["induced_spend"] for channel in channels],
dim=pd.Index(channels, name="channel"),
)
total_post = direct_post + indirect_post
We plot the posterior mean of the direct and indirect contributions as a stack, with the \(94\%\) HDI of the total and the true direct and total contributions. We print the minimum to check that both components are non-negative, which stackplot needs:
min component: 0.03384489044858784
The top of the stack is the posterior mean of the total, and the two thin lines bound its \(94\%\) HDI. These are two different summaries, so the top of the stack does not have to sit midway between the lines.
For tv_spend, the stack follows the true total (dashed) closely in all three geos. For social_spend, the HDI contains the true total in every geo, but in south the stack sits well above it, because the direct part is too large.
The indirect part (orange) is a large share of the stack in every panel. In south, social_spend, the direct part (blue) alone sits above the true direct contribution (dotted), which explains the high total.
Contribution Bias#
Next, we compare the six models on the quantity that a budget decision depends on: the total contribution of each upper-funnel channel. Here, bias means the error of the posterior mean against the truth on this dataset.
For the naive models, the channel contribution is already the counterfactual against zero spend. For the funnel model, we add the counterfactual indirect effect to the direct contribution.
We group the rows by estimand, following the d-separation analysis. The models above the separator answer the total-effect question, or claim to. Naive A sits below the separator with the direct path of the funnel model, because it answers the direct-effect question.
Every bias annotation names its benchmark. The figure is a forest plot of posterior means and intervals against the true values. The ROAS comparison below uses the same panel, so we define a helper function first:
Next, we collect the mean weekly contribution of every model and channel, and the true direct and total contributions:
def total_contribution(model: MMM, channel: str) -> xr.DataArray:
"""Mean-over-time total contribution of ``channel`` for a naive model."""
return (
model.idata.posterior["channel_contribution_original_scale"]
.sel(channel=channel)
.mean(dim="date")
)
estimates = {
"Naive B (mediator omitted)": lambda ch: total_contribution(naive_b, ch),
"Naive B+ (demand as control)": lambda ch: total_contribution(naive_bplus, ch),
"Naive B++ (trend as control)": lambda ch: total_contribution(naive_bpp, ch),
"Naive C (indicator as channel)": lambda ch: total_contribution(naive_c, ch),
"Funnel (total)": lambda ch: total_post.sel(channel=ch).mean(dim="date"),
"Naive A (mediator as channel)": lambda ch: total_contribution(naive_a, ch),
"Funnel (direct only)": lambda ch: direct_post.sel(channel=ch).mean(dim="date"),
}
CONTRIBUTION_BENCHMARKS = {
"Naive A (mediator as channel)": "direct",
"Funnel (direct only)": "direct",
}
truth_mean = data["total_true"].mean("date")
direct_truth_mean = data["direct_true"].mean("date")
We plot the estimates for each channel and geo, with the total-effect rows above the separator and the direct-effect rows below it:
Every model that ignores the mediator is biased, and the two channels show different failure mechanisms.
Naive B overstates tv_spend in every geo, by \(+16\%\) to \(+31\%\). The correlation of tv_spend with category demand stays at \(+0.63\) to \(+0.67\) after the seasonal adjustment, so the model credits the shared growth trend to a channel that also trends.
Both controls reduce the tv_spend bias, and the trend control reduces it more. Naive B++ moves tv_spend to \(-6\%\) to \(-2\%\) of the truth, while Naive B+, the valid adjustment set of the graph, leaves \(+8\%\) to \(+20\%\). In west, Naive B+ (\(+18\%\)) is even slightly further from the truth than Naive B (\(+16\%\)).
The graph tells us what to condition on, but not how strongly a linear control competes with flexible media curves for shared variation. The coefficient table below shows this competition.
On social_spend, the missing shared mediated response matters. Naive B overstates social_spend by \(+89\%\) in south and \(+35\%\) in west, while north is almost exact (\(-2\%\)). Naive B+ and Naive B++ make the social_spend errors larger in south and west.
The backdoor of social_spend is open on the graph, but its residual correlation with demand is at noise level (\(|r| \le 0.12\)). The shared noise is negative for social_spend in this draw (pooled \(-0.126\)), so it cannot explain a positive bias either.
A model with a separate response curve per channel cannot represent one shared, saturating mediated pool, and an adjustment column cannot add this structure. The tables in target units below show where the omitted demand goes.
Naive C partly blocks the mediator and leaves the backdoor of tv_spend open. Its tv_spend errors range from \(+10\%\) to \(+21\%\), and its social_spend errors from \(-14\%\) in north to \(+69\%\) in south. With mixed signs and sizes, no single correction factor could fix it.
Naive A answers the direct-effect question approximately. As a total-effect estimate, it understates tv_spend, because conditioning on the mediator blocks the indirect path. Against the true direct effect, its errors still range from \(+7\%\) to \(+86\%\), because conditioning on the noisy M is not the same as conditioning on the structural \(M^{*}\).
The funnel model has the right structure, but it still has posterior error in this run, in particular in south, social_spend (\(+27\%\)). Its direct and mediated paths share the upper-funnel spend, so the split between them can stay loose even when the total is well recovered.
We use the error of the funnel model as a reference below, not as an amount that we can subtract from every naive model.
Control Coefficients#
Why does the valid control of Naive B+ correct only part of the tv_spend bias? To answer this question, we first derive what its coefficient should be.
In the true DGP, a sustained increase of one unit in category demand raises latent demand by \(\gamma_g\), and therefore raises \(M^{*}\) by \(\gamma_g\). The geometric adstock is normalized, so a sustained shift passes through with a total weight of one.
The target then rises by \(\gamma_g\) times the slope of the true conversion curve at its operating point. LogisticSaturation is \(\beta \tanh(\lambda x / 2)\), so this slope is \(\beta\,(\lambda/2)/\cosh^{2}(\lambda x / 2)\), which we average over the adstocked \(M^{*}\) series. We call the result the implied slope.
The control of Naive B++ is the calendar trend itself, so its implied slope goes one step further back. A sustained increase of one unit in \(t\) raises the category index by its trend loading times the softplus derivative, times the Jacobian of the rescaling of the index.
The increase then reaches the target through the same chain from demand to \(M^{*}\) to conversion, with the same averaging. This makes the two implied slopes comparable.
Note
Two caveats on the implied slopes: first, the normalized adstock passes a sustained shift with a total weight of one only from week l_max on. The zero-padded start passes slightly less, which changes the value by less than half a percent here.
Second, we average the pointwise slope, so we get the mean of the slope and not the slope at the mean, and there is no Jensen gap. This is still not the estimand of the regression, which weights the slope by where the variance is. The two have the same order of magnitude but are not identical.
adstock_weights = np.power(
true["adstock_contribution_alpha"][None, :], np.arange(l_max)[:, None]
)
adstock_weights /= adstock_weights.sum(axis=0)
mstar = data["lf_spend_true"].transpose("date", "geo").to_numpy()
adstocked_mstar = np.stack(
[
np.convolve(mstar[:, g], adstock_weights[:, g])[: len(mstar)]
for g in range(len(geos))
],
axis=1,
)
lam_lf = float(true["sat_contribution_lam"])
pointwise_slope = (
true["sat_contribution_beta"]
* (lam_lf / 2)
/ np.cosh(lam_lf * adstocked_mstar / 2) ** 2
)
conversion_slope = pointwise_slope.mean(axis=0)
# The softplus derivative is recovered as 1 - exp(-softplus_output), on the
# index with the geo rescaling undone, the same division the noise cell
# applies before inverting.
softplus_out = x_data[..., 2] / geo_size[None, :]
softplus_slope = (1.0 - np.exp(-softplus_out)).mean(axis=0) * (
geo_size / x_data[..., 2].max()
)
implied_slopes = {
"Naive B+": true["budget_dual_parent_coef"] * conversion_slope,
"Naive B++": true["budget_dual_parent_coef"]
* trend[drivers.index("category_demand")]
* softplus_slope
* conversion_slope,
}
control_of = {
"Naive B+": (naive_bplus, "category_demand"),
"Naive B++": (naive_bpp, "t"),
}
Next, we compare the posterior of each control coefficient, on the original target scale, with its implied slope. The table also gives the distance from the implied slope and from zero in posterior sds, and the prior sd relative to the implied slope:
# The default control prior is Normal(0, 2) on the *scaled* target, so its sd in
# the original units the implied slope lives in is 2 * target_scale.
control_blocks = {}
for control_label, (control_model, control_name) in control_of.items():
control_scale = control_model.idata.constant_data["target_scale"]
gamma_scaled = (
control_model.idata.posterior["gamma_control"].sel(control=control_name)
* control_scale
)
post_mean = gamma_scaled.mean(dim=["chain", "draw"]).to_series()
post_sd = gamma_scaled.std(dim=["chain", "draw"]).to_series()
implied = pd.Series(implied_slopes[control_label], index=pd.Index(geos, name="geo"))
prior_sd = (2.0 * control_scale).to_series()
control_blocks[control_label] = pd.DataFrame(
{
"posterior mean (target scale)": post_mean,
"implied true slope": implied,
"posterior sd (target scale)": post_sd,
"sd from implied": (post_mean - implied).abs() / post_sd,
"sd from zero": post_mean.abs() / post_sd,
"prior sd (target scale)": prior_sd,
"prior sd / implied": prior_sd / implied,
}
)
summary = pd.concat(control_blocks, names=["model"])
summary.style.format(
{col: "{:.3f}" for col in summary.columns} | {"prior sd / implied": "{:.0f}x"}
)
| posterior mean (target scale) | implied true slope | posterior sd (target scale) | sd from implied | sd from zero | prior sd (target scale) | prior sd / implied | ||
|---|---|---|---|---|---|---|---|---|
| model | geo | |||||||
| Naive B+ | north | 0.136 | 0.118 | 0.057 | 0.317 | 2.379 | 4.367 | 37x |
| south | 0.210 | 0.119 | 0.088 | 1.035 | 2.379 | 3.530 | 30x | |
| west | 0.027 | 0.098 | 0.112 | 0.633 | 0.242 | 2.561 | 26x | |
| Naive B++ | north | 0.144 | 0.031 | 0.037 | 3.021 | 3.859 | 4.367 | 139x |
| south | 0.093 | 0.022 | 0.039 | 1.816 | 2.373 | 3.530 | 161x | |
| west | 0.058 | 0.012 | 0.026 | 1.755 | 2.211 | 2.561 | 216x |
The implied slope of Naive B+ is small, about \(0.1\) target units per unit of category demand (\(0.098\) to \(0.119\)), because a unit of demand converts through the saturating lower-funnel curve. The posterior means scatter around this value: above it in north and south (\(0.136\) and \(0.210\)) and below it in west (\(0.027\)).
Each posterior mean is within about one posterior sd of the implied slope, so the distance from zero is more informative. In north and south, the coefficient is more than two posterior sds from zero (\(2.38\) in both geos), so the control has a measurable effect there.
In west, the coefficient is only \(0.24\) sds from zero, so the data cannot tell this control apart from an omitted column. With this sample size, we cannot tell a control that does its full job from one that does half of it, or, in west, from one that does nothing.
Note
Shrinkage is not informative here: we do not read the shrinkage from prior to posterior as evidence, as we did for \(\lambda\). The default \(\text{Normal}(0, 2)\) control prior is on the scaled target, so in original units its sd is \(2 \cdot \texttt{target\_scale}\).
This is \(26\) to \(37\) times larger than the implied slope of Naive B+ (last column). A posterior that is a few percent of that width is still wide compared with the slope that it tries to measure.
The rows of Naive B++ show the other side. Its implied slope is even smaller, because it goes through the trend loading of the category index, the softplus derivative and the rescaling.
However, the posterior means of Naive B++ are far above its implied slope in all three geos. They are \(0.144\), \(0.093\) and \(0.058\) against \(0.031\), \(0.022\) and \(0.012\), about \(4\) to \(5\) times larger. Their posterior sds are between about \(2/3\) and \(1/4\) of those of Naive B+.
Therefore, the linear term in \(t\) absorbs more of the trend than the demand path alone implies. This is consistent with the term also absorbing trend co-movement that the response curves cannot mimic. In west, the category index competes with the same curves for shared variance and recovers only a small part of its slope (\(0.027\) against \(0.098\)).
The control and the flexible tv_spend response curve compete for the same trend variance, and the posterior splits it between them.
In north and south, the tv_spend bias of Naive B+ lies between those of Naive B and Naive B++. In north, the three values are \(+31\%\), \(+20\%\) and \(-6\%\), and in south they are \(+24\%\), \(+8\%\) and \(-6\%\).
In west, where the control coefficient is \(0.24\) posterior sds from zero, Naive B+ hardly moves, and it moves in the wrong direction (\(+16\%\) to \(+18\%\)).
The trend control of Naive B++ takes more of this variance, because a response curve has to bend to mimic a trend, while \(t\) enters linearly. This is why the conditionally valid control closes the tv_spend gap more fully than the valid one.
In other words, a valid control can be estimated about right and still have too small an effect. Validity says nothing about how strongly the likelihood uses a linear term when nonlinear terms can partly mimic it.
Neither control fixes the social_spend side, as the tables below show. There the backdoor was open on the graph but at noise level in the measurement, so the controls had little confounding to remove.
Return on Ad Spend (ROAS)#
To have a more quantitative metric of efficiency, we turn to ROAS. All estimates use the same definition: the contribution summed over the observed weeks, divided by the spend over the same weeks. This way, any difference between the estimates comes from the models, not from the definition.
The standard rows divide by upper-funnel spend only. They report the return per upper-funnel dollar and treat the lower-funnel spend that upper-funnel activity induces as free.
The coefficient on \(D\) in the \(M^{*}\) equation is one, so every unit of demand that a channel creates is a unit of extra lower-funnel spend that someone has to approve. About a third of the effect of each channel travels this path, so the induced spend is large.
Only the funnel model can price this induced spend. Therefore, next to Funnel (total), we also report a row that adds the induced spend to the denominator: \(\Delta y \,/\, (\text{channel spend} + \Delta M^{*})\).
This row assumes that all spend series use the same currency, so that we can add upper-funnel and induced lower-funnel spend in one denominator. On real data, this usually holds.
In this synthetic DGP, the scale of each series sets the exchange rate between a unit of tv_spend and a unit of lower-funnel spend, so this rate is arbitrary. Read the induced-spend row as a demonstration of the accounting, not as a meaningful price.
The next cell computes the ROAS of every model with the incrementality module, the induced-spend row of the funnel model and the true ROAS values:
spend = data[channels].to_array("channel").sum("date")
def incrementality_roas(model: MMM) -> xr.DataArray:
"""All-time ROAS via the incrementality module's spend counterfactual.
For a model whose effects opted in, this is the total effect.
"""
return model.incrementality.contribution_over_spend(frequency="all_time").sel(
channel=channels
)
def roas_from(contribution: xr.DataArray) -> xr.DataArray:
"""All-time ROAS: summed contribution over summed spend, per channel and geo."""
return contribution.sel(channel=channels).sum(dim="date") / spend
# Denominator including each channel's induced lower-funnel spend. The denominator
# is itself a posterior quantity: the model estimates how much M* the channel moved.
denom_incl = spend + induced_spend_post.sum(dim="date")
# `BudgetDualIndicatorEffect` opted in through `incrementality_spec`, so for the funnel model
# the built-in follows the spend counterfactual down both paths and returns the
# total (direct + mediated). The direct-only row is assembled from the posterior
# channel contributions instead, and only the induced-spend denominator still
# needs the manual counterfactual: the module's increment is in target units,
# while the induced spend is a spend series it does not report.
roas = {
"Naive B (mediator omitted)": incrementality_roas(naive_b),
"Naive B+ (demand as control)": incrementality_roas(naive_bplus),
"Naive B++ (trend as control)": incrementality_roas(naive_bpp),
"Naive C (indicator as channel)": incrementality_roas(naive_c),
"Funnel (total)": incrementality_roas(funnel),
"Funnel (total / all incremental spend)": (
total_post.sel(channel=channels).sum(dim="date") / denom_incl
),
"Naive A (mediator as channel)": incrementality_roas(naive_a),
"Funnel (direct only)": roas_from(
funnel.idata.posterior["channel_contribution_original_scale"]
),
}
ROAS_BENCHMARKS = {
"Funnel (total / all incremental spend)": "total incl. LF",
"Naive A (mediator as channel)": "direct",
"Funnel (direct only)": "direct",
}
true_total_roas = data["total_true"].sum("date") / spend
true_direct_roas = data["direct_true"].sum("date") / spend
true_denom_incl = spend + data["induced_spend_true"].sum("date")
true_incl_roas = data["total_true"].sum("date") / true_denom_incl
Note
One call for every model: every row comes from the same built-in call, mmm.incrementality.contribution_over_spend(frequency="all_time"), which runs an on/off spend counterfactual per channel. For the naive models, this is their channel ROAS.
For the funnel model, BudgetDualIndicatorEffect opted in through incrementality_spec, so the counterfactual follows both paths and Funnel (total) is the module output itself. Without the opt-in, the call would raise an error instead of reporting the direct path as if it were the total.
Two rows use other inputs. Funnel (direct only) divides the posterior channel contributions by spend, and the induced-spend row needs the posterior counterfactual of the previous section.
The table below compares the module output with the direct-plus-indirect ROAS that we assembled by hand. With frequency="all_time", the evaluation window is the whole series, so the carryover extension of the module has nothing to extend into.
Hence, the windowed attribution and the whole-series zeroing compute the same difference, draw by draw, and the two columns must agree up to rounding. The carryover extension only matters when a window boundary falls inside the data, as the yearly breakdown below shows.
manual_total = roas_from(total_post)
pd.DataFrame(
{
"incrementality, all_time (direct + mediated)": roas["Funnel (total)"]
.mean(dim=["chain", "draw"])
.transpose("channel", "geo")
.to_series(),
"manual counterfactual (direct + indirect)": manual_total.mean(
dim=["chain", "draw"]
)
.transpose("channel", "geo")
.to_series(),
}
).style.format("{:.4f}")
| incrementality, all_time (direct + mediated) | manual counterfactual (direct + indirect) | ||
|---|---|---|---|
| channel | geo | ||
| tv_spend | north | 0.4220 | 0.4220 |
| south | 0.5053 | 0.5053 | |
| west | 0.6707 | 0.6707 | |
| social_spend | north | 0.5277 | 0.5277 |
| south | 0.8352 | 0.8352 | |
| west | 0.8472 | 0.8472 |
The two columns agree to four decimals in all six cells, as expected.
To see include_carryover at work, we ask for yearly windows. The 2021 and 2022 windows end inside the data, and the carryover option credits each year with the response that its spend causes after the year ends.
This reach is l_max weeks of direct adstock plus the extra weeks of the second adstock of the mediator, which is the window that the declared additional_carryover_lags sets. Because the effect opted in, these yearly values are totals (direct plus mediated), like the all-time rows above:
# The same integer seed for both calls, so they subsample the same posterior
# draws and the difference column is not polluted by Monte Carlo noise. A shared
# `rng` would not do here: the subsampling helper consumes the generator it is
# handed, so the second call would land on different draws.
yearly_with = funnel.incrementality.contribution_over_spend(
frequency="yearly",
num_samples=500,
random_state=seed,
)
yearly_without = funnel.incrementality.contribution_over_spend(
frequency="yearly",
include_carryover=False,
num_samples=500,
random_state=seed,
)
yearly = pd.DataFrame(
{
"with carryover": yearly_with.mean(dim=["chain", "draw"])
.sel(channel=channels)
.to_series(),
"without": yearly_without.mean(dim=["chain", "draw"])
.sel(channel=channels)
.to_series(),
}
)
yearly["difference"] = yearly["with carryover"] - yearly["without"]
yearly.style.format("{:.4f}")
| with carryover | without | difference | |||
|---|---|---|---|---|---|
| date | channel | geo | |||
| 2021-12-31 00:00:00 | tv_spend | north | 0.4817 | 0.4688 | 0.0129 |
| south | 0.5640 | 0.5493 | 0.0147 | ||
| west | 0.7456 | 0.7273 | 0.0183 | ||
| social_spend | north | 0.5620 | 0.5571 | 0.0049 | |
| south | 0.8794 | 0.8722 | 0.0072 | ||
| west | 0.8924 | 0.8872 | 0.0052 | ||
| 2022-12-31 00:00:00 | tv_spend | north | 0.4297 | 0.4209 | 0.0088 |
| south | 0.5190 | 0.5074 | 0.0116 | ||
| west | 0.6777 | 0.6652 | 0.0125 | ||
| social_spend | north | 0.5251 | 0.5204 | 0.0047 | |
| south | 0.8464 | 0.8384 | 0.0080 | ||
| west | 0.8728 | 0.8646 | 0.0082 | ||
| 2023-12-31 00:00:00 | tv_spend | north | 0.3444 | 0.3444 | 0.0000 |
| south | 0.4100 | 0.4100 | 0.0000 | ||
| west | 0.5491 | 0.5491 | 0.0000 | ||
| social_spend | north | 0.4579 | 0.4579 | 0.0000 | |
| south | 0.7304 | 0.7304 | 0.0000 | ||
| west | 0.7501 | 0.7501 | 0.0000 |
The 2021 and 2022 values differ between the two columns by \(0.0047\) to \(0.0183\), as expected.
Note two details. First, the window of the final, partial year ends at the calendar boundary, after the last observation. Its two columns are identical for the same reason as the all-time columns: after the data there is nothing to extend into.
Second, with carryover on, the yearly windows overlap, because the carryover weeks of each year lie inside the window of the next year. Therefore, the yearly values do not add up to the all-time value.
National ROAS#
First, we aggregate to a spend-weighted national ROAS per model. Each row keeps its own denominator, so we weight the induced-spend row by spend plus induced spend:
# National = total contribution over total spend; each row keeps its own
# denominator, so the induced-spend row is weighted by spend *plus* induced spend.
national = {
label: (da * spend).sum("channel").sum("geo") / spend.sum()
for label, da in roas.items()
if label != "Funnel (total / all incremental spend)"
}
national["Funnel (total / all incremental spend)"] = (
roas["Funnel (total / all incremental spend)"] * denom_incl
).sum(["channel", "geo"]) / denom_incl.sum(["channel", "geo"])
national = {label: national[label] for label in roas} # restore row order
national_true_total = float((true_total_roas * spend).sum() / spend.sum())
national_true_direct = float((true_direct_roas * spend).sum() / spend.sum())
national_true_incl = float(
(true_incl_roas * true_denom_incl).sum() / true_denom_incl.sum()
)
fig, ax = plt.subplots(figsize=(11, 6.5))
forest_panel(
ax,
national,
reference={
"total": (national_true_total, TRUE_STYLE["total"]),
"direct": (national_true_direct, TRUE_STYLE["direct"]),
"total incl. LF": (
national_true_incl,
{"color": "black", "linestyle": "-.", "linewidth": 1.5, "alpha": 0.6},
),
},
benchmarks=ROAS_BENCHMARKS,
)
# Estimand separator: total-effect rows above, direct-effect rows below.
ax.axhline(len(national) - 2.5, color="gray", linestyle=":", lw=1)
ax.set_xlabel("ROAS")
ax.legend(loc="lower right")
ax.set_title(
"Spend-weighted national ROAS by model (94% HDI)",
fontsize=15,
fontweight="bold",
)
fig.set_layout_engine("tight")
The rows above the separator estimate a total effect, and we compare the first five with the true total ROAS (dashed line). Naive B and Naive B+ overstate it by \(31\%\) and \(30\%\), Naive B++ by \(20\%\) and Naive C by \(18\%\), while Funnel (total) is \(6\%\) above the truth.
The sixth row divides the same total contribution by all the incremental money: channel spend plus induced lower-funnel spend. We compare it with its own benchmark (dash-dotted line), and it is \(7\%\) above that benchmark.
The gap between the two funnel rows is the price of the induced spend. The return per upper-funnel dollar and the return per incremental dollar answer different questions, and only a model that sees the funnel can answer the second one.
Below the separator are the two direct-effect rows, which we compare with the true direct ROAS (dotted line). Following the d-separation analysis, Naive A belongs here.
Both direct rows are above the truth, Naive A by \(29\%\) and Funnel (direct only) by \(14\%\). The direct path inherits the amplitude trade-off of the recovery table.
A correct decomposition targets the right estimand, but one dataset still leaves estimation error, as the \(+14\%\) of Funnel (direct only) shows.
Warning
Critical Assessment: as budget inputs, the naive models miss the truth by different amounts and in different directions. Nothing in the models themselves tells you which error you are looking at.
Geo-Level ROAS#
The national number averages over the geo differences that a geo panel can show, so we also break it out. contribution_over_spend returns (chain, draw, channel, geo), so we only need a different view of the same draws. We leave the induced-spend row out of these panels to keep them readable:
The per-cell picture repeats the contribution figure. The Funnel (total) interval contains the true total ROAS in every panel.
The largest misses of the naive models are in south, social_spend, where Naive B, Naive B+ and Naive B++ are \(89\%\) to \(104\%\) above the truth. The direct rows also miss most in this cell.
Next, we show the same comparison as a table of levels:
roas_summary = pd.DataFrame(
{
label: da.mean(dim=["chain", "draw"]).transpose("channel", "geo").to_series()
for label, da in roas.items()
}
)
roas_summary["true direct"] = true_direct_roas.transpose("channel", "geo").to_series()
roas_summary["true total"] = true_total_roas.transpose("channel", "geo").to_series()
roas_summary["true total incl. LF"] = true_incl_roas.transpose(
"channel", "geo"
).to_series()
roas_summary.style.format("{:.3f}")
| Naive B (mediator omitted) | Naive B+ (demand as control) | Naive B++ (trend as control) | Naive C (indicator as channel) | Funnel (total) | Funnel (total / all incremental spend) | Naive A (mediator as channel) | Funnel (direct only) | true direct | true total | true total incl. LF | ||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| channel | geo | |||||||||||
| tv_spend | north | 0.519 | 0.476 | 0.375 | 0.480 | 0.422 | 0.306 | 0.329 | 0.290 | 0.259 | 0.397 | 0.294 |
| south | 0.612 | 0.534 | 0.466 | 0.546 | 0.505 | 0.343 | 0.334 | 0.285 | 0.309 | 0.496 | 0.349 | |
| west | 0.726 | 0.741 | 0.615 | 0.688 | 0.671 | 0.419 | 0.482 | 0.423 | 0.394 | 0.627 | 0.406 | |
| social_spend | north | 0.552 | 0.593 | 0.572 | 0.486 | 0.528 | 0.363 | 0.390 | 0.371 | 0.365 | 0.565 | 0.378 |
| south | 1.242 | 1.317 | 1.345 | 1.113 | 0.835 | 0.610 | 0.784 | 0.667 | 0.421 | 0.659 | 0.428 | |
| west | 1.129 | 1.175 | 1.237 | 1.025 | 0.847 | 0.558 | 0.733 | 0.637 | 0.555 | 0.836 | 0.501 |
In south, social_spend, every model overstates the true total ROAS of \(0.659\), with estimates from \(0.784\) (Naive A) to \(1.345\) (Naive B++).
Next, we compute the bias of every model against the true total ROAS. For five of the six models, this is their own estimand.
We include Naive A on the same scale on purpose, because analysts use a mediator-as-channel model as a total-effect model in practice. Its column shows the cost of reading it that way, and we score it against its own estimand separately below.
First, the mean absolute error per channel:
total_estimates = [
"Naive A (mediator as channel)",
"Naive B (mediator omitted)",
"Naive B+ (demand as control)",
"Naive B++ (trend as control)",
"Naive C (indicator as channel)",
"Funnel (total)",
]
bias = roas_summary[total_estimates].div(roas_summary["true total"], axis=0) - 1
bias.loc[("mean |error|", ""), :] = bias.abs().mean()
# The summary row averages across both channels, while the comparison this
# notebook is organized around is per channel: one failure mechanism each.
by_channel_error = bias.loc[channels].abs().groupby(level="channel", sort=False).mean()
by_channel_error.T.style.format("{:.1%}").set_caption("Mean |error| by channel")
| channel | tv_spend | social_spend |
|---|---|---|
| Naive A (mediator as channel) | 24.3% | 20.8% |
| Naive B (mediator omitted) | 23.3% | 42.0% |
| Naive B+ (demand as control) | 15.3% | 48.4% |
| Naive B++ (trend as control) | 4.4% | 51.1% |
| Naive C (indicator as channel) | 13.5% | 35.2% |
| Funnel (total) | 5.0% | 11.6% |
The funnel model has the smallest social_spend error (\(11.6\%\)). On tv_spend, Naive B++ (\(4.4\%\)) is slightly more accurate than the funnel model (\(5.0\%\)).
Every naive model has a social_spend error between \(20.8\%\) (Naive A) and \(51.1\%\). Next, the full table per geo and channel, with the mean absolute error over all six cells in the last row:
bias.style.background_gradient(
cmap="RdBu_r", vmin=-0.5, vmax=0.5, subset=pd.IndexSlice[channels, :]
).format("{:+.1%}")
| Naive A (mediator as channel) | Naive B (mediator omitted) | Naive B+ (demand as control) | Naive B++ (trend as control) | Naive C (indicator as channel) | Funnel (total) | ||
|---|---|---|---|---|---|---|---|
| channel | geo | ||||||
| tv_spend | north | -17.1% | +30.6% | +19.9% | -5.5% | +20.8% | +6.3% |
| south | -32.6% | +23.5% | +7.7% | -5.9% | +10.1% | +1.9% | |
| west | -23.2% | +15.8% | +18.2% | -1.9% | +9.7% | +6.9% | |
| social_spend | north | -31.1% | -2.4% | +4.9% | +1.2% | -14.0% | -6.7% |
| south | +19.0% | +88.6% | +100.0% | +104.3% | +69.0% | +26.8% | |
| west | -12.4% | +34.9% | +40.4% | +47.9% | +22.5% | +1.3% | |
| mean |error| | +22.5% | +32.6% | +31.9% | +27.8% | +24.3% | +8.3% |
The funnel model has the smallest mean absolute error over all six cells (\(8.3\%\)), followed by Naive A against the total (\(22.5\%\)) and Naive C (\(24.3\%\)). Its largest miss is south, social_spend (\(+26.8\%\)).
Read the summary row of Naive A with care. A percentage bias divides by its benchmark, and the true direct effect is about two thirds of the true total (the indirect shares are \(33.7\%\) to \(37.7\%\)). As a result, the same absolute error looks smaller against the total than against the direct effect.
Bias in Target Units#
A percentage bias is a ratio, and the two channels do not share a denominator. Before we compare tv_spend and social_spend, we measure the errors in target units. This is the estimated minus the true contribution, which is what a budget discussion would get wrong.
The first table sums the signed errors over geos:
spend_series = spend.transpose("channel", "geo").to_series()
true_contribution = (
data["total_true"].sum("date").transpose("channel", "geo").to_series()
)
abs_bias = (
roas_summary[total_estimates]
.sub(roas_summary["true total"], axis=0)
.mul(spend_series, axis=0)
)
by_channel = abs_bias.groupby(level="channel", sort=False).sum().loc[channels]
true_by_channel = (
true_contribution.groupby(level="channel", sort=False).sum().loc[channels]
)
by_channel.T.style.format("{:+.1f}").set_caption(
"Signed (estimated - true) contribution, summed over geos; true totals: "
+ ", ".join(f"{ch} {true_by_channel[ch]:.1f}" for ch in channels)
)
| channel | tv_spend | social_spend |
|---|---|---|
| Naive A (mediator as channel) | -38.8 | -13.8 |
| Naive B (mediator omitted) | +39.0 | +55.0 |
| Naive B+ (demand as control) | +24.8 | +66.9 |
| Naive B++ (trend as control) | -7.5 | +69.8 |
| Naive C (indicator as channel) | +22.8 | +34.0 |
| Funnel (total) | +8.1 | +9.5 |
In truth, the two channels contribute \(161.8\) and \(145.5\) target units. Naive B overstates tv_spend by \(39.0\) and social_spend by \(55.0\) target units.
Note
Signed and absolute sums: when we sum signed errors over geos, the errors of a model can cancel. A model that is wrong in both directions can then net to a small value. The uncancelled sum below adds the absolute error of each cell, so a lucky sign cannot hide an error.
In both tables, we sum at full precision and round at the end, so a total can differ from the sum of the printed cells by \(0.1\).
uncancelled = abs_bias.abs().groupby(level="channel", sort=False).sum().loc[channels]
uncancelled_table = uncancelled.T
uncancelled_table["both channels"] = abs_bias[total_estimates].abs().sum()
uncancelled_table.style.format("{:.1f}").set_caption(
"Misstated contribution, summed as |bias| per cell over geos"
)
| channel | tv_spend | social_spend | both channels |
|---|---|---|---|
| Naive A (mediator as channel) | 38.8 | 32.0 | 70.8 |
| Naive B (mediator omitted) | 39.0 | 57.8 | 96.8 |
| Naive B+ (demand as control) | 24.8 | 66.9 | 91.7 |
| Naive B++ (trend as control) | 7.5 | 69.8 | 77.3 |
| Naive C (indicator as channel) | 22.8 | 50.1 | 72.9 |
| Funnel (total) | 8.1 | 17.2 | 25.3 |
The funnel model misstates \(25.3\) target units over both channels, against \(70.8\) to \(96.8\) for the naive models. This ranking of the models is the same as the ranking by mean absolute percentage error, so the good rank of Naive A does not come from the denominator.
Next, the signed error of every cell:
abs_bias.style.background_gradient(cmap="RdBu_r", vmin=-30, vmax=30).format("{:+.1f}")
| Naive A (mediator as channel) | Naive B (mediator omitted) | Naive B+ (demand as control) | Naive B++ (trend as control) | Naive C (indicator as channel) | Funnel (total) | ||
|---|---|---|---|---|---|---|---|
| channel | geo | ||||||
| tv_spend | north | -10.6 | +19.0 | +12.4 | -3.4 | +12.9 | +3.9 |
| south | -17.7 | +12.8 | +4.2 | -3.2 | +5.5 | +1.0 | |
| west | -10.5 | +7.1 | +8.2 | -0.8 | +4.4 | +3.1 | |
| social_spend | north | -17.9 | -1.4 | +2.8 | +0.7 | -8.1 | -3.8 |
| south | +9.1 | +42.5 | +48.0 | +50.1 | +33.1 | +12.9 | |
| west | -4.9 | +13.9 | +16.0 | +19.0 | +8.9 | +0.5 |
The contrast between the channels does not come from the denominators, because the true contributions of the two channels are similar. Naive B misstates social_spend by \(57.8\) target units (uncancelled), against \(39.0\) on tv_spend.
Naive B++ corrects tv_spend (\(7.5\) units) but leaves \(69.8\) units of social_spend error. This is consistent with the posterior putting the extra demand on social_spend when the shared mediated response is missing.
Signed totals can hide errors that cancel across geos. The social_spend errors of the funnel model are \(-3.8\), \(+12.9\) and \(+0.5\) target units, so they sum to \(+9.5\) signed but \(17.2\) uncancelled. Naive C has \(+34.0\) signed against \(50.1\) uncancelled.
The three tv_spend errors of Naive B++ are all negative, so its signed and uncancelled tv_spend totals have the same size, \(7.5\). Use the uncancelled table to judge how wrong a single allocation cell can be.
Adjustment also moves attribution between channels. From Naive B to Naive B+, about \(14\) target units of tv_spend overstatement go away and about \(12\) units are added to social_spend. From Naive B to Naive B++, about \(46.5\) units go away from tv_spend and about \(15\) are added to social_spend.
Both controls address a large tv_spend backdoor, but neither can represent the shared mediated pool, so the naive model still has to attribute this demand to some channel. Closing the backdoor moves part of the misattributed demand from tv_spend to social_spend, but it does not remove the error.
Carryover Absorption#
If a channel curve absorbs the mediated response, its carryover should look longer than the direct carryover alone. The table shows the posterior of the pooled direct-path parameters of social_spend in every model, next to the true values:
tilt_params = ["adstock_alpha", "saturation_lam"]
social_index = channels.index("social_spend")
tilt_cols = pd.MultiIndex.from_product([tilt_params, ["mean", "sd"]])
tilt_rows = {
"true": [
float(true["adstock_alpha"][social_index]),
np.nan,
float(true["saturation_lam"][social_index]),
np.nan,
]
}
for label, model in models.items():
posterior = model.idata.posterior
row = []
for param in tilt_params:
draws = posterior[param].sel(channel="social_spend")
row += [float(draws.mean()), float(draws.std())]
tilt_rows[label] = row
pd.DataFrame.from_dict(tilt_rows, orient="index", columns=tilt_cols).style.format(
"{:.3f}", na_rep=""
).set_caption("Pooled direct-path parameters on social_spend")
| adstock_alpha | saturation_lam | |||
|---|---|---|---|---|
| mean | sd | mean | sd | |
| true | 0.300 | 4.000 | ||
| Naive A | 0.353 | 0.074 | 4.137 | 0.731 |
| Naive B | 0.404 | 0.054 | 4.444 | 0.606 |
| Naive B+ | 0.423 | 0.052 | 4.523 | 0.603 |
| Naive B++ | 0.458 | 0.057 | 4.644 | 0.612 |
| Naive C | 0.388 | 0.060 | 4.329 | 0.651 |
| Funnel | 0.344 | 0.080 | 4.110 | 0.785 |
Every model overestimates the carryover adstock_alpha of social_spend (true \(0.300\)). The funnel model is the closest (\(0.344\)), and Naive B++ is the furthest (\(0.458\)).
To judge these values, we compute a mean-lag ruler from the true kernels. It asks which single geometric alpha has the same average delay as two references: the direct path alone, and a curve that must also absorb the mediated response. The carryover of the mediated response is the upper-funnel adstock followed by the lower-funnel adstock.
Note
What the ruler measures: the ruler compares kernel shapes, so it predicts an order of size, not a fitted value. Only in the funnel model is adstock_alpha a pure direct-path parameter. In a naive model, it measures whatever the curve has to absorb.
kernel_lags = np.arange(l_max)
alpha_true = float(true["adstock_alpha"][social_index])
def normalized_kernel(alpha):
weights = alpha**kernel_lags
return weights / weights.sum()
def kernel_mean_lag(kernel):
return float((kernel_lags[: len(kernel)] * kernel).sum())
def alpha_with_mean_lag(target_lag):
low, high = 0.001, 0.999
for _ in range(60):
mid = (low + high) / 2
if kernel_mean_lag(normalized_kernel(mid)) < target_lag:
low = mid
else:
high = mid
return (low + high) / 2
direct_lag = kernel_mean_lag(normalized_kernel(alpha_true))
indirect_share = (
(data["indirect_true"].mean("date") / data["total_true"].mean("date"))
.sel(channel="social_spend")
.to_numpy()
)
uf_kernel = normalized_kernel(float(true["adstock_hub_input_alpha"][social_index]))
mediated_lags = []
for lf_alpha in true["adstock_contribution_alpha"]:
mediated_kernel = np.convolve(uf_kernel, normalized_kernel(lf_alpha))[:l_max]
mediated_lags.append(kernel_mean_lag(mediated_kernel / mediated_kernel.sum()))
blend_lags = [
(1 - s) * direct_lag + s * m
for s, m in zip(indirect_share, mediated_lags, strict=True)
]
pd.DataFrame(
[
[alpha_true] * len(geos),
[alpha_with_mean_lag(lag) for lag in blend_lags],
[alpha_with_mean_lag(lag) for lag in mediated_lags],
],
index=[
"direct path alone",
"blend at the true indirect share",
"mediated path alone",
],
columns=pd.Index(geos, name="geo"),
).style.format("{:.3f}").set_caption(
"Mean-lag ruler for the carryover column (alpha equivalents; "
"the direct row is geo-independent)"
)
| geo | north | south | west |
|---|---|---|---|
| direct path alone | 0.300 | 0.300 | 0.300 |
| blend at the true indirect share | 0.397 | 0.384 | 0.367 |
| mediated path alone | 0.527 | 0.498 | 0.470 |
The ruler goes from \(0.300\) for the direct path to \(0.470\) to \(0.527\) for the mediated path alone. A blend at the true indirect share of each geo gives \(0.397\), \(0.384\) and \(0.367\). A single channel curve that must represent both paths should therefore move toward a larger alpha.
Important
The ruler cannot measure absorption: the \(94\%\) HDI of the funnel model for adstock_alpha of social_spend goes from \(0.197\) to \(0.489\). This interval includes the direct reference, all three blends and part of the mediated range.
We cannot place any model on one row of the ruler. We also cannot convert the distance between a naive model and the funnel model into an amount of absorbed mediated response.
The ordering of the posterior means still fits the mechanism: Funnel \(0.344\), Naive A \(0.353\), Naive C \(0.388\), then Naive B \(0.404\), Naive B+ \(0.423\) and Naive B++ \(0.458\). The curvature saturation_lam has the same ordering (\(4.110\), \(4.137\), \(4.329\), \(4.444\), \(4.523\) and \(4.644\)).
However, the differences within the Naive B family cannot all come from mediated absorption, because these three models omit the mediator in the same way and differ only in their controls. Their increasing alpha goes together with the shift of trend variation that we saw above.
The posterior sds in the table also warn against reading small differences between model means as firm rankings.
The \(17.2\) units of uncancelled social_spend error of the funnel model are posterior-mean errors inside its own intervals, not a measured common offset. Hence, we cannot split the \(57.8\) units of Naive B or the \(69.8\) units of Naive B++ into a part from the draw and a part from the specification.
Their much larger errors, together with the missing shared mediated response, show that the missing mediation matters. They do not measure how much bias it causes alone.
Interval Coverage#
A point estimate that is off is one problem, and an estimate that is confidently off is a worse one. The tables below check, for every geo and channel, whether the \(94\%\) HDI of each model contains the true total ROAS, and how wide the intervals are.
Note that these are six correlated cells from a single simulated dataset. Treat the counts as a sanity check, not as a calibration study, which would need many replicate datasets.
coverage_rows = []
for c in channels:
for g in geos:
truth = float(true_total_roas.sel(channel=c, geo=g))
row = {}
for label in total_estimates:
lo, hi = hdi_bounds(roas[label].sel(channel=c, geo=g))
row[(label, "covers")] = bool(float(lo) <= truth <= float(hi))
row[(label, "width")] = float(hi) - float(lo)
coverage_rows.append(row)
coverage = pd.DataFrame(
coverage_rows,
index=pd.MultiIndex.from_product([channels, geos], names=["channel", "geo"]),
)
coverage.columns = pd.MultiIndex.from_tuples(coverage.columns)
pd.DataFrame(
{
"covers (of 6)": {
label: int(coverage[(label, "covers")].sum()) for label in total_estimates
},
"mean 94% HDI width": {
label: coverage[(label, "width")].mean() for label in total_estimates
},
}
).style.format({"mean 94% HDI width": "{:.3f}"}).set_caption(
"Coverage of the true total ROAS, per model"
)
| covers (of 6) | mean 94% HDI width | |
|---|---|---|
| Naive A (mediator as channel) | 5 | 0.471 |
| Naive B (mediator omitted) | 4 | 0.544 |
| Naive B+ (demand as control) | 5 | 0.596 |
| Naive B++ (trend as control) | 5 | 0.675 |
| Naive C (indicator as channel) | 6 | 0.567 |
| Funnel (total) | 6 | 0.422 |
The funnel model contains the true total ROAS in all six cells with the narrowest intervals on average (mean width \(0.422\)). Naive C also covers all six cells, but with wider intervals (\(0.567\)), and Naive B covers only four of six.
Next, we score Naive A against the true direct ROAS, which the d-separation analysis says it estimates:
# Naive A, scored against the quantity the d-separation analysis says it
# actually estimates, the true *direct* ROAS:
label = "Naive A (mediator as channel)"
a_bias = (roas_summary[label] / roas_summary["true direct"] - 1).loc[channels]
a_covers = sum(
bool(float(lo) <= float(true_direct_roas.sel(channel=c, geo=g)) <= float(hi))
for c in channels
for g in geos
for lo, hi in [hdi_bounds(roas[label].sel(channel=c, geo=g))]
)
pd.DataFrame(
{"mean |error|": [a_bias.abs().mean()], "covers (of 6)": [a_covers]},
index=["Naive A vs its own estimand (true direct)"],
).style.format({"mean |error|": "{:.1%}"})
| mean |error| | covers (of 6) | |
|---|---|---|
| Naive A vs its own estimand (true direct) | 30.5% | 6 |
Naive A contains the true direct ROAS in all six cells, but its mean absolute error against it is \(30.5\%\). Finally, the coverage and the width of every cell:
coverage.style.format({(m, "width"): "{:.3f}" for m in total_estimates}).map(
lambda v: (
f"background-color: {'#c8e6c9' if v else '#ffcdd2'}"
if isinstance(v, bool)
else ""
)
)
| Naive A (mediator as channel) | Naive B (mediator omitted) | Naive B+ (demand as control) | Naive B++ (trend as control) | Naive C (indicator as channel) | Funnel (total) | ||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| covers | width | covers | width | covers | width | covers | width | covers | width | covers | width | ||
| channel | geo | ||||||||||||
| tv_spend | north | True | 0.163 | False | 0.167 | True | 0.190 | True | 0.253 | True | 0.207 | True | 0.138 |
| south | False | 0.255 | True | 0.263 | True | 0.332 | True | 0.505 | True | 0.319 | True | 0.202 | |
| west | True | 0.319 | True | 0.276 | True | 0.342 | True | 0.445 | True | 0.292 | True | 0.247 | |
| social_spend | north | True | 0.420 | True | 0.558 | True | 0.604 | True | 0.623 | True | 0.614 | True | 0.430 |
| south | True | 0.802 | False | 0.993 | False | 1.064 | False | 1.074 | True | 0.990 | True | 0.718 | |
| west | True | 0.868 | True | 1.006 | True | 1.045 | True | 1.149 | True | 0.983 | True | 0.795 | |
Warning
Critical Assessment: read the bias and coverage tables together. A narrow interval can still miss the relevant truth, and a wide interval can cover the truth with an inaccurate posterior mean.
Naive A answers a different question. Against the true total ROAS, it misses tv_spend in south with an interval width of \(0.255\) and an estimate \(32.6\%\) too low. Its narrowest tv_spend interval, in north (\(0.163\)), covers the truth even though its mean is \(17.1\%\) too low.
Against its own direct-effect estimand, Naive A covers all six cells but has a mean error of \(30.5\%\). Coverage does not turn a rough direct-effect estimate into a good total-effect estimate.
Naive B has an open backdoor and a missing mediator. tv_spend is too high by \(+15.8\%\) to \(+30.6\%\), and social_spend reaches \(+88.6\%\) in south. It covers only four of six true totals.
The mean interval width of Naive B is \(0.544\), about \(29\%\) wider than the \(0.422\) of the funnel model. The wider intervals still miss two of the six true totals.
On average, both controls improve tv_spend and make social_spend worse. Naive B++ brings tv_spend to \(-5.9\%\) to \(-1.9\%\), while Naive B+ leaves \(+7.7\%\) to \(+19.9\%\). Both cover all three tv_spend cells and miss social_spend in south, so each covers five of six.
The mean absolute social_spend errors of Naive B+ and Naive B++ rise to \(48.4\%\) and \(51.1\%\), from \(42.0\%\) for Naive B. On tv_spend, Naive B++ is slightly more accurate than the funnel model in this run (\(4.4\%\) against \(5.0\%\) mean error, and \(7.5\) against \(8.1\) target units).
The mean error of Naive C hides large differences between cells. Its \(24.3\%\) combines cells from \(-14.0\%\) to \(+69.0\%\), with no stable sign or size that a correction factor could use. It offers neither the decomposition of the funnel model nor its induced-spend accounting.
The funnel model covers all six true totals with the narrowest mean interval (\(0.422\), against \(0.567\) for Naive C and \(0.544\) for Naive B). Its mean absolute error is \(8.3\%\), against \(22.5\%\) for Naive A scored against the total and \(24.3\%\) for Naive C.
The funnel model also has extra spend and search likelihoods that give it information that the naive models do not have. Hence, we cannot credit the narrower intervals to structure alone.
The funnel model is not the best in every cell. Its largest miss is south, social_spend (\(+26.8\%\)), the cell whose direct loading is high in the recovery table. Among the tv_spend cells, west has the widest interval (\(0.247\)) and, with north, the largest error (\(+6.9\%\) and \(+6.3\%\)).
For social_spend, west has the widest interval (\(0.795\)) and a small error (\(+1.3\%\)). north has the narrowest interval (\(0.430\)) and an error of \(-6.7\%\). In other words, the width follows the uncertainty of the channel more than the accuracy of each cell.
Business Value: the funnel model has the smallest overall ROAS error and the narrowest mean interval width, and it is the only model that can price the induced lower-funnel spend.
Predictive Fit#
Can predictive fit tell the models apart? We compute the in-sample \(R^2\) of the posterior predictive mean of every model:
r2 = {}
for label, model in models.items():
pp = (
model.idata.posterior_predictive["y"]
* model.idata.constant_data["target_scale"]
)
y_hat = pp.mean(dim=["chain", "draw"])
y_true, y_hat = xr.align(data["y"], y_hat, join="exact")
ss_res = float(((y_true - y_hat) ** 2).sum())
ss_tot = float(((y_true - y_true.mean()) ** 2).sum())
r2[label] = 1 - ss_res / ss_tot
pd.Series(r2, name="in-sample $R^2$").to_frame().style.format("{:.4f}")
| in-sample $R^2$ | |
|---|---|
| Naive A | 0.9857 |
| Naive B | 0.9791 |
| Naive B+ | 0.9798 |
| Naive B++ | 0.9802 |
| Naive C | 0.9790 |
| Funnel | 0.9871 |
The six models span less than one point of in-sample \(R^2\) (\(0.9790\) to \(0.9871\)), while their ROAS estimates differ by tens of percent.
The four models that do not observe the mediator are almost identical (\(0.9790\) to \(0.9802\), a spread of \(0.0012\)). Conditioning on \(C\) or \(t\) adds almost nothing to the fit, even though the trend control corrects the tv_spend estimate.
The two models that observe the mediator are a small but clear step higher (\(0.9857\) and \(0.9871\)). Note that this grouping by fit differs from the carryover ordering above, where Naive C sits closer to the funnel model. Predicting the target and absorbing the mediated response are different things.
Fit does carry some information: it shows that the mediator predicts the target. But it cannot separate Naive A from the funnel model, which are \(0.0014\) apart in \(R^2\) and tens of percent apart in the total value of tv_spend.
In-sample fit measures the predictive content of the mediator, not the causal structure around it. Two remarks apply: this is in-sample \(R^2\), and the funnel model has two extra likelihood terms, so this is not a like-for-like model-selection comparison.
The Nürnberger Versicherung case study, which faced the same funnel problem in production, makes the same point:
“Trust is not created by \(R^2\) values. It is created when business reality matches model expectations.” (Philip Herp, Nürnberger Versicherung)
Built-In Decomposition#
Almost every causal number in this notebook comes from post-processing the posterior. The exception is the incrementality module, which ran the per-channel spend counterfactual for us once the effect opted in. Let’s see what the decomposition of the library reports for the same fitted model.
decomposition.contributions_over_time and decomposition.waterfall know a fixed set of components: channels, baseline, controls and seasonality. A custom MuEffect is none of these.
The counterfactual decomposition of MMM does include a registered effect as its own part, but the Model Limitations section explains why this does not solve the attribution either. The waterfall for north:
funnel.plot.decomposition.waterfall(dims={"geo": ["north"]}, figsize=(11, 5));
The waterfall shows the baseline (\(43.7\%\)), tv_spend (\(30.3\%\)), social_spend (\(25.4\%\)) and seasonality (\(0.6\%\)), with a total of about \(1.2\). The channel bars are the direct paths only.
The waterfall does not move the mediated contribution into the baseline either. It leaves it out, so the bars do not add up to the target, as the next table shows:
scale = funnel.idata.constant_data["target_scale"]
post = funnel.idata.posterior
components = {
"intercept": post["intercept_contribution"] * scale,
"channels (direct)": post["channel_contribution_original_scale"].sum("channel"),
"seasonality": post["yearly_seasonality_contribution"] * scale,
}
# Bars are rounded for display and summed after rounding, so the printed column
# adds up exactly.
bars = {
name: round(float(da.sel(geo="north").mean()), 3) for name, da in components.items()
}
mediated = float(
post["budget_dual_effect_contribution_original_scale"].sel(geo="north").mean()
)
two_dp_rows = ["mediated, missing from the bars", "observed target mean"]
waterfall_summary = pd.DataFrame(
{
"value": [
*bars.values(),
sum(bars.values()),
mediated,
float(data["y"].sel(geo="north").mean()),
]
},
index=[*bars, "sum of the waterfall bars", *two_dp_rows],
)
waterfall_summary.style.format("{:.3f}").format(
"{:.2f}", subset=pd.IndexSlice[two_dp_rows, :]
).set_caption("Waterfall completion (geo=north, original scale)")
| value | |
|---|---|
| intercept | 0.503 |
| channels (direct) | 0.641 |
| seasonality | 0.007 |
| sum of the waterfall bars | 1.151 |
| mediated, missing from the bars | 0.59 |
| observed target mean | 1.74 |
The bars sum to \(1.151\), while the observed target mean is \(1.74\), and the mediated term of \(0.59\) is missing.
Note that the posterior mediated term (\(0.59\)) is lower than the true mediated term in north (\(0.79\), see the next table). In contrast, the posterior intercept (\(0.503\)) is higher than the true intercept (\(0.35\)). Their sums are close (\(1.09\) against \(1.14\)), which suggests that the model shifts part of the mediated level into the intercept.
How much of the mediated term belongs to the channels? The next two tables come from the generative model for north, so they are comparable with each other but not with the posterior table above:
north = geos.index("north")
pd.DataFrame(
{
"target units": [
float(indirect_raw[:, north].mean()),
float(indirect_both[:, north].mean()),
]
},
index=["true mediated term (north)", "of which the channels created"],
).style.format("{:.2f}")
| target units | |
|---|---|
| true mediated term (north) | 0.79 |
| of which the channels created | 0.40 |
On the true model, the channels created \(0.40\) of the \(0.79\) target units per week of the mediated term in north. Next, we split the mean structural lower-funnel spend \(M^{*}\) in north into its four sources:
mstar_parts = {
"demand baseline": float(true["budget_dual_baseline"][north]),
"category demand (gamma C)": float(
true["budget_dual_parent_coef"][north] * category_demand[:, north].mean()
),
"lf budget (lambda B)": float(
true["budget_dual_budget_coef"][north] * lf_budget[:, north].mean()
),
}
mstar_mean = float(data["lf_spend_true"].sel(geo="north").mean())
mstar_parts["demand the channels created"] = mstar_mean - sum(mstar_parts.values())
pd.DataFrame(
{
"value (demand units)": list(mstar_parts.values()),
"share of M*": [v / mstar_mean for v in mstar_parts.values()],
},
index=list(mstar_parts),
).style.format({"value (demand units)": "{:.3f}", "share of M*": "{:.1%}"}).set_caption(
"Mean structural lower-funnel spend M* (north)"
)
| value (demand units) | share of M* | |
|---|---|---|
| demand baseline | 0.200 | 14.9% |
| category demand (gamma C) | 0.171 | 12.7% |
| lf budget (lambda B) | 0.164 | 12.2% |
| demand the channels created | 0.811 | 60.2% |
The channels created \(60.2\%\) of the mean structural lower-funnel spend in north. The demand baseline (\(14.9\%\)), the category demand term \(\gamma C\) (\(12.7\%\)) and the budget term \(\lambda B\) (\(12.2\%\)) make up the rest.
About a third of the target is missing from the waterfall. The missing \(0.59\) is the entire mediated term, and only one of its four sources, the demand that the channels created, belongs to media.
Compare like with like when you size it. The \(0.59\) is a posterior quantity, while the two values of the mediated-term table are generative truth. On the true model, the channels created about half of the mediated term, and the other three sources share the rest.
The \(M^{*}\) shares are in demand units, and the conversion from demand to target is concave. Therefore, the \(60.2\%\) share of the channels in \(M^{*}\) is not their share of the mediated term. Their share is the joint indirect effect, \(0.40\) of \(0.79\) target units on the true model.
If we read the bars as an attribution, we make two errors. We under-credit both channels by their mediated contributions, and we drop the demand baseline, the category index and the budget term, which make up the rest of the missing mass.
This is a limit of the framework. The DAG lives in your create_effect code, and only you can tell the library that an additive term means “demand that this channel created”.
The one built-in that follows the funnel, the incrementality module through the opt-in in incrementality_spec, answers the per-channel spend counterfactual, not this attribution question.
Funnel-Aware Budget Optimization#
Up to this point, we have measured effects. An MMM most often informs where the money should go, so we now turn the fitted model into a budget decision. In this section, we put the fitted funnel model through BudgetOptimizer, in-sample.
We optimize the weekly allocation across the six geo \(\times\) channel cells at a fixed total, turn the optimum into per-cell multipliers and re-weight the historical spend. Then we score the original and the re-weighted plans with the posterior and against the known truth.
Three choices make the exercise funnel-aware. If we get any of them wrong, the optimizer treats the funnel model as an ordinary model without raising an error.
First, the intervention reaches the effect. BudgetDualIndicatorEffect reads mmm.channel_data_scaled, so the candidate budgets of the optimizer also move the mediated path. For every candidate plan, the model recomputes latent demand, lower-funnel spend and the conversion term from the fitted graph.
An effect that read only its own pm.Data variables would not see the budgets. The optimizer would then price the media as if they created nothing downstream.
Second, the objective is a choice. The default response_variable, total_media_contribution_original_scale, comes from channel_contribution alone, so it scores the direct path only.
A model with mu_effects also registers total_response_original_scale, the total mean response on the original scale summed over dates and geos, which includes the mediated term. Its other parts (the intercept and the seasonality) do not depend on the budgets. Under the mean utility, its argmax is therefore the argmax of the direct plus the mediated response.
A risk-adjusted utility would also see the variance of the baseline, so it should use a contribution variable instead. We run both objectives on the same posterior, and the default objective shows what changes when we price the mediated path.
Third, the optimization window is the training window. We optimize over the training dates, so we can score the plans against the posterior and the known truth on the data that the model saw.
create_optimization_model() fills in everything that the effect needs over this window: category_demand, lf_budget and the two proxies keep their training values. We choose the end date so that the window plus the effective carryover of the model (effective_carryover_lags()) covers the training range.
Two conventions complete the setup. total_budget and the decision variables of the optimizer are weekly levels per cell, so the sum constraint holds on per-period budgets. We also pass budget_distribution_over_period equal to the historical spend pattern of each cell, instead of a uniform spread over the window.
This way, the optimizer scores the plan “multiplier times history”, which is the same object that we re-weight and evaluate later. The bounds are \(\pm 50\%\) around the historical weekly mean of each cell. The end of the window has a truncated carryover, so we find an in-sample optimum, not a steady-state one.
The next cell builds the optimization model and checks that the window covers the training dates:
carryover_lags = funnel.effective_carryover_lags() # l_max + the effect's declared lags
opt_model = funnel.create_optimization_model(
start_date=date_range[0], end_date=date_range[-(carryover_lags + 1)]
)
# Window plus effective carry-over is exactly the training range;
# create_optimization_model restores the observed values of every variable the
# effects read over those dates.
if len(opt_model.coords["date"]) != n_dates:
raise ValueError("the optimization window does not match the training dates")
We use the historical budget: the mean weekly spend of each cell, summed over the six cells. The week-by-week pattern stays fixed, so the optimizer moves money across cells, not across time. The next cell prints the total weekly budget and the mean weekly spend per cell:
spend_hist = data[channels].to_array("channel").transpose("date", "geo", "channel")
mean_plan = spend_hist.mean("date")
total_budget = float(mean_plan.sum())
budget_pattern = spend_hist / spend_hist.sum("date")
budget_bounds = xr.concat(
[0.5 * mean_plan, 1.5 * mean_plan],
dim=pd.Index(["lower", "upper"], name="bound"),
)
print(f"total weekly budget: {total_budget:.3f}")
mean_plan.to_pandas().style.format("{:.3f}").set_caption(
"Historical mean weekly spend per geo x channel"
)
total weekly budget: 4.316
| channel | tv_spend | social_spend |
|---|---|---|
| geo | ||
| north | 1.205 | 0.786 |
| south | 0.844 | 0.561 |
| west | 0.555 | 0.365 |
The total weekly budget is \(4.316\), and north has the largest spend in both channels (\(1.205\) for tv_spend and \(0.786\) for social_spend).
Next, a helper builds one BudgetOptimizer per objective, on the same optimization model and posterior:
def make_optimizer(response_variable: str) -> BudgetOptimizer:
"""Build a BudgetOptimizer on the funnel model for one objective.
Parameters
----------
response_variable : str
Name of the scalar model variable to maximize on average over the posterior.
Returns
-------
BudgetOptimizer
Optimizer over the six geo x channel weekly budgets, with the historical
spend pattern fixed over the training window.
"""
return BudgetOptimizer(
model=opt_model,
idata=funnel.idata,
num_periods=n_dates,
adstock_periods=0,
response_variable=response_variable,
budget_distribution_over_period=budget_pattern,
)
direct_label, funnel_label = "direct-only objective", "funnel-aware objective"
optimizers = {
direct_label: make_optimizer("total_media_contribution_original_scale"),
funnel_label: make_optimizer("total_response_original_scale"),
}
Optimizer Checks#
Before we optimize, we check that the optimizer sees the funnel. We evaluate and differentiate the responses that it works with through its public API, and three identities must hold:
Historical plan: at the historical plan, the direct, the mediated and the total response each equal the in-sample posterior mean of the matching variable. This confirms that the fixed pattern reproduces history and that the effect uses the training data.
Mediated marginals: the marginal of the mediated response with respect to the weekly spend of every cell is positive. Hence, the mediated path moves with the intervention, cell by cell.
Total marginals: the marginal of the total response equals the marginal of the direct response plus the marginal of the mediated response. Hence, the objective is the direct plus the mediated response plus a baseline that does not depend on the budget.
We let the optimizers move money only after all three checks pass.
post = funnel.idata.posterior
response_vars = {
"direct": "channel_contribution_original_scale",
"mediated": "budget_dual_effect_contribution_original_scale",
"total": "total_response_original_scale",
}
optimizer = optimizers[funnel_label]
flat = optimizer.optimization_variables.flat
def mean_over_posterior(name: str):
"""Sum a response graph over its dims and average it over the posterior draws."""
response = optimizer.extract_response_distribution(name)
reduce_dims = [dim for dim in response.dims if dim != "sample"]
if reduce_dims:
response = response.sum(dim=reduce_dims)
return response.mean(dim="sample")
def gradient_of(name: str):
"""Exact gradient of a response's posterior mean with respect to the decision vector.
xtensor Ops have no gradients yet, so the graph is lowered to plain tensor
Ops first, as the optimizer does for its own objective.
"""
scalar = rewrite_graph(
mean_over_posterior(name).values,
include=("lower_xtensor", "canonicalize", "stabilize"),
)
return pytensor.grad(scalar, flat)
# The gradient of the structural lower-funnel spend is the demand a cell creates
# per unit of its spend (the coefficient on D in M* is one); it splits the
# mediated marginal into demand created and demand converted further below.
gradient_vars = {**response_vars, "induced": "budget_dual_spend"}
evaluate = pytensor.function(
[flat], [mean_over_posterior(name) for name in response_vars.values()]
)
gradients = pytensor.function(
[flat], [gradient_of(name) for name in gradient_vars.values()]
)
def response_at(x: np.ndarray) -> pd.Series:
"""Posterior-mean direct, mediated and total response of a flat weekly plan."""
return pd.Series(
[float(value) for value in evaluate(x)], index=list(response_vars), dtype=float
)
x0 = optimizer.optimization_variables.pack(mean_plan)
at_history = response_at(x0)
in_sample = pd.Series(
{
"direct": post[response_vars["direct"]].sum(("date", "geo", "channel")).mean(),
"mediated": post[response_vars["mediated"]].sum(("date", "geo")).mean(),
"total": post[response_vars["total"]].mean(),
},
dtype=float,
)
np.testing.assert_allclose(at_history, in_sample, rtol=1e-6)
# The same cell labels serve every per-cell table below; everything is indexed
# by label, never by position in the flat vector.
mask = optimizer.budgets_to_optimize
cell_index = list(itertools.product(mask["geo"].values, mask["channel"].values))
cell_labels = [f"{geo}, {channel}" for geo, channel in cell_index]
def marginals_at(x: np.ndarray) -> pd.DataFrame:
"""Marginal response per unit of weekly spend in each cell, summed over the window.
Parameters
----------
x : np.ndarray
Flat weekly plan, in the optimizer's decision-vector layout.
Returns
-------
pd.DataFrame
One row per geo x channel cell; the direct, mediated and total response
and the induced lower-funnel spend, differentiated with respect to that
cell's weekly spend.
"""
columns = {}
for part, gradient in zip(gradient_vars, gradients(x), strict=True):
labelled = optimizer.optimization_variables.unpack(gradient)[
optimizer.channel_data_var
]
columns[part] = [
float(labelled.sel(geo=geo, channel=channel)) for geo, channel in cell_index
]
return pd.DataFrame(columns, index=cell_labels)
marginals = marginals_at(x0)
if not (marginals["mediated"] > 0).all():
raise RuntimeError("the mediated term does not increase with every cell's spend")
np.testing.assert_allclose(
marginals["total"], marginals["direct"] + marginals["mediated"], rtol=1e-6
)
if not all(opt.budgets_to_optimize.all() for opt in optimizers.values()):
raise RuntimeError("some geo x channel cell was left out of the optimization")
np.testing.assert_array_equal(
x0, optimizers[direct_label].optimization_variables.pack(mean_plan)
)
print(f"direct response at the historical plan: {at_history['direct']:.2f}")
print(f"mediated response at the historical plan: {at_history['mediated']:.2f}")
print(
f"total response at the historical plan: {at_history['total']:.2f}"
" (the funnel-aware objective; it includes the baseline)"
)
direct response at the historical plan: 223.98
mediated response at the historical plan: 192.65
total response at the historical plan: 536.22 (the funnel-aware objective; it includes the baseline)
All three checks pass. At the historical plan, the posterior mean direct response is \(223.98\) and the mediated response is \(192.65\), summed over the window. The total response of \(536.22\) also includes the baseline.
Marginal Response#
The marginals are the quantity that the decision depends on: the marginal response of each cell per extra unit of weekly spend, and its split between the direct and the mediated path. An interior optimum equalizes this marginal across cells.
Therefore, the ordering of the total column says where the funnel-aware optimizer moves money, and the ordering of the direct column says where the default objective moves it. The share column shows how much of the marginal the default objective sees in each cell.
Because the data are synthetic, we can also compute the same marginals on the generative model, clamped to the true parameters, with central differences of \(\pm 0.5\%\) around the historical plan. The two halves of the table check that the optimizer sees the funnel with about the right slope, cell by cell:
def true_response_parts(multiplier: xr.DataArray) -> pd.Series:
"""Direct and media-created mediated response, and lower-funnel spend, of a plan on the truth.
Parameters
----------
multiplier : xr.DataArray
Per-cell multipliers with dims ``("geo", "channel")`` applied to the
historical media plan.
Returns
-------
pd.Series
``direct`` and ``mediated`` response summed over the window (the mediated
part net of the both-channels-off pass), and the structural lower-funnel
spend ``lf_spend`` the plan induces, all on the clamped generative model.
"""
media = (
media_raw
* multiplier.sel(geo=geos, channel=channels)
.transpose("geo", "channel")
.to_numpy()[None]
)
names = [
"channel_contribution_original_scale",
"budget_dual_effect_contribution_original_scale",
"budget_dual_spend",
]
# Every parameter is clamped, so this is one deterministic forward pass;
# the historical media are restored whatever happens in between.
with gen.model:
try:
pm.set_data({"channel_data": media})
values = pm.draw([gen.model[name] for name in names])
finally:
pm.set_data({"channel_data": media_raw})
draw = {
name: xr.DataArray(value, dims=gen.model.named_vars_to_dims[name])
for name, value in zip(names, values, strict=True)
}
mediated = (
draw["budget_dual_effect_contribution_original_scale"]
.transpose("date", "geo")
.to_numpy()
- cf_gen["both"]["budget_dual_effect_contribution_original_scale"]
)
return pd.Series(
{
"direct": float(draw["channel_contribution_original_scale"].sum()),
"mediated": float(mediated.sum()),
"lf_spend": float(draw["budget_dual_spend"].sum()),
},
dtype=float,
)
true_at_history = true_response_parts(xr.ones_like(mean_plan))
# Central differences of +/-0.5% around the historical plan, per unit of the
# cell's mean weekly spend, indexed by the same labels as the posterior gradients.
true_marginal_rows = {}
for label, (geo, channel) in zip(cell_labels, cell_index, strict=True):
up, down = xr.ones_like(mean_plan), xr.ones_like(mean_plan)
up.loc[{"geo": geo, "channel": channel}] = 1.005
down.loc[{"geo": geo, "channel": channel}] = 0.995
per_unit = 0.01 * float(mean_plan.sel(geo=geo, channel=channel))
true_marginal_rows[label] = (
true_response_parts(up) - true_response_parts(down)
) / per_unit
true_marginals = pd.DataFrame(true_marginal_rows).T.rename(
columns={"lf_spend": "induced"}
)
true_marginals["total"] = true_marginals["direct"] + true_marginals["mediated"]
np.testing.assert_allclose(gen.model["channel_data"].get_value(), media_raw)
marginal_table = pd.concat(
{
"posterior": marginals[["direct", "mediated", "total"]].assign(
**{
"share the default objective sees": marginals["direct"]
/ marginals["total"]
}
),
"truth": true_marginals[["direct", "mediated", "total"]],
},
axis=1,
)
marginal_table.style.format("{:.1f}").format(
"{:.0%}", subset=[("posterior", "share the default objective sees")]
).set_caption(
"Marginal response per extra unit of weekly spend in a cell, summed over the window"
" (posterior: exact gradient; truth: central differences)"
)
| posterior | truth | ||||||
|---|---|---|---|---|---|---|---|
| direct | mediated | total | share the default objective sees | direct | mediated | total | |
| north, tv_spend | 25.0 | 10.9 | 35.9 | 70% | 23.2 | 11.1 | 34.3 |
| north, social_spend | 19.5 | 12.5 | 32.0 | 61% | 20.5 | 12.9 | 33.5 |
| south, tv_spend | 24.6 | 19.0 | 43.5 | 56% | 27.8 | 16.0 | 43.8 |
| south, social_spend | 38.4 | 14.9 | 53.3 | 72% | 26.0 | 17.4 | 43.4 |
| west, tv_spend | 33.3 | 20.5 | 53.8 | 62% | 32.5 | 19.6 | 52.1 |
| west, social_spend | 34.2 | 18.3 | 52.5 | 65% | 31.9 | 20.4 | 52.4 |
The default objective sees \(56\%\) to \(72\%\) of the marginal response, depending on the cell. In five of six cells, the posterior total marginal is within about \(5\%\) of the truth.
The exception is south, social_spend, where the posterior total marginal is \(53.3\) against a true \(43.4\), because its direct marginal is too high (\(38.4\) against \(26.0\)).
Optimal Allocations#
Both runs start from the historical plan and share the bounds and the total. Each optimizer maximizes its own objective, so the funnel-aware optimum cannot score below the direct-only optimum on the total response. Likewise, the direct-only optimum cannot score below the funnel-aware optimum on the direct response.
The next cell runs both optimizations, checks both claims and checks that every multiplier stays strictly inside the \(\pm 50\%\) box. The table shows the optimal weekly spend of each cell as a multiple of its historical mean:
minimize_kwargs = {"options": {"ftol": 1e-8, "maxiter": 2_000}}
allocations = {}
for label, opt in optimizers.items():
result = opt.allocate_budget(
total_budget=total_budget,
budget_bounds=budget_bounds,
x0=mean_plan,
minimize_kwargs=minimize_kwargs,
)
if not result.scipy_result.success:
raise RuntimeError(f"{label}: {result.scipy_result.message}")
np.testing.assert_allclose(result.budgets.sum(), total_budget, rtol=1e-6)
allocations[label] = result
print(
f"{label}: {result.scipy_result.message} ({result.scipy_result.nit} iterations)"
)
multipliers = xr.concat(
[allocations[label].budgets / mean_plan for label in allocations],
dim=pd.Index(list(allocations), name="plan"),
)
# Strictly interior: every multiplier at least 2% of the box's width away from
# the bounds, derived from the bounds actually in force.
bound_ratio = budget_bounds / mean_plan
lower_ratio = float(bound_ratio.sel(bound="lower").max())
upper_ratio = float(bound_ratio.sel(bound="upper").min())
margin = 0.02 * (upper_ratio - lower_ratio)
if (
float(multipliers.min()) < lower_ratio + margin
or float(multipliers.max()) > upper_ratio - margin
):
raise RuntimeError("an optimal multiplier sits at the budget bounds")
# Each optimizer wins on its own objective, up to a slack tied to the solver's ftol.
at_optimum = {
label: response_at(allocations[label].scipy_result.x) for label in allocations
}
ftol = minimize_kwargs["options"]["ftol"]
for winner, loser, part in [
(funnel_label, direct_label, "total"),
(direct_label, funnel_label, "direct"),
]:
slack = 10 * ftol * max(abs(at_optimum[winner][part]), 1.0)
if at_optimum[winner][part] < at_optimum[loser][part] - slack:
raise RuntimeError(f"the {winner} optimum loses on its own objective")
comparison = pd.DataFrame(
{
label: multipliers.sel(plan=label).stack(cell=("geo", "channel")).to_pandas()
for label in allocations
}
)
comparison.index = [f"{geo}, {channel}" for geo, channel in comparison.index]
comparison["difference (funnel - direct)"] = (
comparison[funnel_label] - comparison[direct_label]
)
comparison["same direction"] = np.sign(comparison[funnel_label] - 1) == np.sign(
comparison[direct_label] - 1
)
comparison.style.format("{:.3f}", subset=comparison.columns[:3]).set_caption(
"Optimal weekly spend as a multiple of the historical mean"
)
direct-only objective: Optimization terminated successfully (8 iterations)
funnel-aware objective: Optimization terminated successfully (11 iterations)
| direct-only objective | funnel-aware objective | difference (funnel - direct) | same direction | |
|---|---|---|---|---|
| north, tv_spend | 0.936 | 0.851 | -0.085 | True |
| north, social_spend | 0.838 | 0.843 | 0.004 | True |
| south, tv_spend | 0.918 | 1.040 | 0.121 | False |
| south, social_spend | 1.207 | 1.161 | -0.046 | True |
| west, tv_spend | 1.195 | 1.228 | 0.033 | True |
| west, social_spend | 1.134 | 1.144 | 0.011 | True |
Both optimizers converge (\(8\) and \(11\) iterations), and every multiplier stays inside the box, between \(0.838\) and \(1.228\). The two objectives agree on the direction of five of six moves.
They disagree on south, tv_spend: the direct-only objective cuts it to \(0.918\), while the funnel-aware objective grows it to \(1.040\). The other large difference is north, tv_spend (\(0.936\) against \(0.851\)).
The marginal table explains why. A cell with a modest direct marginal can still have a large mediated marginal, and pricing that path changes its recommendation. The figure shows the historical mean and both optima side by side:
The funnel-aware plan cuts north, tv_spend more than the direct-only plan and grows south, tv_spend instead of cutting it. Both plans cut north, social_spend and grow south, social_spend and both west cells.
Marginal Equalization#
At an interior optimum with one budget constraint, the optimizer equalizes the marginal that it sees per unit of weekly spend across the cells. This common value is the multiplier of the budget constraint, and the next cell reads it off both optima. The default objective equalizes the direct marginal, and the funnel-aware objective equalizes the total marginal.
Two consequences follow. First, the position of the historical marginal of a cell relative to this level decides whether the optimizer cuts or grows the cell. With diminishing returns, a marginal below the level rises only if the cell is cut, and a marginal above the level falls only if the cell grows.
The shared demand pool couples the two channels of a geo, so this rule holds one cell at a time and is not an exact identity. Second, how far a cell moves also depends on its curvature.
As a result, we can only compare the multipliers of cells that share their curvature parameters and their position on the response curve. The north and south tv_spend cells are such a pair, because the carryover and saturation parameters of tv_spend are pooled across geos. We check their position on the curve below.
optimum_marginals = {
label: marginals_at(allocations[label].scipy_result.x) for label in allocations
}
# The equalization is checked at a relative 1e-2, the sharpness a first-order
# condition supports at the solver tolerance used above.
equalization_rtol = 1e-2
equalized = {}
for label, part in [(direct_label, "direct"), (funnel_label, "total")]:
values = optimum_marginals[label][part]
spread = float(values.max() - values.min()) / float(values.mean())
if spread > equalization_rtol:
raise RuntimeError(f"{label}: the {part} marginals are not equalized")
equalized[label] = float(values.mean())
marginals_by_plan = pd.concat(
{
"historical plan": marginals[["direct", "mediated", "total"]],
**{
label: pd.concat(
[
comparison.loc[cell_labels, label].rename("multiplier"),
optimum_marginals[label][["direct", "total"]],
],
axis=1,
)
for label in allocations
},
},
axis=1,
)
marginals_by_plan.style.format("{:.1f}").format(
"{:.3f}", subset=[(label, "multiplier") for label in allocations]
).set_caption(
"Marginal response per unit of weekly spend at the historical plan and at each"
f" optimum. The direct-only optimum equalizes the direct marginal at"
f" {equalized[direct_label]:.1f}, the funnel-aware optimum the total marginal at"
f" {equalized[funnel_label]:.1f}"
)
| historical plan | direct-only objective | funnel-aware objective | |||||||
|---|---|---|---|---|---|---|---|---|---|
| direct | mediated | total | multiplier | direct | total | multiplier | direct | total | |
| north, tv_spend | 25.0 | 10.9 | 35.9 | 0.936 | 26.6 | 38.5 | 0.851 | 28.8 | 41.6 |
| north, social_spend | 19.5 | 12.5 | 32.0 | 0.838 | 26.6 | 41.6 | 0.843 | 26.3 | 41.6 |
| south, tv_spend | 24.6 | 19.0 | 43.5 | 0.918 | 26.6 | 46.6 | 1.040 | 23.6 | 41.6 |
| south, social_spend | 38.4 | 14.9 | 53.3 | 1.207 | 26.6 | 39.1 | 1.161 | 28.9 | 41.6 |
| west, tv_spend | 33.3 | 20.5 | 53.8 | 1.195 | 26.6 | 43.2 | 1.228 | 25.6 | 41.6 |
| west, social_spend | 34.2 | 18.3 | 52.5 | 1.134 | 26.6 | 42.4 | 1.144 | 26.0 | 41.6 |
The direct-only optimum equalizes the direct marginal at \(26.6\) in every cell, and the funnel-aware optimum equalizes the total marginal at \(41.6\).
The next figure shows the reading rule in one picture: the marginal of each cell at the historical plan, split into its direct and mediated parts, against the two levels. A cell whose direct bar ends left of the first line is cut by the direct-only optimizer. A cell whose full bar ends left of the second line is cut by the funnel-aware optimizer:
Two cells sit below both levels (north, both channels), and both objectives cut them. Three cells sit above both levels (south, social_spend and both west cells), and both objectives grow them.
south, tv_spend is the only cell between the levels: its direct marginal is below the direct level, and its total marginal is above the total level. Therefore, its direction flips between the two objectives.
Why do the mediated marginals differ between cells? The mediated marginal is the product of two factors.
The first factor is the demand that a unit of spend in the cell creates over the window. Because \(M^{*} = D + \lambda B\) has a unit coefficient on \(D\), this demand equals the lower-funnel spend that the cell induces per unit.
The second factor is how much response a unit of this demand buys. It tracks the slope of the shared conversion curve at the level of structural lower-funnel spend of the geo. Strictly, it is a ratio of two window-summed marginals, so it is a demand-weighted average of this slope over the window.
The next cell prints the mean structural lower-funnel spend per geo and the position of each cell on its own spend scale. Then it splits each mediated marginal into the two factors:
def conversion_split(table: pd.DataFrame) -> pd.DataFrame:
"""Split a mediated marginal into demand created per unit of spend and its conversion."""
return pd.DataFrame(
{
"induced lower-funnel spend per unit": table["induced"],
"response per unit of induced spend": table["mediated"] / table["induced"],
"mediated": table["mediated"],
}
)
mediated_split = pd.concat(
{
"posterior": conversion_split(marginals),
"truth": conversion_split(true_marginals),
},
axis=1,
)
lf_level = pd.DataFrame(
{
"posterior mean": post["budget_dual_spend"]
.mean(("chain", "draw", "date"))
.to_pandas(),
"truth": data["lf_spend_true"].mean("date").to_pandas(),
}
)
print("structural lower-funnel spend per week, mean over the window:")
print(lf_level.round(3).to_string())
# The channel scale is each cell's maximum week of spend, so a unit of raw spend
# is a smaller step up the (pooled) response curves of a larger cell.
channel_scale = funnel.idata.constant_data["channel_scale"].transpose("geo", "channel")
spend_position = pd.concat(
{
"channel scale (max weekly spend)": channel_scale.to_pandas(),
"mean week as a share of the scale": (mean_plan / channel_scale).to_pandas(),
},
axis=1,
)
print("\nwhere each cell's mean week sits on its own scale:")
print(spend_position.round(3).to_string())
mediated_split.style.format("{:.1f}").format(
"{:.2f}",
subset=[
(half, "response per unit of induced spend") for half in ("posterior", "truth")
],
).set_caption(
"The mediated marginal as demand created times demand converted, per cell"
)
structural lower-funnel spend per week, mean over the window:
posterior mean truth
geo
north 1.350 1.347
south 1.031 1.029
west 0.776 0.777
where each cell's mean week sits on its own scale:
channel scale (max weekly spend) mean week as a share of the scale
channel tv_spend social_spend tv_spend social_spend
geo
north 2.607 1.295 0.462 0.607
south 1.791 0.997 0.471 0.562
west 1.046 0.611 0.531 0.597
| posterior | truth | |||||
|---|---|---|---|---|---|---|
| induced lower-funnel spend per unit | response per unit of induced spend | mediated | induced lower-funnel spend per unit | response per unit of induced spend | mediated | |
| north, tv_spend | 36.2 | 0.30 | 10.9 | 34.5 | 0.32 | 11.1 |
| north, social_spend | 41.1 | 0.30 | 12.5 | 38.3 | 0.34 | 12.9 |
| south, tv_spend | 45.0 | 0.42 | 19.0 | 41.6 | 0.38 | 16.0 |
| south, social_spend | 35.1 | 0.42 | 14.9 | 44.1 | 0.39 | 17.4 |
| west, tv_spend | 53.4 | 0.38 | 20.5 | 50.7 | 0.39 | 19.6 |
| west, social_spend | 47.6 | 0.38 | 18.3 | 52.2 | 0.39 | 20.4 |
The two channels of a geo have nearly the same conversion rate: \(0.30\) in north, \(0.42\) in south and \(0.38\) in west. This makes the conversion rate a property of the geo.
The mean week of tv_spend sits at nearly the same point of its scaled curve in north and south (\(0.462\) and \(0.471\) of the maximum), so these two cells are comparable.
Let’s read the north and south tv_spend cells off the tables. The direct-only objective sees nearly equal direct marginals (\(25.0\) and \(24.6\) per unit of weekly spend). Both are below the level of \(26.6\), so it cuts both cells (multipliers \(0.936\) and \(0.918\)).
The funnel-aware objective sees total marginals of \(35.9\) and \(43.5\). It cuts north, tv_spend to \(0.851\) and grows south, tv_spend to \(1.040\), until both total marginals reach \(41.6\). The true total marginals (\(34.3\) and \(43.8\)) support this distinction.
The mediated marginal is \(19.0\) in south against \(10.9\) in north, about \(1.74\) times larger (\(16.0\) against \(11.1\) in truth). Two factors produce the gap.
First, south creates more demand per unit of tv_spend: \(45.0\) units of lower-funnel spend over the window against \(36.2\) in north (truth: \(41.6\) against \(34.5\)). A raw unit of tv_spend is a smaller step up the north curve, because its maximum weekly spend is \(2.607\) instead of \(1.791\).
Second, south converts the induced spend more steeply: \(0.42\) target units per induced unit against \(0.30\) in north (truth: \(0.38\) against \(0.32\)). north has the highest structural lower-funnel spend (about \(1.35\) per week, against \(1.03\) in south and \(0.78\) in west), so it sits further up the shared saturating curve.
The posterior sat_contribution_beta in south is \(1.108\) against a true \(0.90\), covered but high, which adds to its conversion slope. Both factors point in the same direction in the DGP and in the posterior, but the posterior makes their combined gap larger.
In-Sample Evaluation#
We turn each optimum into per-cell multipliers, re-weight the historical spend and score the plans with the posterior. The total spend stays the same by construction, because the sum constraint holds on weekly levels and the pattern is fixed. The assertion below confirms this.
The one funnel-specific choice is what to score. The quantity that a budget decision depends on is the media-driven response: the direct contributions plus the part of the mediated term that the media created.
We follow the both-channels-off convention of the true decomposition above, which gives the joint indirect effect and not the leave-one-out sum. Per posterior draw, we take the mediated term under a plan minus the mediated term with both channels switched off. Hence, the demand baseline, the category index and the lower-funnel budget cancel.
This needs one more posterior predictive pass, with both channels at zero and the exogenous data of the effect at their training values. Then the both-off pass differs from a plan only in the channel spend.
Under the mean utility, this yardstick and the objective have the same argmax, because they differ by a constant per draw. We recompute both from the model instead of reading them off the optimizer. The check passes when the value of each optimizer at its optimum equals the posterior predictive value of the re-weighted plan:
def sample_plan(X: pd.DataFrame) -> xr.Dataset:
"""Posterior-predictive contributions of the funnel model under a spend plan.
Parameters
----------
X : pd.DataFrame
Long frame with the date, geo and upper-funnel channel columns.
Returns
-------
xr.Dataset
Direct, mediated and total response on the original scale, and the
structural lower-funnel spend, per chain and draw.
"""
return funnel.sample_posterior_predictive(
X,
var_names=[*response_vars.values(), "budget_dual_spend"],
combined=False,
extend_idata=False,
progressbar=False,
random_seed=rng,
)
def plan_frame(multiplier: xr.DataArray) -> pd.DataFrame:
"""Re-weight the historical spend by per-cell multipliers into a long frame."""
return (spend_hist * multiplier).to_dataset("channel").to_dataframe().reset_index()
plan_frames = {label: plan_frame(multipliers.sel(plan=label)) for label in allocations}
for X_plan in plan_frames.values():
np.testing.assert_allclose(
X_plan[channels].to_numpy().sum(), float(spend_hist.sum()), rtol=1e-6
)
X_zero = X_upper.copy()
X_zero[channels] = 0.0
pp = {label: sample_plan(X_plan) for label, X_plan in plan_frames.items()}
pp_zero = sample_plan(X_zero)
pp["historical"] = xr.Dataset(
{var: post[var] for var in [*response_vars.values(), "budget_dual_spend"]}
)
for samples in pp.values():
xr.align(samples, pp_zero, join="exact")
mediated_zero = pp_zero[response_vars["mediated"]]
if float(mediated_zero.mean()) >= float(
pp["historical"][response_vars["mediated"]].mean()
):
raise RuntimeError("switching both channels off did not reduce the mediated term")
# The both-channels-off yardstick is the optimizer's own graph at zero spend.
np.testing.assert_allclose(
mediated_zero.sum(("date", "geo")).mean(),
response_at(np.zeros_like(x0))["mediated"],
rtol=1e-6,
)
def media_response(samples: xr.Dataset) -> xr.DataArray:
"""Media-driven response per date: direct plus media-created mediated term."""
direct = samples[response_vars["direct"]].sum(("geo", "channel"))
mediated = (samples[response_vars["mediated"]] - mediated_zero).sum("geo")
return direct + mediated
response = {label: media_response(samples) for label, samples in pp.items()}
response_total = {label: r.sum("date") for label, r in response.items()}
direct_total = {
label: samples[response_vars["direct"]].sum(("date", "geo", "channel"))
for label, samples in pp.items()
}
lf_spend_total = {
label: samples["budget_dual_spend"].sum(("date", "geo"))
for label, samples in pp.items()
}
# Loop closure: each optimizer's reported optimum equals the posterior-predictive
# recomputation of its own objective under the re-weighted plan.
recomputed = {
funnel_label: pp[funnel_label][response_vars["total"]].mean(),
direct_label: direct_total[direct_label].mean(),
}
for label, part in [(funnel_label, "total"), (direct_label, "direct")]:
np.testing.assert_allclose(at_optimum[label][part], recomputed[label], rtol=1e-6)
np.testing.assert_allclose(
-allocations[label].scipy_result.fun, at_optimum[label][part], rtol=1e-6
)
All checks pass. Because the data are synthetic, we can also price each plan on the generative model, with the same helper that produced the true marginals above. One forward pass per plan gives the truth against which we check the posterior uplifts:
true_parts = {"historical": true_at_history}
for label in allocations:
true_parts[label] = true_response_parts(multipliers.sel(plan=label))
np.testing.assert_allclose(gen.model["channel_data"].get_value(), media_raw)
true_response = {
label: parts["direct"] + parts["mediated"] for label, parts in true_parts.items()
}
true_direct_hist = float(data["direct_true"].sum())
np.testing.assert_allclose(true_parts["historical"]["direct"], true_direct_hist)
The next three figures show the result of the funnel-aware optimization. First, the plans themselves: the historical spend against the funnel-aware optimum, per channel and geo (solid against dashed lines):
Each dashed line is the historical line scaled by its multiplier, from \(0.84\) for north, social_spend to \(1.23\) for west, tv_spend. The week-by-week pattern of every cell stays the same.
Next, we plot the media-driven response over time under the original plan and the funnel-aware plan, summed over geos and channels, with \(94\%\) HDIs:
The two bands overlap almost completely, and the optimized mean is only slightly above the original mean in most weeks. The reallocation is a small change to a plan that was already close to the posterior optimum within the \(\pm 50\%\) box.
Finally, the uplift: the posterior distribution of the ratio of the media-driven response under each optimized plan to that under the historical plan. The dashed line is the true uplift from the generative model. We score the direct-only plan with the same funnel-aware yardstick, so the two panels compare like with like:
uplift = xr.Dataset(
{
label: response_total[label] / response_total["historical"] - 1
for label in allocations
}
)
true_uplift = {
label: true_response[label] / true_response["historical"] - 1
for label in allocations
}
pc = azp.plot_dist(
uplift,
ci_kind="hdi",
ci_prob=0.94,
col_wrap=1,
figure_kwargs={"figsize": (10, 9), "sharex": True},
)
fig = pc.viz["figure"].item()
xlim = (float(uplift.to_array().min()), float(uplift.to_array().max()))
for label in allocations:
ax = pc.viz["plot"][label].item()
ax.axvline(0, color="C3", linestyle=":", linewidth=2, label="no change")
ax.axvline(
true_uplift[label],
color="black",
linestyle="--",
linewidth=2,
label="true uplift",
)
ax.set_xlim(xlim)
ax.xaxis.set_major_formatter(plt.matplotlib.ticker.PercentFormatter(xmax=1))
ax.set_title(f"{label}: plan re-weighted by its optimum")
ax.legend(loc="upper right")
fig.suptitle(
"Uplift in media-driven response after optimization", fontsize=18, fontweight="bold"
);
Both posterior uplifts are positive in most draws, and both true uplifts fall inside the posterior distributions. The tables below give the numbers.
The first table also shows what the direct-only optimizer itself expects to gain, on its own objective. Every uplift here is the mean of the per-draw ratios, while the optimizer objective is a posterior mean, so the reported optimum of the optimizer is not a row of these tables:
def summarise(da: xr.DataArray) -> tuple[float, float, float]:
"""Return the posterior mean and 94% HDI bounds of a scalar posterior quantity."""
lo, hi = hdi_bounds(da)
return float(da.mean()), float(lo), float(hi)
rows = {}
for label in allocations:
mean, lo, hi = summarise(uplift[label])
rows[label] = {
"posterior mean uplift": mean,
"94% HDI lower": lo,
"94% HDI upper": hi,
"true uplift": true_uplift[label],
}
rows[direct_label]["self-assessed uplift (direct objective)"] = float(
(direct_total[direct_label] / direct_total["historical"] - 1).mean()
)
uplift_table = pd.DataFrame(rows).T
uplift_table.style.format("{:+.2%}", na_rep="").set_caption(
"Uplift in media-driven response versus the historical plan"
)
| posterior mean uplift | 94% HDI lower | 94% HDI upper | true uplift | self-assessed uplift (direct objective) | |
|---|---|---|---|---|---|
| direct-only objective | +0.68% | +0.06% | +1.33% | +0.28% | +0.81% |
| funnel-aware objective | +0.82% | +0.23% | +1.45% | +0.54% |
The funnel-aware plan raises the posterior media-driven response by \(0.82\%\) (\(94\%\) HDI \(0.23\%\) to \(1.45\%\)), and the direct-only plan raises it by \(0.68\%\) (\(0.06\%\) to \(1.33\%\)). On its own objective, the direct-only optimizer expects \(0.81\%\).
Next, the absolute uplifts with their true direct and mediated parts, and the change in the lower-funnel spend that each plan induces:
parts_rows = {}
for label in allocations:
delta_true = true_parts[label] - true_parts["historical"]
parts_rows[label] = {
"absolute uplift (posterior mean)": float(
(response_total[label] - response_total["historical"]).mean()
),
"true absolute uplift": delta_true["direct"] + delta_true["mediated"],
"true direct part": delta_true["direct"],
"true mediated part": delta_true["mediated"],
"induced lower-funnel spend (posterior mean)": float(
(lf_spend_total[label] / lf_spend_total["historical"] - 1).mean()
),
"induced lower-funnel spend (true)": (
true_parts[label]["lf_spend"] / true_parts["historical"]["lf_spend"] - 1
),
}
parts_table = pd.DataFrame(parts_rows).T
parts_table.style.format("{:+.2f}").format(
"{:+.2%}", subset=parts_table.columns[-2:]
).set_caption("What the uplift is made of, and the lower-funnel spend the plan induces")
| absolute uplift (posterior mean) | true absolute uplift | true direct part | true mediated part | induced lower-funnel spend (posterior mean) | induced lower-funnel spend (true) | |
|---|---|---|---|---|---|---|
| direct-only objective | +2.24 | +0.91 | +0.21 | +0.70 | -0.12% | +0.21% |
| funnel-aware objective | +2.71 | +1.74 | +0.68 | +1.06 | +0.13% | +0.38% |
In truth, most of each gain travels through the mediated path: \(1.06\) of \(1.74\) target units for the funnel-aware plan, and \(0.70\) of \(0.91\) for the direct-only plan. In truth, both plans also induce more lower-funnel spend (\(+0.38\%\) and \(+0.21\%\)).
The posterior gets the direction of the induced spend right for the funnel-aware plan (\(+0.13\%\)) but not for the direct-only plan (\(-0.12\%\)).
Finally, the media-driven response of the historical plan, under the default objective and under the full yardstick:
level_table = pd.DataFrame(
{
"posterior mean": [
float(direct_total["historical"].mean()),
float(response_total["historical"].mean()),
],
"true value": [true_direct_hist, true_response["historical"]],
},
index=["direct only (default objective)", "direct + media-created mediated"],
)
level_table.loc["share the default objective sees"] = (
level_table.iloc[0] / level_table.iloc[1]
)
level_table.style.format("{:.1f}").format(
"{:.1%}", subset=pd.IndexSlice[["share the default objective sees"], :]
).set_caption("Media-driven response of the historical plan, summed over the window")
| posterior mean | true value | |
|---|---|---|
| direct only (default objective) | 224.0 | 197.1 |
| direct + media-created mediated | 333.7 | 324.2 |
| share the default objective sees | 67.1% | 60.8% |
Optimization Results#
Let’s start with the level of the historical plan. The default objective values the direct media response at \(224.0\) target units, while the direct plus media-created mediated response is \(333.7\). In other words, the default objective sees \(67.1\%\) of what media do under the posterior, and \(60.8\%\) in truth (\(197.1\) of \(324.2\)).
These shares use the both-channels-off counterfactual, so they exclude the demand from the baseline, the category index and the lower-funnel budget. For a budget decision, the marginal table is more relevant: depending on the cell, the default objective sees \(56\%\) to \(72\%\) of the marginal response.
The reallocations are modest. The funnel-aware plan raises the posterior media-driven response by \(0.82\%\), and the direct-only plan raises it by \(0.68\%\) on the same yardstick.
In truth, the uplifts are \(+0.54\%\) and \(+0.28\%\), that is, \(1.74\) and \(0.91\) target units, so the funnel-aware plan is worth more. Both uplifts are gross of the extra lower-funnel spend that the plans induce, which the ROAS section prices and this optimizer does not.
The true gains fall inside the posterior \(94\%\) intervals, in their lower halves. The posterior overstates the absolute gains: \(2.71\) against \(1.74\) for the funnel-aware plan, and \(2.24\) against \(0.91\) for the direct-only plan.
Three effects can contribute to this overstatement: selection on the posterior mean, the over-credited south, social_spend marginal, and extrapolation beyond the largest fitted weeks of some cells. One synthetic dataset cannot separate them.
This is an in-sample optimum with a truncated carryover at the end of the window. It compares two objectives on the same fitted funnel model, not the plan that a model without the mediator would recommend.
Business Value: on this dataset, the funnel-aware objective gives a plan with a larger true gain than the default objective, and it reverses the move of one cell (south, tv_spend). The default objective misses \(28\%\) to \(44\%\) of the marginal response that the decision depends on, depending on the cell.
Lower-Funnel Co-Optimization#
Up to this point, we treated lf_budget as given: the optimizer moved media across the six geo \(\times\) channel cells, and the lower-funnel budget kept its historical values. Lower-funnel spend is also money, so a planner may want to split one pot between the upper funnel and the lower funnel.
The argument spend_vars names additional monetary pm.Data nodes to optimize jointly with the media. Each one is decided over its own non-date dimensions, so lf_budget, which has dimensions (date, geo), becomes one budget per geo. Its spend counts toward the same budget-sum constraint as the media budgets.
Spend variables differ from optimizable_vars, which are for quantities in their own units, for example a discount depth or a frequency. Those quantities sit outside the sum constraint, while spend variables count inside it.
This run differs from the runs above in three ways:
Pot: we add the historical lower-funnel budget to the total. Otherwise, the joint plan would have to fund the lower funnel from nothing.
Time profile: media follow the historical weekly pattern through
budget_distribution_over_period, while a spend variable is spread uniformly across the window unless you pass a pattern for it.Bounds: we pass no
budget_bounds, so every budget can move between \(0\) and the total, as theUserWarningin the output states. The media cells also lose the \(\pm 50\%\) box of the runs above.
The next cell runs the joint optimization and checks that media and lower-funnel spend add up to the joint pot:
# The joint pot is what was actually spent on both, so the two plans are
# comparable: the optimizer chooses the split rather than being handed one.
lf_mean_weekly = data["lf_budget"].mean("date")
joint_budget = total_budget + float(lf_mean_weekly.sum())
joint_optimizer = BudgetOptimizer(
model=opt_model,
idata=funnel.idata,
num_periods=n_dates,
adstock_periods=0,
response_variable="total_response_original_scale",
budget_distribution_over_period=budget_pattern,
spend_vars=["lf_budget"],
)
joint_result = joint_optimizer.allocate_budget(
total_budget=joint_budget, minimize_kwargs=minimize_kwargs
)
if not joint_result.scipy_result.success:
raise RuntimeError(joint_result.scipy_result.message)
# `spend_var_allocations` is the monetary subset of `optimized_vars`:
# levers live there too.
media_spend = float(joint_result.budgets.sum())
lower_funnel_spend = sum(
float(allocation.sum())
for allocation in joint_result.spend_var_allocations.values()
)
np.testing.assert_allclose(media_spend + lower_funnel_spend, joint_budget, rtol=1e-6)
split = pd.DataFrame(
{
"historical": [total_budget, float(lf_mean_weekly.sum())],
"joint optimum": [media_spend, lower_funnel_spend],
},
index=pd.Index(["upper funnel (media)", "lower funnel (lf_budget)"], name="pot"),
)
split.loc["total"] = split.sum()
split.style.format("{:.3f}").set_caption(
f"Weekly spend split, one pot of {joint_budget:.3f}"
)
/var/folders/cm/3dzy9rdd5s3672z0s1brjkvh0000gn/T/ipykernel_90584/3328158137.py:16: UserWarning: No budget bounds provided. Using default bounds (0, total_budget) for each channel.
joint_result = joint_optimizer.allocate_budget(
| historical | joint optimum | |
|---|---|---|
| pot | ||
| upper funnel (media) | 4.316 | 5.026 |
| lower funnel (lf_budget) | 0.710 | 0.000 |
| total | 5.026 | 5.026 |
The joint optimum moves the whole lower-funnel budget into media. The weekly lower-funnel spend falls from \(0.710\) to \(0\), and the media spend rises from \(4.316\) to \(5.026\). The same split, per geo:
pd.concat(
{
"historical": lf_mean_weekly.to_pandas(),
"joint optimum": joint_result.spend_var_allocations["lf_budget"].to_pandas(),
},
axis=1,
).style.format("{:.3f}").set_caption(
"Lower-funnel weekly budget per geo: given, and decided"
)
| historical | joint optimum | |
|---|---|---|
| geo | ||
| north | 0.329 | 0.000 |
| south | 0.249 | 0.000 |
| west | 0.133 | 0.000 |
The lower-funnel budget falls to \(0\) in every geo. The constraint check of the previous cell confirms that the optimizer adds media and lower-funnel spend together and matches them against one total. The default sum constraint adds every decision variable with a monetary contribution, and a spend variable is one of them.
Warning
Critical Assessment: do not read this split as a recommendation to stop lower-funnel spend. It is a corner solution of a setup with two weaknesses.
First, no bounds keep the budget in the range of the data. The data never go below the always-on level of \(0.12\) times the geo size, so a budget of \(0\) is an extrapolation. Second, the objective treats the lower-funnel spend that the media induce as free, which favors media.
To use this feature for planning, pass budget_bounds that keep each budget in the range of the data, and price the induced spend separately, as the ROAS section does.
Requirements for Funnel Models#
Budget optimization with a custom effect has three requirements:
Channel data: the effect must read the channel data of the model.
create_effectreadsmmm.channel_data_scaled, so the intervention of the optimizer onchannel_datareaches the mediated path. An effect that read only its ownpm.Datavariables would not move, and the optimizer would price the direct path alone without any warning.Objective: the objective must be
total_response_original_scale, because the defaultresponse_variablecounts the direct path only.budget_optimizer()warns when it detects a mediated effect and no objective is named. ABudgetOptimizerbuilt directly on a model, as in the in-sample flow above, cannot warn, so name the objective.Exogenous data: the exogenous data of the effect must cover the window.
create_optimization_modeltakes care of this: it uses the observed values on the training dates and zeros after them, with the carryover sized byeffective_carryover_lags().
In-sample, this reproduces the training data, as above. For a future window, the zeros are a scenario choice (“no committed activity”). To plan against committed activity, set the variables of the effect with pm.set_data to the values that you expect before you pass the model to BudgetOptimizer.
A window that starts after the beginning of the training data also needs carry_in_periods. Then the carry-in, decision and adstock periods together cover the date range of the model. budget_optimizer() sets all three itself.
The optimizer does not count the lower-funnel spend that a plan induces. \(M^{*} = D + \lambda B\) rises with the demand that a plan creates, so the media-driven uplift above is gross of that spend. The optimizer only needs the total media-driven response, so the attribution problem of the shared saturating pool does not affect it.
Conclusion#
In this notebook, we compared six ways to treat a lower-funnel mediator in a geo-level MMM, on synthetic data with a known causal structure. We also optimized the budget of the fitted funnel model with two objectives.
Key Findings#
Mediated share: a bit over a third of the true effect of each upper-funnel channel (\(33.7\%\) to \(37.7\%\)) travels through the mediator. A model that cannot see the funnel has to attribute this mass to some channel.
Fit does not show causal bias: all six models fit the target almost equally well (in-sample \(R^2\) from \(0.9790\) to \(0.9871\)), but their ROAS estimates differ by tens of percent.
Funnel model accuracy: the funnel model covers the true total ROAS in all six geo-channel cells, with the narrowest mean interval width (\(0.422\)). Its mean absolute error is \(8.3\%\). Its largest miss is
south, social_spend(\(+26.8\%\)).Naive failure modes:
Naive Aanswers the direct-effect question, andNaive Bis confounded and misplaces the mediated demand.Naive B+andNaive B++reduce thetv_spendbias on average (Naive B+only partly, and not inwest), but neither fixessocial_spend.Naive Cblocks part of the effect that it tries to measure.Graph analysis: the d-separation analysis predicts most of these failures before any fit. However,
Naive B++correctstv_spendmore than the validNaive B+, andNaive B+has a larger mean ROAS error thanNaive AandNaive C, mostly onsocial_spend. A valid adjustment set says which variables you need, not how they enter the model or which structure is missing.Optimizer objective: the budget optimizer prices the funnel only if you ask it to. At the historical plan, the default objective sees \(67.1\%\) of the media-driven response under the posterior (\(60.8\%\) in truth), and \(56\%\) to \(72\%\) of the marginal response per cell.
Funnel-aware plan: the funnel-aware objective agrees with the default on the direction of five of six moves, and it reverses the move of
south, tv_spend.southcreates more demand per unit oftv_spendthannorth, and its lower funnel converts this demand at a steeper point of the shared curve.Uplift: on the truth, the funnel-aware plan gains \(+0.54\%\) in media-driven response against \(+0.28\%\) for the default plan, gross of the induced lower-funnel spend. Both gains are small, so the historical plan was already close to the posterior optimum within the \(\pm 50\%\) bounds.
Model Limitations#
True functional form: in this synthetic exercise, the funnel model is the DGP, with the right response families, pooling, anchoring and exclusion restrictions. On real data, nobody gives you the functional form, and misspecified structural equations would reduce the advantage of the funnel model.
Structure and information: the funnel model observes two likelihood series and two inputs that the naive models never see. Hence, its narrower intervals do not come from structure alone.
One dataset: every comparison uses one simulated dataset. Treat the coverage counts and the close rankings as a sanity check, not as a calibration study.
Graph assumptions: the verdicts depend on the graph, including the parts that are easy to leave out. Without the noise strand \(N\), the graph would accept
Naive B++as valid, and without the trend strand, it would acceptNaive B. The graph is also a single time slice, so a confounder that acts only with a lag would need a graph with explicit time.Latent confounders: on real data, the market driver can have strands that no regressor spans, for example taste shifts or competitor activity. Then even a trend control cannot close the backdoor, and a measured demand index helps only as far as the problem is confounding.
No organic path: we route the whole mediated effect through lower-funnel spend. A model with an organic path \(D \to Y\) next to \(M^{*} \to Y\) would be more realistic, but the two paths differ only by \(\lambda B\), so this sample would separate them only weakly.
Framework support: no built-in decomposition knows the funnel structure, because it lives in the code of
BudgetDualIndicatorEffect. The data-layer decomposition behind the plot suite andmmm.summary.contributionsomits a customMuEffectunder the identity link (underlink="log", it folds the effect into the baseline). The counterfactual decomposition shows the effect as onebudget_dual_effectbar, which is not the causal share of the channels.Pooling: we pool carryover and curvature across geos and leave the amplitudes free. Freeing every parameter per geo with about \(130\) weeks each would widen the posteriors considerably and can cause divergences.
Unscaled effect inputs: the base
MMMscales its channel and target data, but aDataVarMuEffectreads its variables unscaled. Here the data are close to unit scale, so theHalfNormal(1)priors of the demand equation fit the data.In-sample optimum: the optimization uses the training window with a truncated carryover at its end, so it finds an in-sample optimum, not a steady-state one.
Recommendations#
Model the mediator explicitly: when lower-funnel spend responds to upper-funnel demand, encode the mediation in a structural equation. Do not list the mediator as a channel. Drop it only if you need the total effect alone and you adjust for the backdoor.
No indicators as channels: ask “if we doubled it, would sales move?” Branded search fails this test, so it must not enter the model as a channel.
Anchor the latent scale: fix one coefficient on the latent demand (here the coefficient on \(D\) in the \(M^{*}\) equation), because otherwise the model is not identified. This choice decides what one unit of demand means.
Independent budget variation: the budget restriction identifies \(\lambda\) only if the budget varies on its own. Check its testable implication, \(B\) independent of \(S\) per geo on raw data.
Check the structural equations: use parameter recovery on simulated data and posterior predictive checks of the proxies to check that the structural equations are close enough to the data.
Use the funnel-aware objective: read the channel data in the effect and set
response_variable="total_response_original_scale"when you optimize a funnel model. To optimize the lower-funnel budget too, name it inspend_varswith bounds in the range of the data (see Lower-Funnel Co-Optimization).Build from a
Dataset: build the model from anxarray.Dataset, because aDataFramedrops the funnel columns, and fit it with the long frame. Scale the inputs of the effect (or its priors) yourself.State the attribution convention: with a shared saturating mediator, leave-one-out indirect effects understate the joint effect (here by \(13.3\%\)). Report the convention next to the number, as Mediation Analysis and (In)Direct Effects with PyMC explains with explicit interaction and dependence terms.
Price the induced spend: report the ROAS per incremental dollar, including the induced lower-funnel spend, next to the ROAS per upper-funnel dollar.
Next Steps#
Hierarchical priors: pool the geo parameters partially, instead of choosing between full pooling and no pooling.
Censored lower-funnel spend: model budget caps of lower-funnel channels with a censored likelihood.
Several demand pools: model more than one latent pool, for example branded and non-branded demand.
Time-varying baselines: let the baseline of the mediator change over time.
Information without structure: fit
Naive B++withlower_spendandsearch_volumeas control columns. This model sees the extra series without the structural system, but both columns are descendants of the treatment, so it answers a different causal question.Organic demand path: add a direct path \(D \to Y\), with data or priors that can separate it from the mediated path.
References#
The Nürnberger Versicherung case-study series from PyMC Labs works through the same funnel problem in a production setting:
Cracking the Funnel: How Nürnberger Quantified Cross-Channel Synergies (Part I)
Building a Funnel-Aware Bayesian Marketing Mix Model (Part II)
A PyMC-Marketing Case Study: Campaign Dynamics, Influencer Effectiveness, and Time-Slice CV (Part III)
On the graphical methods of the D-Separation Analysis:
Pearl, J., Glymour, M., and Jewell, N. P. (2016). Causal Inference in Statistics: A Primer. Wiley. Chapters 2 and 3 cover d-separation and the backdoor criterion.
Pearl, J., and Mackenzie, D. (2018). The Book of Why: The New Science of Cause and Effect. Basic Books.
%load_ext watermark
%watermark -n -u -v -iv -w -p pymc_marketing,pytensor
Last updated: Tue, 13 Oct 2026
Python implementation: CPython
Python version : 3.12.12
IPython version : 9.15.0
pymc_marketing: 1.2.0
pytensor : 3.3.2
arviz : 1.2.0
arviz_plots : 1.2.0
graphviz : 0.21
matplotlib : 3.11.1
networkx : 3.6.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.2.0
pytensor : 3.3.2
xarray : 2026.4.0
Watermark: 2.6.0