Calibrating an MMM with ROAS Estimates#

Marketing mix models are fit to sales data, and sales data alone is often not enough to pin down how much each channel contributed. Two very different sets of channel parameters can reproduce the same sales curve almost equally well. When that happens, the model is not wrong so much as under-determined, and the posterior will happily settle on values that a marketing team knows to be implausible.

Experiments are the usual way out. A geo lift test, a conversion lift study, or a switchback experiment gives a direct, causal estimate of what a channel returns. The question this notebook answers is mechanical: how do you get that number into the model?

PyMC-Marketing exposes this as a one-line method, add_cost_per_target_calibration(), which attaches an observed likelihood to the model’s implied return on ad spend (ROAS) or cost per acquisition (CPA).

What you will learn

  1. What likelihood calibration means and why it is not the same as putting a prior on a parameter.

  2. How to express an experiment result as a ROAS calibration table.

  3. How to add it to an MMM with add_cost_per_target_calibration.

  4. How the sigma column controls how strongly the model is pulled.

  5. How to compare a naive model against calibrated ones.

  6. How the likelihood family (Normal, StudentT or LogNormal) sets what a large disagreement costs, and when to switch from the default Normal.

Prerequisites

This notebook assumes you have seen the Media Mix Modeling (MMM) in Python: End-to-End Example notebook and are comfortable building an MMM. No prior exposure to calibration is assumed.

We deliberately reuse the simulated dataset from the Mitigating Unobserved Confounders in MMMs with Lift Test Likelihoods notebook. That dataset contains a known unobserved confounder, so we know the true ROAS and can check whether calibration actually recovers it. Where that notebook calibrates through lift tests (points on the saturation curve), this one calibrates directly on the ROAS quantity itself.

Prepare Notebook#

import warnings

import arviz as az
import arviz_plots as azp
import graphviz as gr
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import seaborn as sns
import xarray as xr
from matplotlib.colors import to_hex
from matplotlib.lines import Line2D
from pymc_extras.prior import Prior
from xarray import DataArray

from pymc_marketing.hsgp_kwargs import HSGPKwargs
from pymc_marketing.metrics import crps
from pymc_marketing.mmm import MMM, GeometricAdstock, LogisticSaturation
from pymc_marketing.paths import data_dir

warnings.filterwarnings("ignore", category=FutureWarning)

az.style.use("arviz-darkgrid")
plt.rcParams["figure.figsize"] = [12, 7]
plt.rcParams["figure.dpi"] = 100
plt.rcParams["figure.facecolor"] = "white"

%load_ext autoreload
%autoreload 2
%config InlineBackend.figure_format = "retina"

seed: int = sum(map(ord, "mmm_roas_calibration"))
rng: np.random.Generator = np.random.default_rng(seed=seed)

What Is Likelihood Calibration?#

Before touching any code it is worth being precise about what we are doing, because the phrase “ROAS prior” is used loosely in the industry and can mean two quite different things.

A Bayesian model has two ingredients: priors over parameters and a likelihood connecting parameters to observed data. In an MMM the parameters are things like the adstock decay alpha, the saturation lam, and the channel coefficient beta. The likelihood says how sales y are generated from them.

ROAS is not one of those parameters. It is a derived quantity: you have to run the adstock transform, the saturation transform, sum the channel contribution over time, and divide by spend before you get a number you can call ROAS. In PyMC terms it is a deterministic function of the parameters, and you cannot put a prior on a deterministic. A prior is a statement about a free variable; once a quantity is a function of other variables its distribution is already implied.

So what do you do when you have external information about a derived quantity? You treat that information as data. Concretely, you add a second observed node to the model:

\[ \text{ROAS}_{\text{experiment}} \sim \mathrm{Normal}\big(\underbrace{f(\theta)}_{\text{model-implied ROAS}},\ \sigma\big) \]

Read this literally. The experiment produced a number. The model, for any given draw of its parameters \(\theta\), implies its own number. We declare that the experiment’s number is a noisy observation of the model’s number, with noise scale \(\sigma\). Parameter draws whose implied ROAS sits far from the experiment now score badly in the log-likelihood, and the sampler stops visiting them. The Normal is the default shape of that noise; a later section replaces it.

That is likelihood calibration: encoding external knowledge as an extra observation rather than as a prior.

Why this is a soft constraint, not a hard one

Nothing here forces the model to reproduce the experimental ROAS exactly. The sales likelihood is still in play and still pulling in its own direction. The posterior is a compromise between the two, and sigma sets the exchange rate.

  • Small sigma → the calibration term dominates. The model is pulled hard onto the experimental value. Use this when you trust the experiment far more than the observational sales data.

  • Large sigma → a gentle nudge. Use this when the experiment was noisy, ran on a different period, or measured a slightly different thing.

sigma is your statement about how much you trust the experiment. It is not a tuning knob to be optimized. Set it from the experiment’s own standard error whenever you have one.

Relationship to lift test calibration

The Lift Test Calibration notebook calibrates using the saturation curve: a lift test gives you two points (pre-spend and post-spend) and the inferred change in sales, which constrains the shape of the curve locally. This approach allows you to pass more information to the model as you are calibrating the model with marginal ROAS at different spend levels.

Tip

Three notebooks in the gallery get ROAS information into an MMM, and they do it in three different ways. It is worth knowing which one you want before you start.

  • Custom MMM with ROAS Parameterization reparametrizes the model. Rather than sampling a channel coefficient and deriving ROAS from it, the model is written so that ROAS is itself a free parameter. Then you really can place a prior on it, because it is no longer a deterministic. This is the most direct route, at the cost of writing the model yourself.

  • Mitigating Unobserved Confounders in MMMs with Lift Test Likelihoods brings experiments in as lift tests. That is also an extra likelihood term, exactly like the one in this notebook, but it carries more information: a lift test says what happened when spend moved from one level to another, so it constrains points on the saturation curve rather than a single whole-period average.

  • This notebook covers the simplest case: one ROAS or CPA number per channel, aggregated over the whole modeling period, added as a single observation.

The three are not competing answers to one question. If you have a full experiment readout, the lift-test route uses more of it. If all that reached you is a headline ROAS figure, this is the method that fits what you have.

The Data#

We use the simulated dataset shipped with PyMC-Marketing and analyzed in the Mitigating Unobserved Confounders in MMMs with Lift Test Likelihoods notebook. It comes from the blog post Media Mix Model and Experimental Calibration: A Simulation Study.

data_path = data_dir / "mmm_roas_data.csv"

raw_df = pd.read_csv(data_path, parse_dates=["date"])

model_df = raw_df.copy().filter(["date", "x1", "x2", "y"])

model_df.head()
date x1 x2 y
0 2021-10-02 0.646554 0.336188 199.329637
1 2021-10-09 1.411917 0.203931 371.237041
2 2021-10-16 0.837610 0.024026 272.215933
3 2021-10-23 0.973612 0.120257 291.104040
4 2021-10-30 1.415985 0.084630 386.243000

For modeling we only use date, the two channels x1 and x2, and the target y. The remaining columns in raw_df describe the data generating process and let us check our answers later.

The generating structure is the following.

g = gr.Digraph()
g.node(name="seasonality", label="seasonality", color="lightgray", style="filled")
g.node(name="trend", label="trend")
g.node(name="z", label="z", color="lightgray", style="filled")
g.node(name="x1", label="x1", color="#2a2eec80", style="filled")
g.node(name="x2", label="x2", color="#fa7c1780", style="filled")
g.node(name="y", label="y", color="#328c0680", style="filled")
g.edge(tail_name="seasonality", head_name="x1")
g.edge(tail_name="z", head_name="x1")
g.edge(tail_name="x1", head_name="y")
g.edge(tail_name="seasonality", head_name="y")
g.edge(tail_name="trend", head_name="y")
g.edge(tail_name="z", head_name="y")
g.edge(tail_name="x2", head_name="y")
g
../../_images/7b90a46e3b5413ba00ba635faf22c8e4946084ec622dfe73eff7969af114ead7.svg

The key feature is z: a variable that drives both spend on x1 and sales y, and which we do not observe. This is a textbook confounder. Because z is missing from the model, part of its effect on y gets attributed to x1, and the model will overstate how effective x1 is.

Channel x2 carries no confounding bias of its own: nothing in the graph drives both x2 and y. That does not mean the model recovers it cleanly. When x1 is overstated, the rest of the decomposition has to give up sales it was explaining, and x2 is part of that remainder, so its estimate can be dragged off target indirectly. The contrast between the two channels is therefore direct bias versus inherited bias, not bias versus a clean control.

Since we simulated the data, we know the true ROAS of each channel.

true_roas_x1 = (raw_df["y"] - raw_df["y01"]).sum() / raw_df["x1"].sum()
true_roas_x2 = (raw_df["y"] - raw_df["y02"]).sum() / raw_df["x2"].sum()

true_roas = {"x1": true_roas_x1, "x2": true_roas_x2}

print(f"True ROAS for x1: {true_roas_x1:.2f}")
print(f"True ROAS for x2: {true_roas_x2:.2f}")

True ROAS for x1: 93.39
True ROAS for x2: 171.41

Note

These ROAS values are in the hundreds because the simulated target is scaled by an amplitude factor. Do not read them as realistic revenue multiples. What matters for this notebook is the ratio between the two channels and whether we recover it.

Let us look at the series before modeling.

fig, ax = plt.subplots(
    nrows=2,
    ncols=1,
    sharex=True,
    sharey=False,
    layout="constrained",
)

sns.lineplot(x="date", y="y", data=model_df, color="black", ax=ax[0])
ax[0].set_title("Target Data")

model_df.melt(
    id_vars=["date"], value_vars=["x1", "x2"], var_name="channel", value_name="x"
).pipe(
    (sns.lineplot, "data"),
    x="date",
    y="x",
    hue="channel",
    ax=ax[1],
)
ax[1].legend(title="Channel", title_fontsize=12)
ax[1].set_title("Channels Spend Data");
../../_images/598ecfefef38687cd8177cf89013970c22d852654f8db142cbf52c4f1856c23c.png

The Naive Model#

We start with a model that sees only date, x1, x2 and y. No experiment information at all. This is our baseline and it uses exactly the same configuration as the baseline model in the Mitigating Unobserved Confounders in MMMs with Lift Test Likelihoods notebook, so the two notebooks are directly comparable.

Without external information, a sensible default is to encode the cost share of each channel into the prior of its coefficient: channels that absorb more of the budget are given more room to contribute.

cost_share = DataArray(
    model_df[["x1", "x2"]].sum() / model_df[["x1", "x2"]].sum().sum(),
    dims="channel",
)

cost_share

<xarray.DataArray (channel: 2)> Size: 16B
array([0.7366047, 0.2633953])
Coordinates:
  * channel  (channel) object 16B 'x1' 'x2'
model_config = {
    "likelihood": Prior("Normal", sigma=Prior("HalfNormal", sigma=2)),
    "gamma_fourier": Prior("Normal", mu=0, sigma=2, dims="fourier_mode"),
    "intercept_tvp_config": HSGPKwargs(
        m=100, L=None, eta_lam=1.0, ls_mu=5.0, ls_sigma=10.0, cov_func=None
    ),
    "adstock_alpha": Prior("Beta", alpha=2, beta=3, dims="channel"),
    "saturation_lam": Prior("Gamma", alpha=2, beta=2, dims="channel"),
    "saturation_beta": Prior("HalfNormal", sigma=cost_share, dims="channel"),
}

A Gaussian process handles the non-linear trend through a time-varying intercept (see MMM with time-varying parameters (TVP)). We collect the constructor arguments and the sampler settings in dictionaries so that every model in this notebook is built identically and only the calibration differs.

mmm_kwargs = dict(
    adstock=GeometricAdstock(l_max=4),
    saturation=LogisticSaturation(),
    date_column="date",
    channel_columns=["x1", "x2"],
    target_column="y",
    time_varying_intercept=True,
    time_varying_media=False,
    yearly_seasonality=5,
    model_config=model_config,
)

sampler_config = {
    "tune": 1_000,
    "chains": 4,
    "draws": 1_000,
    "nuts_sampler": "nutpie",
    "target_accept": 0.95,
    "random_seed": rng,
}

y = model_df["y"]
X = model_df.drop(columns=["y"])

One detail matters for everything that follows. Calibration compares model contribution against spend in original units, so the model has to expose the channel contribution on its original scale, channel_contribution_original_scale. add_original_scale_contribution_variable registers it, and add_cost_per_target_calibration registers it too when it is missing. We call the former explicitly because the plots below also read the other components on the original scale.

%%time

mmm_naive = MMM(**mmm_kwargs)

vars_of_interest = [
    "channel_contribution",
    "fourier_contribution",
    "intercept_contribution",
    "y",
]

mmm_naive.build_model(X, y)
mmm_naive.add_original_scale_contribution_variable(var=vars_of_interest)
_ = mmm_naive.fit(X, y, **sampler_config)

NUTS[nutpie]: [y_sigma, gamma_fourier, adstock_alpha, saturation_lam, saturation_beta, intercept_latent_process_raw_hsgp_coefs_offset, intercept_latent_process_raw_eta, intercept_latent_process_raw_ls, intercept_baseline]





CPU times: user 27.5 s, sys: 421 ms, total: 28 s
Wall time: 13 s













# Number of diverging samples
mmm_naive.idata["sample_stats"]["diverging"].sum().item()








../../_images/ab94bd3b9d42a4fe80b1028678599c507867a8521a469af5018c945a3dba9e72.png





There are no divergences, but divergences do not reveal a chain stuck in a minor mode. A quick check is the mean joint log density of each chain: the log prior of every free parameter plus the log-likelihood of y.






# Joint log density of each chain: log prior of every free parameter plus log-likelihood of y
log_lik_y = pm.stats.compute_log_likelihood(
    mmm_naive.idata,
    model=mmm_naive.model,
    extend_inferencedata=False,
    progressbar=False,
)["y"].sum(dim="date")
log_prior = pm.stats.compute_log_prior(
    mmm_naive.idata,
    model=mmm_naive.model,
    extend_inferencedata=False,
    progressbar=False,
)
log_prior_total = sum(
    lp.sum(dim=[dim for dim in lp.dims if dim not in ("chain", "draw")])
    for lp in log_prior.data_vars.values()
)
chain_logp = (log_prior_total + log_lik_y).mean("draw").rename("logp")

# Keep only the chains within 10 nats of the best one
good_chains = np.flatnonzero((chain_logp > chain_logp.max() - 10).to_numpy())
mmm_naive.idata = mmm_naive.idata.isel(chain=good_chains, missing_dims="ignore")

chain_logp.round(1).to_series()









chain
0    123.7
1    124.2
2    123.8
3     80.2
Name: logp, dtype: float64











Dropping a stuck chain


One chain sits about \(44\) nats below the other three. It is stuck in a minor mode where x1 saturates heavily and its ROAS lands near \(464\), and it fits y worse. A mode that far down holds a negligible share of the posterior, yet pooling all four chains would give it a quarter of the draws, so we keep only the chains within \(10\) nats of the best one. This is the pruning of minor modes recommended in Bayesian Workflow (Gelman et al., 2020, Section 5.6), and what stacking the chains (Yao, Vehtari and Gelman, 2022) does approximately. Only drop a chain once you know why it differs.




Now let us look at what the model believes about ROAS. The incrementality accessor computes contribution over spend for us.


One bookkeeping note before reading any plots. Three closely related quantities
appear in this notebook: the plotted ROAS below comes from
contribution_over_spend, a zero-spend counterfactual that includes carryover
effects; the calibration target added later constrains the ratio of mean
channel contribution to mean spend over the fitted dates; and the true ROAS
was computed above directly from the simulation’s ground-truth columns. For an
additive model like this one they agree in principle, but they are three
different computations, so small numerical differences between them are expected
and not a bug.






def get_roas(mmm: MMM) -> DataArray:
    """Posterior ROAS per channel, aggregated over the whole period."""
    return mmm.incrementality.contribution_over_spend(frequency="all_time").rename(
        "roas"
    )


roas_naive = get_roas(mmm_naive)













pc = azp.plot_dist(
    roas_naive.to_dataset(name="roas"),
    col_wrap=1,
    figure_kwargs={
        "figsize": (12, 7),
        "sharex": True,
        "sharey": False,
        "layout": "constrained",
    },
)
fig = pc.viz["/"]["figure"].values.item()
axes = fig.axes
for ax, channel in zip(axes, ["x1", "x2"], strict=False):
    ax.axvline(
        true_roas[channel],
        color="black",
        linestyle="--",
        linewidth=3,
        label="True ROAS",
    )
    ax.set(title=f"Channel ${channel[0]}_{{{channel[1]}}}$")
axes[-1].set(xlabel="ROAS")
handles, labels = axes[0].get_legend_handles_labels()
fig.legend(handles, labels, loc="outside right upper")
fig.suptitle("ROAS Posterior - Naive Model", fontsize=16, fontweight="bold");








../../_images/49c1aa916892d2fb758b790a443b5f47864ac0c327049db5a6588b1d2902c81a.png






Neither channel is recovered, and x1 is the large failure. The posterior mean puts it at about \(197\) against a true ROAS of \(93\), more than twice the truth. This is the confounder at work.


x2 is the milder case: about \(142\), against a true \(171\), so low by roughly \(17\%\). Note that x2 is not confounded, so its error is inherited rather than direct. When one channel’s contribution is overstated, the rest of the model has to give up some of the sales it was explaining, and x2 is part of that remainder. This notebook does not isolate that mechanism, so treat it as the likely reading rather than a demonstrated one.


The ordering is the consequence that matters most. The model ranks x1 above x2, and the truth is the opposite. A budget recommendation built on this posterior moves money in the wrong direction.


No amount of extra sales data fixes this.
The bias is structural, not statistical.
We need information from outside the sales series.






Adding a ROAS Calibration#


Suppose the measurement team ran geo lift experiments on both channels and
reported a ROAS point estimate with a standard error for each. That is exactly
the input add_cost_per_target_calibration expects.


The calibration table needs one row per experiment, with:





column


meaning






channel


which channel the experiment measured, matching channel_columns




value column


the measured ROAS (or CPA); you choose the column name




sigma


how much deviation you accept, on the same scale as the value




dims columns


one column per model dimension, for example geo, when the model has them









Our model has no extra dimensions, so channel, the value, and sigma are all
we need.






roas_calibration_df = pd.DataFrame(
    {
        "channel": ["x1", "x2"],
        "roas": [true_roas_x1, true_roas_x2],
        "sigma": [0.15 * true_roas_x1, 0.15 * true_roas_x2],
    }
)

roas_calibration_df













  
    

      
      channel
      roas
      sigma
    

  
  
    

      0
      x1
      93.388131
      14.008220
    

    

      1
      x2
      171.408804
      25.711321
    

  









Where these numbers come from in practice


We plug in the true ROAS because this is a simulation and we want to check
recovery. In a real project the value column comes from an experiment readout
and sigma comes from that experiment’s own uncertainty, not from a guess. A
geo lift analysis in CausalPy, for
instance, gives you a posterior over the incremental effect; its standard
deviation is a defensible sigma.


Never calibrate to a number that came out of another MMM. You would be feeding
the model its own opinion back as evidence.




Here sigma is \(15\%\) of each ROAS value, which is a moderate stance: we take the
experiments seriously but leave the sales likelihood some room to disagree.


Now we build the calibrated model. Note the ordering:


    

  1. build_model: the calibration needs a model to attach to.
    2. 

  2. add_original_scale_contribution_variable: optional for the calibration, which registers channel_contribution_original_scale itself when it is missing. We call it so the calibrated model exposes the same original-scale variables as the naive one.
    3. 

  3. add_cost_per_target_calibration: add the observed likelihood.
    4. 

  4. fit.
    5. 







mmm_calibrated = MMM(**mmm_kwargs)

mmm_calibrated.build_model(X, y)
mmm_calibrated.add_original_scale_contribution_variable(var=vars_of_interest)

mmm_calibrated.add_cost_per_target_calibration(
    data=X,
    calibration_data=roas_calibration_df,
    name_prefix="roas_calibration",
    target_column="roas",
    target_per_cost=True,
);









Three arguments deserve comment.


data is a spend table with the same layout as X. In this dataset the channel
columns already are spend, so we pass X directly. If your model is fit on
impressions or GRPs rather than money, this is where you supply the matching
monetary spend, and getting it right is what makes the resulting ratio a real
ROAS.


target_column tells the method which column holds the calibration values. It
defaults to "cost_per_target"; we named ours "roas", so we say so.


target_per_cost=True is the direction switch, and it is the important one:


    

  • target_per_cost=False (the default) calibrates
\(\overline{\text{spend}} \,/\, \overline{\text{contribution}}\), which is
cost per target: CPA, CPI, cost per conversion.
    • 

  • target_per_cost=True calibrates
\(\overline{\text{contribution}} \,/\, \overline{\text{spend}}\), which is
ROAS when contribution is revenue.
    • 



The two are reciprocals, so in principle either could express the same
constraint. In practice they are not interchangeable, because sigma is applied
on whichever scale you pick. A symmetric Normal error on ROAS is not a symmetric
error on CPA. For a lift test the ROAS scale is the one whose error is Normal, so
a CPA from one is converted to a ROAS table before it is calibrated (see the
warning in Calibrating on CPA Instead).




The likelihood family is a choice


The calibration term is a Normal by default. likelihood=Prior("StudentT", nu=...)
and likelihood=Prior("LogNormal") are the alternatives. They differ in what
happens when the model and the experiment disagree by a lot, not near the
experiment, and the direction switch above matters a great deal for the
first two and not at all for the third. The section
Choosing the Likelihood Family shows them side by side.




Let us confirm what actually got added to the model.






observed_rvs = [rv.name for rv in mmm_calibrated.model.observed_RVs]

print(
    f"Naive model observed variables:      {[rv.name for rv in mmm_naive.model.observed_RVs]}"
)
print(f"Calibrated model observed variables: {observed_rvs}")









Naive model observed variables:      ['y']
Calibrated model observed variables: ['y', 'roas_calibration']









The naive model observes only y. The calibrated model observes y and
roas_calibration. That second node is the whole mechanism: two likelihood terms
instead of one, and the posterior has to satisfy both.


Now we fit.






%%time

_ = mmm_calibrated.fit(X, y, **sampler_config)









NUTS[nutpie]: [y_sigma, gamma_fourier, adstock_alpha, saturation_lam, saturation_beta, intercept_latent_process_raw_hsgp_coefs_offset, intercept_latent_process_raw_eta, intercept_latent_process_raw_ls, intercept_baseline]









CPU times: user 25.2 s, sys: 316 ms, total: 25.5 s
Wall time: 10.6 s













# Number of diverging samples
mmm_calibrated.idata["sample_stats"]["diverging"].sum().item()








../../_images/ab94bd3b9d42a4fe80b1028678599c507867a8521a469af5018c945a3dba9e72.png









roas_calibrated = get_roas(mmm_calibrated)

pc = azp.plot_dist(
    roas_calibrated.to_dataset(name="roas"),
    col_wrap=1,
    figure_kwargs={
        "figsize": (12, 7),
        "sharex": True,
        "sharey": False,
        "layout": "constrained",
    },
)
fig = pc.viz["/"]["figure"].values.item()
axes = fig.axes
for ax, channel in zip(axes, ["x1", "x2"], strict=False):
    ax.axvline(
        true_roas[channel],
        color="black",
        linestyle="--",
        linewidth=3,
        label="True ROAS",
    )
    ax.legend(loc="upper right")
    ax.set(title=f"Channel ${channel[0]}_{{{channel[1]}}}$", xlim=(80, 300))
axes[-1].set(xlabel="ROAS")
fig.suptitle("ROAS Posterior - Calibrated Model", fontsize=16, fontweight="bold");








../../_images/9612193bafd01a6975dfd0ad8c12988febca08972dead6f5fd0134c5155a232d.png






This calibration helps the model estimates, but it only partially mitigates the confounding effect.
The x1 posterior mean moves from the naive model’s \(197\) down to \(155\) against a true \(93\), so roughly \(40\%\) of the error is gone, and x2 moves from \(142\) up to \(149\) against a true \(171\).
The true x1 ROAS is still outside its \(94\%\) interval, and the ranking is still inverted: about \(155\) for x1 against \(149\) for x2.


That is a weaker result than “calibration fixes it”, and the reason is sigma.
We asked for \(15\%\) relative uncertainty, which leaves the sales likelihood a lot of room to disagree, and it uses it.
The next section refits with a tighter and a looser calibration so you can see how much of the answer is controlled by that one number.






How sigma Controls the Pull#


sigma is the only real dial and it is worth seeing it move. We refit the model
with a tight calibration (\(5\%\) relative uncertainty) and a loose one (\(60\%\)), keeping
everything else fixed.






%%time

sigma_scenarios = {"tight (5%)": 0.05, "loose (60%)": 0.60}
mmm_sigma = {}

for label, rel_sigma in sigma_scenarios.items():
    calibration_df = roas_calibration_df.assign(
        sigma=lambda df, r=rel_sigma: r * df["roas"]
    )

    mmm_s = MMM(**mmm_kwargs)
    mmm_s.build_model(X, y)
    mmm_s.add_original_scale_contribution_variable(var=vars_of_interest)
    mmm_s.add_cost_per_target_calibration(
        data=X,
        calibration_data=calibration_df,
        name_prefix="roas_calibration",
        target_column="roas",
        target_per_cost=True,
    )
    _ = mmm_s.fit(X, y, **sampler_config)
    mmm_sigma[label] = mmm_s









NUTS[nutpie]: [y_sigma, gamma_fourier, adstock_alpha, saturation_lam, saturation_beta, intercept_latent_process_raw_hsgp_coefs_offset, intercept_latent_process_raw_eta, intercept_latent_process_raw_ls, intercept_baseline]









NUTS[nutpie]: [y_sigma, gamma_fourier, adstock_alpha, saturation_lam, saturation_beta, intercept_latent_process_raw_hsgp_coefs_offset, intercept_latent_process_raw_eta, intercept_latent_process_raw_ls, intercept_baseline]









CPU times: user 51.6 s, sys: 656 ms, total: 52.3 s
Wall time: 21.7 s













Comparing the Models#


A forest plot is the natural way to compare posteriors across models: it puts
every interval on one axis so the shrinkage is immediately visible.






roas_by_model = {
    "naive": roas_naive,
    "calibrated, tight (5%)": get_roas(mmm_sigma["tight (5%)"]),
    "calibrated, moderate (15%)": roas_calibrated,
    "calibrated, loose (60%)": get_roas(mmm_sigma["loose (60%)"]),
}

# The naive model kept three chains and the others four, so concat pads the
# naive fourth chain with NaN; the summaries and plots below skip it.
combined = xr.concat(
    [roas.to_dataset() for roas in roas_by_model.values()], dim="model"
).assign_coords(model=list(roas_by_model))

dt = xr.DataTree.from_dict({"/posterior": xr.DataTree(combined)})

combined

























<xarray.Dataset> Size: 264kB
Dimensions:  (model: 4, chain: 4, draw: 1000, channel: 2)
Coordinates:
  * model    (model) <U26 416B 'naive' ... 'calibrated, loose (60%)'
  * chain    (chain) int64 32B 0 1 2 3
  * draw     (draw) int64 8kB 0 1 2 3 4 5 6 7 ... 993 994 995 996 997 998 999
  * channel  (channel) <U2 16B 'x1' 'x2'
Data variables:
    roas     (model, chain, draw, channel) float64 256kB 221.9 129.2 ... 141.2








def plot_model_forest(
    datatree: xr.DataTree,
    channel: str,
    true_value: float,
    colors: list[str] | None = None,
):
    """Forest plot of ROAS posteriors across models for a single channel.

    ``colors`` sets one color per model, in the order of the ``model``
    coordinate. The default uses the color cycle.
    """
    color_kwargs = {} if colors is None else {"color": colors}
    pc = azp.plot_forest(
        datatree.sel(channel=channel),
        combined=True,
        figure_kwargs={"figsize": (10, 5)},
        # Skip the NaN padding of models with fewer chains
        stats={k: {"skipna": True} for k in ("trunk", "twig", "point_estimate")},
        **color_kwargs,
    )
    ax = pc.viz["plot"].sel(column="forest").item()
    ax.axvline(true_value, color="black", linestyle="--", linewidth=2)
    ax.set(xlabel="ROAS")

    color_map = pc.aes["color"]["mapping"]
    handles = [
        Line2D([0], [0], color=str(color_map.sel(model=m).values), linewidth=3)
        for m in datatree["posterior"].coords["model"].values
    ]
    handles.append(Line2D([0], [0], color="black", linestyle="--", linewidth=2))
    labels = [*list(datatree["posterior"].coords["model"].values), "True ROAS"]
    ax.legend(handles, labels, loc="center left", bbox_to_anchor=(1.02, 0.5))
    return pc, ax


pc, ax = plot_model_forest(dt, channel="x1", true_value=true_roas_x1)
ax.set_title(
    "ROAS Posterior by Model - Channel $x_{1}$ (confounded)",
    fontsize=18,
    fontweight="bold",
    y=1.05,
);








../../_images/0feb38deb57b629996e45a6919eea973144bd8facb875766a35dd98513885c26.png










pc, ax = plot_model_forest(dt, channel="x2", true_value=true_roas_x2)
ax.set_title(
    "ROAS Posterior by Model - Channel $x_{2}$ (unconfounded)",
    fontsize=18,
    fontweight="bold",
    y=1.05,
);








../../_images/0938fc35c06553deee4b94511363b2e4029f41cc5e47edfaf1f5c7a5fbc1142b.png










# Posterior mean and 94% HDI of ROAS per model and channel, the numbers the
# text below quotes.
roas_hdi = az.hdi(combined["roas"], prob=0.94, skipna=True)
roas_summary = xr.Dataset(
    {
        "mean": combined["roas"].mean(dim=("chain", "draw")),
        "hdi_lower": roas_hdi.sel(ci_bound="lower", drop=True),
        "hdi_upper": roas_hdi.sel(ci_bound="upper", drop=True),
    }
).to_dataframe()

roas_summary.round(1)













  
    

      
      
      mean
      hdi_lower
      hdi_upper
    

    

      channel
      model
      
      
      
    

  
  
    

      x1
      naive
      197.0
      171.8
      220.3
    

    

      calibrated, tight (5%)
      107.7
      99.6
      116.6
    

    

      calibrated, moderate (15%)
      154.6
      138.3
      170.8
    

    

      calibrated, loose (60%)
      192.8
      168.8
      216.7
    

    

      x2
      naive
      142.2
      120.3
      164.8
    

    

      calibrated, tight (5%)
      164.7
      150.7
      179.0
    

    

      calibrated, moderate (15%)
      149.0
      127.0
      169.7
    

    

      calibrated, loose (60%)
      143.0
      120.6
      164.7
    

  











# Posterior probability that the ranking is correct (true ROAS: x2 > x1),
# computed per model from the joint posterior rather than from point estimates.
prob_correct_ranking = (
    (combined["roas"].sel(channel="x2") > combined["roas"].sel(channel="x1"))
    .where(combined["roas"].notnull().all("channel"))
    .mean(dim=("chain", "draw"))
    .rename("P(ROAS_x2 > ROAS_x1)")
    .to_dataframe()
)

prob_correct_ranking













  
    

      
      P(ROAS_x2 > ROAS_x1)
    

    

      model
      
    

  
  
    

      naive
      0.000667
    

    

      calibrated, tight (5%)
      1.000000
    

    

      calibrated, moderate (15%)
      0.342250
    

    

      calibrated, loose (60%)
      0.002000
    

  







The x1 panel, read against the dashed line at the true ROAS of \(93\), tells a more nuanced story than “calibration works”.


Two things work as advertised:


    

  • Every calibrated model moves toward the truth, and the ordering follows sigma exactly.
The tight (\(5\%\)) calibration lands nearest with a posterior mean near \(108\), the moderate (\(15\%\)) one at about \(155\), and the loose (\(60\%\)) one at about \(193\), against the naive model’s \(197\).
    • 

  • The loose calibration stays closest to the naive fit. It settles almost exactly on the naive fit, with an interval of about the same width.
That is the expected behavior when you have told the model not to trust the experiment.
    • 



One thing does not work.
Not even the tight calibration recovers the true value: its \(94\%\) interval runs from about \(100\) to \(116\), which excludes \(93\) at the low end.
The calibration enters as a single observation that competes against every date in the sales series.
So even at \(5\%\) relative uncertainty, the posterior mean of about \(108\) settles roughly three sigma (here \(5\%\) of \(93\), about \(4.7\)) above the calibration target of \(93\) and not on it.
Tightening sigma buys a large reduction in bias here, not exact recovery.
If you need the experiment to constrain the model more firmly than that, a whole-period ratio is the wrong instrument.
See the note on matching the time window below, and Lift Test Calibration for a constraint at the spend levels the experiment tested.


The x2 panel shows the inherited bias described earlier.
The naive model was not already right about x2: its mean sits near \(142\) against a true \(171\), low by about \(17\%\), with the truth outside its interval.
Calibration helps here too. The tight interval covers the true value, and the moderate one ends just short of it (its upper end is about \(170\)).


The summary that matters is the ranking, because that is what moves budget.
The truth is x2 above x1, and the table above gives the posterior probability of that ordering under each model.
This probability is the right way to read a ranking off two overlapping posteriors. The difference between the posterior means is small relative to their spread, so the means alone would overstate the precision.
The naive model and the loose calibration put that probability below \(0.1\%\) and at about \(0.2\%\), so they bet confidently on the wrong order.
The moderate calibration moves it to about \(35\%\), closer to a coin flip but still in favor of the wrong order.
Only the tight calibration commits to the correct order, at essentially \(100\%\).


Note what made the difference: the tight calibration describes a more precise experiment, not a better choice of sigma for the same experiment.
If the experiment’s standard error is \(15\%\), the moderate fit is the right one, and on this dataset it does not overturn the confounding.
Also, every experiment in this notebook reports the true ROAS, because we set it that way.
What happens to the tight calibration when the experiment itself is biased?




When the Experiment Is Biased#


Let’s repeat the tight (\(5\%\)) calibration with an x1 experiment that overstates the true ROAS, once by \(50\%\) and once by \(100\%\). The x2 experiment stays exact. We build the two tables first.






# Overstatement factor of the x1 experiment
x1_bias = {"x1 experiment +50%": 1.5, "x1 experiment +100%": 2.0}

# label -> calibration table with the overstated x1 value and a 5% sigma
biased_tables = {
    label: roas_calibration_df.assign(
        roas=[factor * true_roas_x1, true_roas_x2],
        sigma=lambda df: 0.05 * df["roas"],
    )
    for label, factor in x1_bias.items()
}

pd.concat(biased_tables, names=["scenario", None])













  
    

      
      
      channel
      roas
      sigma
    

    

      scenario
      
      
      
      
    

  
  
    

      x1 experiment +50%
      0
      x1
      140.082197
      7.004110
    

    

      1
      x2
      171.408804
      8.570440
    

    

      x1 experiment +100%
      0
      x1
      186.776262
      9.338813
    

    

      1
      x2
      171.408804
      8.570440
    

  







Now we fit one model per table.






%%time

mmm_biased = {}

for label, calibration_df in biased_tables.items():
    mmm_b = MMM(**mmm_kwargs)
    mmm_b.build_model(X, y)
    mmm_b.add_original_scale_contribution_variable(var=vars_of_interest)
    mmm_b.add_cost_per_target_calibration(
        data=X,
        calibration_data=calibration_df,
        name_prefix="roas_calibration",
        target_column="roas",
        target_per_cost=True,
    )
    _ = mmm_b.fit(X, y, **sampler_config)
    mmm_biased[label] = mmm_b









NUTS[nutpie]: [y_sigma, gamma_fourier, adstock_alpha, saturation_lam, saturation_beta, intercept_latent_process_raw_hsgp_coefs_offset, intercept_latent_process_raw_eta, intercept_latent_process_raw_ls, intercept_baseline]









NUTS[nutpie]: [y_sigma, gamma_fourier, adstock_alpha, saturation_lam, saturation_beta, intercept_latent_process_raw_hsgp_coefs_offset, intercept_latent_process_raw_eta, intercept_latent_process_raw_ls, intercept_baseline]









CPU times: user 52.6 s, sys: 825 ms, total: 53.4 s
Wall time: 23.3 s









We compare the three tight calibrations, on the exact experiment and on the two overstated ones, with one forest plot per channel.








Hide code cell source





roas_off = {
    "exact experiment": roas_by_model["calibrated, tight (5%)"],
    **{label: get_roas(m) for label, m in mmm_biased.items()},
}

combined_off = xr.concat(
    [roas.to_dataset() for roas in roas_off.values()], dim="model"
).assign_coords(model=list(roas_off))

dt_off = xr.DataTree.from_dict({"/posterior": xr.DataTree(combined_off)})

for channel in ["x1", "x2"]:
    pc, ax = plot_model_forest(dt_off, channel=channel, true_value=true_roas[channel])
    ax.set_title(
        f"Tight Calibration on an Overstated Experiment - Channel ${channel[0]}_{{{channel[1]}}}$",
        fontsize=18,
        fontweight="bold",
        y=1.05,
    )










../../_images/78585c67e878300d0a44fcd9f2928049d7051897332d118a277f686ccfb3e212.png

../../_images/0c48f8bdce5175bcf5f1930ab2a3c78204fffddaeaf5373727b29d582d1c33ad.png






Next, the probability of the correct ranking under each of the three experiments.








Hide code cell source





# Posterior probability that x2 has the higher ROAS: (model,)
prob_correct_off = (
    combined_off["roas"].sel(channel="x2") > combined_off["roas"].sel(channel="x1")
).mean(dim=("chain", "draw"))

fig, ax = plt.subplots(figsize=(10, 4), layout="constrained")

bars = ax.barh(
    prob_correct_off.coords["model"].to_numpy(), prob_correct_off.to_numpy(), color="C0"
)
ax.bar_label(
    bars,
    labels=[f"{prob:.1%}" for prob in prob_correct_off.to_numpy()],
    padding=4,
    fontsize=11,
)
ax.invert_yaxis()
ax.set(xlim=(0, 1.15), xlabel=r"$P(\mathrm{ROAS}_{x_2} > \mathrm{ROAS}_{x_1})$")
ax.set_title("Probability of the Correct Ranking", fontsize=18, fontweight="bold");










../../_images/4dc3bdc27e6b7fafdf20612e3cd728ed1b00eab10884bace6dae68baa236a94c.png






With the x1 experiment at \(140\) (\(+50\%\)), the posterior mean moves to about \(155\). That is about \(15\) above the experiment, close to the offset with the exact experiment (\(108\) against \(93\)). Its interval is still short and excludes the true \(93\) by a wide margin, so the model is now confidently wrong about the x1 ROAS. The decision mostly survives: the probability of the correct ranking drops from \(100\%\) to about \(80\%\), because the x1 experiment still sits below the x2 experiment at \(171\).


With the x1 experiment at \(187\) (\(+100\%\)), it sits above the x2 experiment. The posterior follows it to about \(190\), and the probability of the correct ranking falls to \(0.3\%\). The model now commits to the wrong order.


The two comparisons show what decides the result. The precision of the experiment decides how firmly the model follows it. The correctness of the experiment decides where the model goes, and so whether the decision comes out right. sigma should state the precision the experiment has. A tight sigma does not correct an experiment that is wrong.


We can look at the exact-experiment models as densities rather than intervals.






pc = azp.plot_dist(
    dt,
    col_wrap=1,
    visuals={
        "credible_interval": False,
        "point_estimate": False,
        "point_estimate_text": False,
    },
    figure_kwargs={
        "figsize": (12, 8),
        "sharex": True,
        "sharey": False,
        "layout": "constrained",
    },
)
fig = pc.viz["/"]["figure"].values.item()
axes = fig.axes
for ax, channel in zip(axes, ["x1", "x2"], strict=False):
    ax.axvline(true_roas[channel], color="black", linestyle="--", linewidth=2)
    ax.set(title=f"Channel ${channel[0]}_{{{channel[1]}}}$")
axes[-1].set(xlabel="ROAS")
color_map = pc.aes["color"]["mapping"]
handles = [
    Line2D([0], [0], color=str(color_map.sel(model=m).values), linewidth=3)
    for m in dt["posterior"].coords["model"].values
]
handles.append(Line2D([0], [0], color="black", linestyle="--", linewidth=2))
labels = [*list(dt["posterior"].coords["model"].values), "True ROAS"]
axes[0].legend(handles, labels, loc="center left", bbox_to_anchor=(1.02, 0.5))

fig.suptitle(
    "ROAS Posterior Densities by Model", fontsize=18, fontweight="bold", y=1.05
);








../../_images/980dfb203f20313dbe7a69513f9bb2b5c5fbf7f8615cd4777949dfbf5a62ed4a.png






The tight calibration produces a visibly narrower density for x1, the loose one a broader density on top of the naive one, and the moderate one sits between them.








What Calibration Costs You#


Calibration is not free. Forcing x1 down means the rest of the model has to
absorb the sales that x1 was previously explaining, usually through the
intercept and the trend. It is worth checking that the fit to the target has not
degraded in some pathological way.






for name, mmm_ in [("Naive", mmm_naive), ("Calibrated", mmm_calibrated)]:
    _ = mmm_.sample_posterior_predictive(
        X, extend_idata=True, combined=True, random_seed=rng
    )
    fig, _ = mmm_.plot.posterior_predictive()
    fig.suptitle(f"{name} Model - Posterior Predictive", fontsize=16, fontweight="bold")









Sampling: [y]







Sampling: [roas_calibration, y]







../../_images/d60ca21a6954cd8250060749ab4207aa6621b9895cf9be9aca8954fcabeb8478.png

../../_images/ba912cb2da96693936ba405d72496c6d1249b2d426341a1d45ea2de8e97a4f2a.png






Both models track the target closely. Calibration changed which component gets
credit for the sales, not whether the sales are explained.


This is the central lesson, and it is the reason predictive accuracy is a poor
way to choose between MMMs. Two models can fit the observed series equally well
and still give opposite budget advice. Only causal information, from experiments,
separates them. The Mitigating Unobserved Confounders in MMMs with Lift Test Likelihoods notebook makes the same point quantitatively
with cross-validated CRPS, where the uncalibrated model actually scores better
out of sample while being the more wrong of the two.






Calibrating on CPA Instead#


Many teams do not report ROAS. If your target variable is conversions rather than
revenue, the natural quantity is cost per acquisition, and that is what the
method does by default.


The only changes are the direction flag and the column name. We calibrate the
CPA table as written, with the default Normal, to see what the direction of
the table does; the warning below gives the conversion to use for a CPA that
comes from a lift test.






cpa_calibration_df = pd.DataFrame(
    {
        "channel": ["x1", "x2"],
        "cost_per_target": [1 / true_roas_x1, 1 / true_roas_x2],
        # 15% *relative* uncertainty on the CPA scale. A relative error is the
        # same on both scales (0.15 / roas divided by (1 / roas) ** 2 is
        # 0.15 * roas), so this table inverts exactly into the ROAS table.
        "sigma": [0.15 / true_roas_x1, 0.15 / true_roas_x2],
    }
)

mmm_cpa = MMM(**mmm_kwargs)
mmm_cpa.build_model(X, y)
mmm_cpa.add_original_scale_contribution_variable(var=vars_of_interest)

# target_per_cost defaults to False, and target_column defaults to
# "cost_per_target", so the CPA direction needs no extra arguments.
mmm_cpa.add_cost_per_target_calibration(
    data=X,
    calibration_data=cpa_calibration_df,
    name_prefix="cpa_calibration",
)

cpa_calibration_df













  
    

      
      channel
      cost_per_target
      sigma
    

  
  
    

      0
      x1
      0.010708
      0.001606
    

    

      1
      x2
      0.005834
      0.000875
    

  









Do not treat CPA and ROAS as interchangeable under a Normal


The point estimates are reciprocals, but the uncertainty is not. A Normal with
sigma on the ROAS scale implies an asymmetric, non-Normal uncertainty on the
CPA scale, and vice versa, and the next section shows that the direction also
decides on which side of the experiment the calibration can hold the model at all.
For a lift test the ROAS scale is the one to calibrate on, whichever way the
result was reported: the incremental outcome has the Normal error and the spend
is known. Convert a CPA and its standard error with
\(\text{ROAS} = 1 / \text{CPA}\) and
\(\sigma_{\text{ROAS}} = \sigma_{\text{CPA}} / \text{CPA}^2\), which is exact for
the delta-method standard error of a CPA and close while the relative error is
about \(30\%\) or less; a CPA given as bootstrap replicates or posterior draws is
inverted draw by draw. Taking the reciprocal of the standard error itself is
wrong.


With likelihood=Prior("LogNormal") the direction stops mattering for the
posterior, because that family reads sigma as a relative error, which is the
same number whichever way the table is written.




We fit this model too, as written, to show what the direction does. The same
experiment, the same \(15\%\), and a different answer from the ROAS table: the
forest plot in the next section puts the two side by side.






%%time

_ = mmm_cpa.fit(X, y, **sampler_config)









NUTS[nutpie]: [y_sigma, gamma_fourier, adstock_alpha, saturation_lam, saturation_beta, intercept_latent_process_raw_hsgp_coefs_offset, intercept_latent_process_raw_eta, intercept_latent_process_raw_ls, intercept_baseline]









CPU times: user 26.6 s, sys: 272 ms, total: 26.9 s
Wall time: 11.1 s













Choosing the Likelihood Family#


The ROAS calibration and the CPA calibration used the same experiment, the same
\(15\%\) relative sigma and the same Normal family, yet they pull x1 by
different amounts. What differs between them is the shape of the penalty.


Write \(r\) for the model’s ratio and \(t\) for the experiment’s value. A Normal
centered on \(r\) with the experiment observed charges



\[
-\log p(t \mid r) = \frac{(r - t)^2}{2\sigma^2} + \text{const}.
\]


As \(r\) grows the charge grows without bound. As \(r\) shrinks toward zero it
stops at \(t^2 / (2\sigma^2)\): with sigma at \(15\%\) of the target that is
about \(22\) nats, and no amount of further shrinking costs more than that. A
ratio of two positive quantities has zero on one side and infinity on the
other, so a Normal on it is a floor or a ceiling, never both. Which side is
capped is decided by how the table is written:


    

  • A CPA table (target_per_cost=False, the default) caps the charge as CPA goes
to zero, which is ROAS going to infinity. That is exactly where the confounder
pushes x1, and on that side the charge is nearly flat: moving x1 from
\(155\) to \(187\) costs the CPA table about \(2\) nats more, where the ROAS table
would charge about \(13\) nats more for the same move.
    • 

  • A ROAS table (target_per_cost=True) caps the charge as ROAS goes to zero and
is unbounded as ROAS grows. On this dataset the bias runs toward its unbounded
side, which is why the ROAS calibration pulled harder. With a channel the
model understates, the two tables would swap roles.
    • 



StudentT caps the same side lower than the Normal and grows only
logarithmically on the other, on purpose: heavier tails let the model walk away
from an experiment that may itself be wrong. LogNormal measures the error on
the log scale, so a \(2\times\) miss costs the same in either direction and the
cost keeps growing as \((\log(r/t))^2\). The figure below computes the
three charges from the same densities the calibration uses.






sigma_rel = 0.15  # sigma as a fraction of the experiment's value
nu = 4
ratio = np.geomspace(0.05, 20, 300)  # model ROAS / experiment ROAS


def penalty(family: str, model_ratio: np.ndarray, sigma: float) -> np.ndarray:
    """Nats charged relative to the experiment's own value, target normalized to 1."""

    def logp(mu: np.ndarray) -> np.ndarray:
        if family == "Normal":
            dist = pm.Normal.dist(mu=mu, sigma=sigma)
        elif family == "StudentT":
            dist = pm.StudentT.dist(nu=nu, mu=mu, sigma=sigma)
        else:
            # LogNormal: centered on log(ratio), width sigma / target on the log scale
            dist = pm.LogNormal.dist(mu=np.log(mu), sigma=sigma)
        return pm.logp(dist, 1.0).eval()

    return logp(np.ones(1)) - logp(model_ratio)















Hide code cell source





fig, axes = plt.subplots(ncols=2, sharey=True, figsize=(12, 5), layout="constrained")
for ax, (title, model_ratio) in zip(
    axes,
    [
        ("Table in CPA units (target_per_cost=False)", 1 / ratio),
        ("Table in ROAS units (target_per_cost=True)", ratio),
    ],
    strict=True,
):
    for family in ["Normal", "StudentT", "LogNormal"]:
        label = f"StudentT(nu={nu})" if family == "StudentT" else family
        ax.plot(
            ratio, penalty(family, model_ratio, sigma_rel), label=label, linewidth=2.5
        )
    ax.axhline(1 / (2 * sigma_rel**2), color="black", linestyle=":", label="Normal cap")
    ax.set(
        xscale="log",
        title=title,
        xlabel="model ROAS / experiment ROAS",
        ylim=(0, 40),
    )
axes[0].set(ylabel="penalty (nats)")
handles, labels = axes[0].get_legend_handles_labels()
fig.legend(handles, labels, loc="outside lower center", ncols=4)
fig.suptitle(
    "Calibration penalty by family, sigma at 15% of the experiment",
    fontsize=16,
    fontweight="bold",
);










../../_images/d2ab3569d6cb7e15a67aee42bf6d8c28cd743a28fed661804170e9c22d1d7e2b.png






Read the two panels against each other. The Normal and StudentT curves
flatten on the right in the CPA panel and on the left in the ROAS panel; the
LogNormal curve is the same parabola in both. Near the experiment the three
families charge similar amounts (at a factor of \(1.3\) either way the charges run
from about \(1.2\) to \(2\) nats at this sigma), so an existing sigma keeps its
meaning; they part ways in the tails.


sigma comes from the same number under every family: the standard error of the
value in the table, in the table’s units, which is the number an experiment
reports next to its point estimate. Normal uses it as its standard deviation.
StudentT uses it as its scale, which is not its standard deviation: with
nu=4 the spread is \(\sqrt{2} \approx 1.41\) times sigma, it has no finite
standard deviation for nu of \(2\) or less, and it approaches the Normal as
nu grows. Smaller nu means heavier tails, so the model walks away from the
experiment more readily. LogNormal turns it into a width on the log scale, sigma / target, the first-order translation of
an additive standard error into a relative one, so our \(15\%\) sigma becomes a
log-scale width of \(0.15\) in both tables. That is why the LogNormal term does
not care which way the table is written. The translation holds while the
relative error is small, as \(15\%\) is; at \(50\%\) the LogNormal reads the
experiment as allowing ratios from \(0.38\) to \(2.66\) times its value, where the
Normal allows \(0.02\) to \(1.98\) times.


The family is the likelihood argument, a
Prior as in add_lift_test_measurements. Normal
stays the default, so existing calibrated models are unchanged. We refit with
Prior("LogNormal") on both tables, and with Prior("StudentT", nu=4) on the
ROAS table, the same nu as in the figure.






%%time

roas_direction = {"target_column": "roas", "target_per_cost": True}

# label -> (calibration table, direction arguments, likelihood)
family_fits = {
    "CPA, LogNormal": (cpa_calibration_df, {}, Prior("LogNormal")),
    "ROAS, LogNormal": (roas_calibration_df, roas_direction, Prior("LogNormal")),
    "ROAS, StudentT": (roas_calibration_df, roas_direction, Prior("StudentT", nu=nu)),
}
mmm_family = {}

for label, (calibration_df, direction, likelihood) in family_fits.items():
    mmm_f = MMM(**mmm_kwargs)
    mmm_f.build_model(X, y)
    mmm_f.add_original_scale_contribution_variable(var=vars_of_interest)
    mmm_f.add_cost_per_target_calibration(
        data=X,
        calibration_data=calibration_df,
        name_prefix="calibration",
        likelihood=likelihood,
        **direction,
    )
    _ = mmm_f.fit(X, y, **sampler_config)
    mmm_family[label] = mmm_f









NUTS[nutpie]: [y_sigma, gamma_fourier, adstock_alpha, saturation_lam, saturation_beta, intercept_latent_process_raw_hsgp_coefs_offset, intercept_latent_process_raw_eta, intercept_latent_process_raw_ls, intercept_baseline]









NUTS[nutpie]: [y_sigma, gamma_fourier, adstock_alpha, saturation_lam, saturation_beta, intercept_latent_process_raw_hsgp_coefs_offset, intercept_latent_process_raw_eta, intercept_latent_process_raw_ls, intercept_baseline]









NUTS[nutpie]: [y_sigma, gamma_fourier, adstock_alpha, saturation_lam, saturation_beta, intercept_latent_process_raw_hsgp_coefs_offset, intercept_latent_process_raw_eta, intercept_latent_process_raw_ls, intercept_baseline]









CPU times: user 1min 21s, sys: 1.12 s, total: 1min 22s
Wall time: 34.5 s









We put the six posteriors on one forest plot per channel. Blue marks the ROAS table and red the CPA table, and the shade marks the family: full color for Normal, lighter for StudentT, lightest for LogNormal. The naive model is gray.








Hide code cell source





roas_by_family = {
    "naive": roas_naive,
    "ROAS, Normal": roas_calibrated,
    "ROAS, StudentT": get_roas(mmm_family["ROAS, StudentT"]),
    "ROAS, LogNormal": get_roas(mmm_family["ROAS, LogNormal"]),
    "CPA, Normal": get_roas(mmm_cpa),
    "CPA, LogNormal": get_roas(mmm_family["CPA, LogNormal"]),
}

# One color per entry of roas_by_family: gray for the naive model, then blues
# for the ROAS table and reds for the CPA table from matplotlib's "seismic"
# colormap, darkest for Normal, lighter for StudentT, lightest for LogNormal
seismic = plt.colormaps["seismic"]
family_colors = ["0.45", *[to_hex(seismic(x)) for x in (0.1, 0.25, 0.35, 0.9, 0.65)]]

combined_family = xr.concat(
    [roas.to_dataset() for roas in roas_by_family.values()], dim="model"
).assign_coords(model=list(roas_by_family))

dt_family = xr.DataTree.from_dict({"/posterior": xr.DataTree(combined_family)})

for channel, note in [("x1", "confounded"), ("x2", "unconfounded")]:
    pc, ax = plot_model_forest(
        dt_family,
        channel=channel,
        true_value=true_roas[channel],
        colors=family_colors,
    )
    ax.set_title(
        f"ROAS Posterior by Likelihood Family - Channel ${channel[0]}_{{{channel[1]}}}$ ({note})",
        fontsize=18,
        fontweight="bold",
        y=1.05,
    )










../../_images/3f24e75cfa023be7c929041ec2139976071fe70d12479be3279401ba4ed820f5.png

../../_images/63f7a3933a3107b923d6a07c031dfe2a48c6e0b2147e4ae561b77135abbe69a4.png






The posterior means in tabular form:






combined_family["roas"].mean(dim=("chain", "draw")).to_dataframe().unstack("channel")













  
    

      
      roas
    

    

      channel
      x1
      x2
    

    

      model
      
      
    

  
  
    

      naive
      196.983296
      142.227095
    

    

      ROAS, Normal
      154.616027
      149.027180
    

    

      ROAS, StudentT
      188.699727
      148.166818
    

    

      ROAS, LogNormal
      173.484459
      150.000712
    

    

      CPA, Normal
      187.191421
      150.370316
    

    

      CPA, LogNormal
      173.624683
      149.609846
    

  







Four things to read off the x1 panel.


First, the two Normal fits disagree, and the direction of the table is the only thing that differs between them.
The CPA table leaves x1 at about \(187\), close to the naive model’s \(197\), because its charge flattens on the side the data pushes toward.
The ROAS table holds it near \(155\) because that side is unbounded for it.


Next, the two LogNormal fits agree, up to sampling noise, as they must: once sigma is read as a relative error, the term is the same in both directions.
They land at about \(173\), between the two Normal fits.
On the high side the LogNormal charges less than the ROAS-table Normal: a \(1.7\times\) miss costs about \(6\) nats on the log scale against about \(11\) on the ROAS scale.
Against the CPA-table Normal it is far stronger, and it keeps getting stronger the further the model wanders, which the Normal does not.


Now the StudentT on the ROAS table. It lands at about \(189\), almost back at the naive model, and further from the truth than every other calibrated fit in this comparison.
This is what heavy tails are for: the sales data disagree with the experiment, and the StudentT lets the model side with the sales data.
Here the sales data are the ones that are wrong, because of the confounder, so this choice costs us.
A StudentT pays off in the opposite case, when the experiment is the doubtful part.


We do not fit a StudentT on the CPA table, because it would combine the two weak sides. A CPA table caps the charge as ROAS grows, which is where the confounder pushes x1, and the StudentT caps it lower than the Normal. The left panel of the penalty figure shows the result: the StudentT charge levels off near \(6\) nats however far ROAS grows, so it would barely hold x1. Use a StudentT on a ROAS table.


Finally, the x2 panel barely moves across families. Nothing pushes x2 far from its experiment, and within a factor of \(1.3\) of the experiment the families charge similar amounts.


None of this is a statement about which family is closer to the truth. On this dataset the ROAS-table Normal happens to land nearest, and it does so because the bias runs toward its unbounded side.
Pick the family by the error you believe the experiment has:


    

  • Additive error: Normal, the default, for an additive error in the units of the table, with a cap on the side you choose.
A ROAS read off a lift test is this case, since the incremental outcome has a Normal standard error and the test knows its spend.
So is a CPA from a lift test once it is inverted to a ROAS table, roas = 1 / cpa with sigma = sigma_cpa / cpa**2.
Here that turns the CPA table into the ROAS table exactly, and with it the \(155\) fit.
    • 

  • Doubtful experiment: Prior("StudentT", nu=...), for an additive error from an experiment that may itself be wrong, so that the model can walk away from it, as it did above.
    • 

  • Relative error: Prior("LogNormal"), so that a miss by a factor costs the same whichever way the table is written and the pull never flattens.
    • 



A ratio of two positive quantities reported with a relative standard error is the third case, which is why LogNormal is the family that closes the ceiling.
It assumes contributions that stay positive, as the HalfNormal channel priors here keep them.






In-Sample Fit Does Not Rank the Models#


The plots above were read by eye. The same point can be made with a proper score,
and it is worth doing because the result is counter-intuitive.


The Continuous Ranked Probability Score (CRPS) generalizes the mean absolute
error to probabilistic predictions. Instead of comparing a point forecast to the
observation, it compares the whole predictive distribution to it, so a model is
rewarded both for being close and for being honest about its uncertainty. Lower
is better. For the intuition, see the blog post
“Intuition behind CRPS”.


We compute it in sample, on the very dates the models were fit to, for the naive
model and the three calibrations of the sigma comparison.
We score posterior predictive draws taken with a fixed seed, so the scores are the
same on every run.






def in_sample_crps(mmm_: MMM) -> float:
    """In-sample CRPS of the posterior predictive against the observed target."""
    # Seeded posterior predictive draws of the target in original units: (chain, draw, date)
    y_pred = mmm_.sample_posterior_predictive(
        X, extend_idata=False, combined=False, random_seed=seed, progressbar=False
    )["y_original_scale"]

    # crps expects the sample dimension first: (sample, date)
    y_pred = (
        y_pred.stack(sample=("chain", "draw")).transpose("sample", "date").to_numpy()
    )
    return crps(mmm_.y, y_pred)


models_by_name = {
    "naive": mmm_naive,
    "calibrated, tight (5%)": mmm_sigma["tight (5%)"],
    "calibrated, moderate (15%)": mmm_calibrated,
    "calibrated, loose (60%)": mmm_sigma["loose (60%)"],
}

crps_by_model = {name: in_sample_crps(m) for name, m in models_by_name.items()}

pd.Series(crps_by_model, name="in-sample CRPS").round(2)









Sampling: [y]





Sampling: [roas_calibration, y]





Sampling: [roas_calibration, y]





Sampling: [roas_calibration, y]





naive                         6.93
calibrated, tight (5%)        9.06
calibrated, moderate (15%)    7.30
calibrated, loose (60%)       6.93
Name: in-sample CRPS, dtype: float64















Hide code cell source





fig, ax = plt.subplots(figsize=(10, 5), layout="constrained")

names = list(crps_by_model)
values = [crps_by_model[name] for name in names]

ax.axvline(
    values[0],
    color="C0",
    linestyle="--",
    linewidth=2,
    label="naive (reference)",
)

for i, value in enumerate(values):
    ax.scatter(value, i, s=200, color=f"C{i}", zorder=3)
    ax.annotate(
        f"{value:,.2f}",
        xy=(value, i),
        xytext=(0, 14),
        textcoords="offset points",
        ha="center",
        fontsize=11,
    )

ax.set_yticks(range(len(names)), names)
ax.invert_yaxis()
ax.legend(loc="upper left", bbox_to_anchor=(1.01, 1.0))
ax.set(xlabel="In-sample CRPS (lower is better)", ylabel=None)
ax.set_title(
    "In-Sample CRPS by Model",
    fontsize=18,
    fontweight="bold",
    y=1.05,
);










../../_images/9c817e41e731646db7898b3cdaa933bb9675990842c82548e464919bfd5d1d9c.png






The scores land in almost the opposite order to the causal accuracy:


    

  • calibrated, loose (\(60\%\)): \(6.93\)
    • 

  • naive: \(6.93\)
    • 

  • calibrated, moderate (\(15\%\)): \(7.30\)
    • 

  • calibrated, tight (\(5\%\)): \(9.06\)
    • 



The best score is shared by the naive model and the loose calibration, the model that stays closest to the naive fit and shares its inverted channel ranking.
The moderate calibration is close behind.
The tight calibration, the only model that recovers the correct ranking, comes last, about \(30\%\) worse than the best score.


This is not a bug, and it is not a sign that calibration damaged the model.
The calibration term pulls the channel contributions away from the split that best explains the observed sales, and the rest of the model absorbs the difference, so the fit to y degrades a little.


In-sample CRPS measures how well a model reproduces y. It does not measure whether the model splits the credit between x1 and x2 correctly, because y does not identify that split: many different attributions reproduce the same series about equally well. So the score can see the small loss in fit that calibration causes, but it cannot see the large gain in attribution.


Read the direction of this result carefully, because it is the practical warning of the whole notebook.
If you select models on predictive score here, you pick the loose calibration or the naive model and you get the budget advice backwards.
Predictive fit and causal correctness are different objectives, and on this dataset they point in opposite directions.
Mitigating Unobserved Confounders in MMMs with Lift Test Likelihoods shows the same pattern holding out of sample with cross-validated CRPS.


Three caveats on the numbers themselves.
They are in sample, computed on the very dates the models were fit to, so they flatter every model and should not be read as generalization performance.
The loose, moderate and naive scores lie close together, so do not read much into their order.
And the spread is small in absolute terms, roughly \(7\) to \(9\) against a target whose mean is about \(372\), so no model here fits badly.
The robust signal is the gap to the tight calibration: it scores about \(30\%\) worse than the best model, and it is the only one with the correct ranking.






Practical Guidance#


One calibration per experiment. Each row is an observation. Two rows for the
same channel mean two independent experiments, and the model will try to satisfy
both. If you have repeated readouts of the same experiment, pool them first.


Match the time window. The ratio is computed over the model’s full date range. If your experiment ran for four weeks in a year of data, the model-implied ROAS covers the whole year while the experiment covers four weeks.


What if ROAS varies inside the window? Then the calibration constrains only the average over the whole window, not the ROAS of any single period. ROAS can vary for two reasons. The channel’s effectiveness can change over time, and ROAS changes with the spend level, because the saturation curve flattens as spend grows. An experiment that ran in a peak quarter, or at an unusual spend level, measures the ROAS of that period and that spend, which can be far from the whole-window average. Calibrating it as a whole-window ratio pulls the average toward a number that was never an average. In that case use lift tests (Lift Test Calibration): each lift test states the spend level and the spend change it tested, and with time_varying_media=True each row can also carry a date column, so the model compares the experiment with its own response at that spend and at that date. The model in this notebook uses time_varying_media=False, and then the date column is ignored.


Set sigma from the experiment. Use the standard error the analysis produced. In this notebook the tightest sigma gave the best result, but only because every experiment here reports the true ROAS and the tight scenario describes a very precise one. In the When the Experiment Is Biased subsection, the same tight sigma on an experiment that overstated the true ROAS held the model about as firmly near the wrong value. If you have no uncertainty estimate at all, prefer a deliberately loose sigma over a tight guess, because a loose sigma does not add certainty that the experiment does not have. A loose sigma has a cost: on this dataset the loose fit kept the wrong ranking. So when sigma is a guess, refit at a few values, as we did here, and report whether the decision changes between them.


Watch for conflict. If the calibration and the sales data disagree sharply,
the sampler will tell you: divergences, poor r_hat, or a posterior that sits
between the two and satisfies neither. Treat that as a signal to investigate,
not something to sample harder through. The usual culprits are a spend column
that is not really spend, a channel naming mismatch, or an experiment measuring a
different outcome from the model target. Poor r_hat can also have nothing to do
with calibration: the naive model in this notebook had a stuck chain without any
calibration term. Check the uncalibrated model first, so you know which problem
you are looking at.


Combining with lift tests. add_lift_test_measurements and
add_cost_per_target_calibration can both be applied to the same model. Just be
aware that if a lift test and a ROAS estimate come from the same experiment,
adding both double-counts the evidence and will overstate your certainty.


Pick the family by the error you believe in. The default Normal reads
sigma as an additive error and caps the charge on one side of the experiment;
which side depends on target_per_cost, so decide it on purpose. A lift test’s
readout has its additive error on the ROAS scale, so a CPA from one is inverted
to a ROAS table (roas = 1 / cpa, sigma = sigma_cpa / cpa**2) and calibrated
with the Normal. Prior("StudentT", nu=...)
caps that side lower still and grows only logarithmically on the other, so the
model can walk away from an experiment that may be wrong.
Prior("LogNormal") reads sigma as a relative error, charges a miss by
a factor the same whichever way the table is written, and never flattens; it is
the family for a ratio reported with a relative standard error. See
Choosing the Likelihood Family.


Multiplicative models (link="log"). The same call works for a log-link model, with one difference in what it constrains.
Under link="log" a channel’s channel_contribution_original_scale is the factor \(s\,e^{m_{t,c}}\), not a level change, so the calibration targets the counterfactual increment \(\Delta_{t,c} = s \left(e^{\mu_t} - e^{\mu_t - m_{t,c}}\right)\) instead, which it registers as the deterministic channel_incremental_contribution.
For a channel that no mu_effect reads, the calibrated ratio is exactly mmm.incrementality.contribution_over_spend(frequency="all_time", central_tendency="mean"), draw for draw.
That "mean" is the calibration’s default central_tendency, and it matters under this link: \(\Delta\) is on the conditional-median scale, while an experiment measures the expected realized lift, which is larger by \(e^{\sigma^2/2}\) (\(1.03\) at \(\sigma = 0.25\), \(1.13\) at \(0.5\)), so the calibration applies that factor per cell of the likelihood scale.
central_tendency="median" drops the factor and matches the default scale of contribution_over_spend, but no experiment measures that scale.
Two more things differ from the identity link.
The increment scales with the baseline \(e^{\mu_t}\), so the time-window mismatch above matters more: a lift test run in a peak season measures a larger ROAS than the full-window average.
And the per-channel increments do not add up to the total media contribution, because each removes one channel while the total removes all of them at once, so do not calibrate per-channel figures obtained by splitting a total-media number.
Under both links the calibration constrains a channel’s direct path only: an experiment on a channel that also feeds a mu_effect (a funnel mediator, say) measures the direct plus the mediated effect, so calibrate direct-only channels or read the target as a direct-path ROAS.






Conclusion#


Calibration turns an experiment result into an extra likelihood term.
The model keeps fitting sales, and simultaneously has to explain the experimental observation.
Where the sales data underdetermines the answer, that second term is what moves it.


Concretely, in this notebook:


    

  • The naive model got channel x1 badly wrong, about \(197\) against a true \(93\), because of an unobserved confounder, and it inverted the ranking of the two channels.
No amount of extra sales data would have fixed either problem.
    • 

  • add_cost_per_target_calibration with target_per_cost=True removed roughly \(40\%\) of that error at a moderate sigma, and at a tight sigma it restored the correct ranking, which is the part that changes a budget decision.
It did not recover the true ROAS exactly: even the tight calibration left the truth just outside the \(94\%\) interval for x1.
One experimental observation against a full sales series pulls the posterior, it does not overrule it.
    • 

  • sigma is the trust dial, and here it decided the outcome.
At \(15\%\) the ranking stayed inverted, and at \(5\%\) it flipped to the correct order.
The same \(5\%\) on an experiment that overstated x1 moved the posterior with it: by \(50\%\), to about \(155\) with an interval that excludes the true \(93\), and by \(100\%\), far enough to commit to the wrong ranking.
That is a lot of leverage for one number, which is exactly why it should come from the experiment’s own uncertainty rather than from tuning until the answer looks right.
    • 

  • In-sample CRPS ranked the models in roughly the reverse order of their causal accuracy.
You cannot choose between these models on predictive fit.
    • 

  • The same experiment at the same sigma pulled x1 by different amounts depending on whether the table was written as CPA or as ROAS, because a Normal on a ratio caps the charge on one side and the direction picks which.
likelihood=Prior("LogNormal") removed that dependence: its charge is symmetric and unbounded on the log scale, and the two tables gave the same posterior.
    • 

  • likelihood=Prior("StudentT", nu=4) on the ROAS table let the model walk away from the experiment, back to about \(189\), close to the naive model.
Heavy tails help when the experiment is the doubtful part. Here the confounded sales data were the doubtful part, so the StudentT gave up most of what the calibration gained.
    • 



If your experiment needs to constrain the model more firmly than it does here, the whole-period ratio may be the wrong instrument: a lift test anchored to the dates the experiment actually ran is a sharper constraint than an average over the full history.


Where to go next:







%load_ext watermark
%watermark -n -u -v -iv -w -p pymc_marketing,pytensor









Last updated: Thu, 05 Nov 2026

Python implementation: CPython
Python version       : 3.12.11
IPython version      : 9.15.0

pymc_marketing: 1.3.0
pytensor      : 3.3.7

arviz         : 1.2.0
arviz_plots   : 1.2.0
graphviz      : 0.21
matplotlib    : 3.11.1
numpy         : 2.4.6
pandas        : 3.0.5
pymc          : 6.4.0
pymc_extras   : 0.16.0
pymc_marketing: 1.3.0
seaborn       : 0.13.2
xarray        : 2026.4.0

Watermark: 2.6.0