Splitting MMM Channel Impact Across Campaigns

Splitting a channel’s Marketing Mix Model impact across campaigns, using a 100-year-old idea from insurance

Executive summary

The question. A Marketing Mix Model (MMM) tells you how much incremental revenue the Paid Social channel earned, day by day, with uncertainty. It cannot tell you which campaign earned it. The platform dashboard reports conversions per campaign, but those numbers over-claim, they are not incremental, and the small campaigns among them swing wildly from week to week.

The idea. Treat each campaign’s dashboard-based estimate as a noisy witness and decide how much to believe each one. The rule comes from insurance pricing (Bühlmann credibility) and is the same thing as Bayesian partial pooling: blend the campaign’s own estimate with the channel average, and trust the campaign’s own estimate more the more unique people it reached. Reach is evidence.

The three steps.

1. Build each campaign’s own noisy estimate \(y_c\): split the MMM’s channel total by click share and express it per reach-day.

2. Let the data set the trust dial. A seven-row hierarchical model learns how different campaigns really are and how noisy one campaign’s estimate is, and returns a trust weight \(Z_c\) between 0 and 1 for each campaign, with uncertainty.

3. Hand the MMM’s daily credit back to the campaigns. Each day’s channel credit is split in proportion to trust-adjusted effectiveness times that day’s reach. The pieces add back to the MMM total exactly, every day.

What you need. A fitted PyMC-Marketing MMM (its channel_contribution posterior), the daily per-campaign platform export (spend, impressions, reach, clicks), and one de-duplicated reach number per campaign over its flight. The expanding-window reach series described in the data section is useful for the audience growth curve, but the model itself only needs the flight total.

The results. The data are synthetic, so we can grade the split against a hidden truth. Against the click split it disciplines, the credibility split improves the campaign ranking (Spearman 0.79 to 0.89), the share error (MAE 0.043 to 0.035) and the probabilistic score (CRPS 0.043 to 0.025), and 5 of the 7 true shares fall inside its 94% interval. It is the best probabilistic forecast of every method tested. We also say plainly what it does not win: on point accuracy alone, a plain delivery split is still slightly better (MAE 0.032).

Checks and extensions. We inherit the MMM’s adstock so a campaign’s credit can outlive its flight, sweep the priors, test unique reach against reach-days as the evidence measure, fit two compositional variants of the likelihood, and close with what to do when you cannot get unique reach at all.

What it cannot do. It cures noise, not bias: if clicks are a poor proxy for incrementality, the split inherits that slant. It is a report card on credit already earned, not advice on where to move budget.

Symbols are collected in the notation table at the start of Step 1, and every symbol is defined again where it first appears.

The problem

Your Marketing Mix Model has done its job. It tells you, with honest uncertainty, that Paid Social drove this much incremental revenue, day by day. You bring it to the review, and the first question is the one the MMM cannot answer:

“Which campaign?”

So you open the platform dashboard. Every campaign claims conversions. Add them up and they total far more than the MMM says the whole channel earned. The always-on brand campaign looks steady, a weekend burst commands attention, a tiny niche campaign swings from hero to flop week to week, and the broad prospecting push that only launched in April has barely four weeks of history. On top of that, most of those platform conversions are rule-based and non-incremental: plenty of those people would have converted anyway.

Strip the noise away and the problem is simple. You have one number you trust (the MMM’s channel total) and seven shaky signals (the per-campaign platform numbers). The question is not whether to use the shaky signals, because they are all you have. It is how much to believe each one.

The old idea we are going to borrow

This is the problem an insurance actuary solves every day.

Picture pricing car insurance. A driver comes in with one claim-free year. Do you give them the flawless-driver discount? Of course not, because one year is too little to tell skill from luck. So you blend their record with the average driver, and you lean on their own history more and more as it piles up:

\[\text{best guess}=Z\cdot(\text{their record})+(1-Z)\cdot(\text{average driver}), \qquad Z=\frac{n}{n+k}.\]

That blend has a name, credibility, and a single dial, Z, that runs from 0 (“I will trust the average”) to 1 (“I will trust you”). It turns on one thing: how much evidence you have (\(n\)), measured against a threshold \(k\) that says how much evidence counts as “enough”. A campaign that has reached exactly \(k\) people is believed half-way; one that has reached far more than \(k\) people is mostly believed; one that has reached far fewer is mostly replaced by the average.

In this notebook the threshold comes out at about 37,000 people, roughly the size of an average campaign here (the fitted value is reported exactly in Step 2). The campaign that reached 62,000 people ends up 61% believed, so we take most of its own story and fill the rest with the group average. The campaign that reached 18,000 ends up 34% believed, mostly the average until it earns more trust. Nobody is ignored, and nobody is blindly believed.

The whole idea of this notebook, in one sentence: swap “years of driving history” for “how many people the campaign reached”, and use that dial to decide how much of each campaign’s own story to believe. Reach is evidence.

What you get at the end is a per-campaign, per-day split of the MMM’s channel credit that (a) always adds back up to the MMM total, (b) cannot be hijacked by small, loud campaigns, and (c) carries an uncertainty band that says how sure the method itself is. To be upfront about the scale of the win: on this data the measurable gains are a better campaign ranking, better calibrated uncertainty and a trust weight you can justify, not a jump in point accuracy over a plain delivery split. And to be upfront about what the answer is: it is a report card on credit already earned, not reallocation advice; the last section spells out that boundary and everything else the method cannot do.

Before you start

What you need. A fitted PyMC-Marketing MMM whose posterior contains the channel contribution (channel_contribution, see the MMM example notebook), the daily per-campaign export from the ad platform (spend, impressions, reach, clicks), and one de-duplicated reach number per campaign over its flight. The data section shows how to get that number, and the audience growth curve behind it, with an expanding-window query. The section “What if you cannot get unique reach?” near the end covers the case where you cannot.

What you will do. Three steps (build each campaign’s noisy own story, let the data set the trust dial, hand the MMM’s daily credit back to the campaigns), then grade the result against a hidden truth, extend it with the MMM’s own carryover, and stress-test it before reading the limits.

What it costs. About twenty small PyMC fits of a seven-row model. A few minutes on a laptop.

The data. Everything here is synthetic. A ground-truth file of per-campaign contributions is loaded alongside the inputs and used only to grade the result at the end, never to fit anything.

For the statistically minded: the exact mapping

This is Bühlmann credibility, which blends a policyholder’s own experience with the collective mean, $\(\text{premium}=Z\cdot(\text{own})+(1-Z)\cdot(\text{collective}),\qquad Z=\frac{n}{n+k},\)$ and it is exactly hierarchical-Bayes partial pooling.

An actuary routinely uses one quantity to scale the premium and a different one to set the credibility weight, and that is what happens here. A person re-reached on a second day adds delivery volume but no new information about the campaign, so reach-days carry the exposure and unique people reached carry the evidence.

Bühlmann credibility (actuarial)	This MMM disaggregation
Policyholder	Campaign \(c\)
Portfolio, the collective	The whole social channel
Own loss experience (noisy)	Naive contribution per reach-day (social split by click share)
Earned exposure, in policy-years	Reach-days \(R_c=\sum_t \text{reach}_{ct}\), the volume the credit is spent over
Volume measure that sets credibility	Unique people reached \(n_c\), from the de-duplicated reach report
Loss per unit exposure	Contribution per reach-day
Collective mean	Channel-average log-effectiveness \(\mu\)
Credibility factor \(Z_c=\frac{n_c}{n_c+k}\)	Trust in a campaign’s own signal
Credibility premium	Trust-adjusted rate per reach-day, spent over daily reach-days

Writing the sampling variance of campaign \(c\)’s direct estimate as \(\psi_c = s^2/n^{\text{rel}}_c\), where \(n^{\text{rel}}_c\) is its unique reach relative to the average campaign and \(s^2\) is the noise variance at average reach (both defined properly in Step 1), $\(Z_c=\frac{n^{\text{rel}}_c}{n^{\text{rel}}_c+k}=\frac{\tau^2}{\tau^2+\psi_c},\qquad k=\frac{s^2}{\tau^2},\)$ which is one minus the shrinkage factor of a Fay-Herriot small-area model. Fay-Herriot is the workhorse of small-area estimation: each area’s noisy direct estimate is shrunk toward a model-based mean, more strongly the larger its known sampling variance. The only twist here is that the sampling variances are tied to reach rather than supplied by hand.

from dataclasses import dataclass

import arviz as az
import arviz_plots as azp
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import pytensor.tensor as pt
import pytensor.xtensor as ptx
import seaborn as sns
import xarray as xr
from IPython.display import display
from matplotlib.lines import Line2D
from scipy.stats import norm, spearmanr

from pymc_marketing.metrics import crps, per_observation_crps
from pymc_marketing.mmm.transformers import geometric_adstock
from pymc_marketing.paths import data_dir

az.style.use("arviz-darkgrid")
plt.rcParams["figure.figsize"] = [10, 6]
plt.rcParams["figure.dpi"] = 100
plt.rcParams["figure.facecolor"] = "white"
az.rcParams["stats.ci_kind"] = "hdi"
az.rcParams["stats.ci_prob"] = 0.94

CI_PROB = az.rcParams["stats.ci_prob"]
seed = 42
rng: np.random.Generator = np.random.default_rng(seed=seed)

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

Before touching any data, here is the whole method in one picture. The trust dial \(Z\) is the share of a campaign’s own story that we believe. It rises with the number of people the campaign has reached, \(n\), and the threshold \(k\) decides how fast: at \(n = k\) a campaign is believed exactly half-way. Everything that follows is about measuring \(n\) properly, letting the data choose \(k\), and then spending the MMM’s credit according to the resulting blend.

people_grid = np.geomspace(5e3, 5e5, 400)

fig, ax = plt.subplots()
for threshold, color in zip((10_000, 40_000, 160_000), ("C0", "C1", "C2"), strict=True):
    ax.plot(
        people_grid,
        people_grid / (people_grid + threshold),
        color=color,
        lw=2.2,
        label=f"$k$ = {threshold / 1e3:.0f}k people",
    )
    ax.plot(threshold, 0.5, "o", color=color, ms=9, zorder=4)
ax.axhline(0.5, ls=":", color="gray", lw=1)
ax.annotate(
    "half-trust point $n = k$:\na campaign this big is\nbelieved 50/50",
    (40_000, 0.5),
    xytext=(20, -80),
    textcoords="offset points",
    ha="center",
    fontsize=10,
    arrowprops={"arrowstyle": "->", "color": "gray"},
)
ax.text(6e3, 0.06, "trust the\naverage", fontsize=11, ha="left", style="italic")
ax.text(8e4, 0.90, "trust the\ncampaign", fontsize=11, ha="right", style="italic")
ax.set(
    xscale="log",
    xlim=(5e3, 5e5),
    ylim=(0, 1),
    xlabel="people reached by the campaign, $n$ (log scale)",
    ylabel="trust dial $Z = n / (n + k)$",
    title="How much of a campaign's own story to believe",
)
ax.legend(loc="lower right");

Fig. The credibility dial. The x axis is the number of people a campaign has reached (log scale) and the y axis is the share of its own story we believe, \(Z\). Each curve is one value of the evidence threshold \(k\); the dot marks the half-trust point where a campaign has reached exactly \(k\) people. Small campaigns sit on the left and are mostly replaced by the group average, large ones sit on the right and are mostly believed. Step 2 fits \(k\) from the data and puts the seven real campaigns on this curve.

Load the data

Seven campaigns ran on one social channel across a 120-day window: six veterans plus a broad prospecting push that only kicked off on April 1st. We load five inputs and one grading file:

• credibility_campaign_daily_metrics.csv: the everyday delivery export any analyst has, one row per campaign per day with spend, impressions, reach and clicks.

• credibility_campaign_cumulative_reach.csv: the expanding-window unique reach per campaign, a second report explained in the data section below. Its last row is the de-duplicated reach over the whole flight, the one number the model needs.

• credibility_audience_overlap_matrix.csv: the pairwise audience-overlap matrix between campaigns.

• credibility_mmm_input_data.csv: the dates the MMM was fit on.

• credibility_channel_contribution_posterior.nc: the MMM’s social-contribution posterior, the full range of plausible values for the channel’s daily impact rather than a point estimate. That posterior is the trustworthy number we are going to divide up. In your own work this is mmm.idata.posterior["channel_contribution"] selected on the channel.

• credibility_ground_truth_campaign_contributions.csv: the simulated ground truth, used only at the end to grade our work, never to fit.

@dataclass
class CampaignData:
    """Delivery telemetry, the MMM social posterior, and the held-out ground truth."""

    metrics: pd.DataFrame
    cum_reach: pd.DataFrame
    overlap: pd.DataFrame
    mmm: pd.DataFrame
    social: xr.DataArray
    truth: pd.DataFrame  # oracle, used only for grading

    @property
    def campaigns(self) -> list[str]:
        """Campaign names, in the order reported by the expanding-window query."""
        return list(self.cum_reach.columns)


def load_campaign_data() -> CampaignData:
    """Read the campaign telemetry, the MMM social posterior, and the ground truth."""
    contribution = xr.open_dataset(
        data_dir / "credibility_channel_contribution_posterior.nc"
    )["channel_contribution"]
    return CampaignData(
        metrics=pd.read_csv(
            data_dir / "credibility_campaign_daily_metrics.csv", parse_dates=["date"]
        ),
        cum_reach=pd.read_csv(
            data_dir / "credibility_campaign_cumulative_reach.csv",
            index_col=0,
            parse_dates=True,
        ),
        overlap=pd.read_csv(
            data_dir / "credibility_audience_overlap_matrix.csv", index_col=0
        ),
        mmm=pd.read_csv(
            data_dir / "credibility_mmm_input_data.csv", parse_dates=["date"]
        ),
        social=contribution.sel(channel="social"),
        truth=pd.read_csv(
            data_dir / "credibility_ground_truth_campaign_contributions.csv",
            index_col=0,
            parse_dates=True,
        ),
    )


data = load_campaign_data()
campaigns = data.campaigns
colors = {campaign: f"C{i}" for i, campaign in enumerate(campaigns)}
channel_total = float(data.social.sum("date").mean())

print(
    f"{len(campaigns)} campaigns, {data.metrics['date'].nunique()} days from "
    f"{data.metrics['date'].min():%Y-%m-%d} to {data.metrics['date'].max():%Y-%m-%d}.\n"
    f"MMM social contribution over the window: {channel_total:,.0f} "
    "(posterior mean of the total)."
)

7 campaigns, 120 days from 2025-01-01 to 2025-04-30.
MMM social contribution over the window: 281,777 (posterior mean of the total).

Exploratory data analysis

The data you already have

Every ad platform, and its reporting API, hands you a table like this: one row per campaign per day, with delivery, attention, and cost. Here is one mid-flight day in February, when six of the seven campaigns are live.

Two words recur from here on, so let us pin them down. Delivery is the volume the platform serves: impressions and daily reach, the platform’s side of the story. Exposure is the actuarial term for the volume base a rate is charged over; in this notebook that base is reach-days, so “exposure” always means reach-days below. They are close cousins (both count volume, not people) but they play different roles: delivery is what the dashboard reports, exposure is what the credit gets spent over.

platform_export = (
    data.metrics.loc[
        data.metrics["impressions"] > 0,
        ["date", "campaign", "reach", "impressions", "clicks", "spend"],
    ]
    .assign(
        ctr_pct=lambda df: (100 * df["clicks"] / df["impressions"]).round(2),
        cpm=lambda df: (1000 * df["spend"] / df["impressions"]).round(2),
    )
    .drop(columns="spend")
)
mid_flight = sorted(data.metrics["date"].unique())[55]
platform_export.query("date == @mid_flight").reset_index(drop=True)

	date	campaign	reach	impressions	clicks	ctr_pct	cpm
0	2025-02-25	brand_always_on	13911	26919	408	1.52	13.43
1	2025-02-25	spring_promo	9437	18430	371	2.01	15.00
2	2025-02-25	lookalike_conversions	3386	4006	70	1.75	17.13
3	2025-02-25	video_product_launch	7426	14472	150	1.04	10.43
4	2025-02-25	flash_sale_burst	11227	41245	769	1.86	15.30
5	2025-02-25	niche_interest	1859	2342	52	2.22	19.26

Aggregated over each campaign’s flight, the same export gives the summary any media report starts from. Two columns deserve care because the whole notebook turns on the difference between them:

• reach_days is the daily reach column summed over days. It counts the same person again every day the campaign touches them.

• unique_reach is the number of distinct people reached by the end of the window. It comes from a second report, the expanding-window query, which we explain a few cells down.

Their ratio, reach_days_per_person, is how many days the average reached person saw the campaign.

def campaign_totals(metrics: pd.DataFrame, campaigns: list[str]) -> pd.DataFrame:
    """Flight totals per campaign: clicks, impressions, spend, reach, active days."""
    return (
        metrics.groupby("campaign")
        .agg(
            clicks=("clicks", "sum"),
            impressions=("impressions", "sum"),
            spend=("spend", "sum"),
            reach_days=("reach", "sum"),
            active_days=("active", "sum"),
            peak_daily_reach=("reach", "max"),
        )
        .reindex(campaigns)
    )


def final_unique_reach(cum_reach: pd.DataFrame, campaigns: list[str]) -> pd.Series:
    """Unduplicated people reached by the end of the window, the evidence measure."""
    return cum_reach.ffill().iloc[-1].reindex(campaigns).astype(float)


def campaign_summary(data: CampaignData) -> pd.DataFrame:
    """Flight-level delivery, cost and efficiency, one row per campaign."""
    totals = campaign_totals(data.metrics, data.campaigns)
    unique_reach = final_unique_reach(data.cum_reach, data.campaigns)
    return pd.DataFrame(
        {
            "active_days": totals["active_days"],
            "spend": totals["spend"],
            "impressions": totals["impressions"],
            "clicks": totals["clicks"],
            "reach_days": totals["reach_days"],
            "unique_reach": unique_reach,
            "ctr_%": 100 * totals["clicks"] / totals["impressions"],
            "cpm": 1000 * totals["spend"] / totals["impressions"],
            "cpc": totals["spend"] / totals["clicks"],
            "impressions_per_person": totals["impressions"] / unique_reach,
            "reach_days_per_person": totals["reach_days"] / unique_reach,
        }
    ).sort_values("unique_reach", ascending=False)


summary = campaign_summary(data)
summary.round(2)

	active_days	spend	impressions	clicks	reach_days	unique_reach	ctr_%	cpm	cpc	impressions_per_person	reach_days_per_person
brand_always_on	120	36945.66	2952795	40073	1541137	62386.0	1.36	12.51	0.92	47.33	24.70
lookalike_conversions	120	14175.89	825661	16891	597676	42733.0	2.05	17.17	0.84	19.32	13.99
spring_promo	66	9775.35	699815	15985	438462	41719.0	2.28	13.97	0.61	16.77	10.51
flash_sale_burst	15	5184.37	370310	10084	139620	34802.0	2.72	14.00	0.51	10.64	4.01
video_product_launch	71	11884.28	1064727	10628	519340	33427.0	1.00	11.16	1.12	31.85	15.54
broad_prospecting	30	283.05	33994	1294	32921	24140.0	3.81	8.33	0.22	1.41	1.36
niche_interest	120	6619.94	354716	7709	254475	17987.0	2.17	18.66	0.86	19.72	14.15

The seven campaigns are deliberately unalike: spend ranges from $283 to $37k, flights from 15 days to the full 120, CTR from 1.0% to 3.8%, and the average reached person saw brand_always_on on 25 different days but broad_prospecting on 1.4. That spread is what makes a single flat trust level indefensible: these campaigns simply do not carry the same amount of evidence.

fig, axes = plt.subplots(
    nrows=3, ncols=1, sharex=True, figsize=(12, 9), layout="constrained"
)
for ax, column, label in zip(
    axes,
    ["impressions", "reach", "clicks"],
    ["daily impressions (k)", "daily reach (k)", "daily clicks (k)"],
    strict=True,
):
    wide = data.metrics.pivot(index="date", columns="campaign", values=column).reindex(
        columns=campaigns
    )
    for campaign in campaigns:
        ax.plot(
            wide.index,
            wide[campaign] / 1e3,
            color=colors[campaign],
            lw=1.4,
            label=campaign,
        )
    ax.set(ylabel=label)
axes[0].set(title="Daily delivery per campaign")
axes[-1].set(xlabel="date")
axes[-1].xaxis.set_major_locator(mdates.MonthLocator())
axes[-1].xaxis.set_major_formatter(mdates.DateFormatter("%b"))
fig.legend(
    *axes[0].get_legend_handles_labels(),
    loc="outside center right",
);

Fig. Daily impressions, reach and clicks per campaign. The flight structure is visible: brand_always_on, lookalike_conversions and niche_interest run all 120 days, spring_promo and video_product_launch run in blocks, flash_sale_burst is a 15-day spike, and broad_prospecting only starts on April 1st. Reach and impressions move together within a campaign, which is the first hint that daily reach is a delivery measure rather than an audience measure.

delivered = data.metrics.query("impressions > 0 and reach > 0").assign(
    ctr_pct=lambda df: 100 * df["clicks"] / df["impressions"],
    daily_frequency=lambda df: df["impressions"] / df["reach"],
)
order = delivered.groupby("campaign")["ctr_pct"].median().sort_values().index.tolist()
palette = [colors[campaign] for campaign in order]

fig, (ax_ctr, ax_freq) = plt.subplots(
    ncols=2, sharey=True, figsize=(12, 7), layout="constrained"
)
for ax, column, title in zip(
    [ax_ctr, ax_freq],
    ["ctr_pct", "daily_frequency"],
    ["Daily CTR (%)", "Daily impressions per person"],
    strict=True,
):
    sns.boxplot(
        data=delivered,
        x=column,
        y="campaign",
        order=order,
        hue="campaign",
        hue_order=order,
        palette=palette,
        legend=False,
        fliersize=2,
        ax=ax,
    )
    ax.set(title=title, xlabel=None, ylabel=None)

Fig. Distribution over delivery days of each campaign’s click-through rate and daily frequency, ordered by median CTR. Days with no impressions are excluded, since 298 of the 840 campaign-day cells carry no delivery. Two things matter for what follows. First, CTR differs by campaign by a factor of nearly four, so a click-share split is not a neutral allocation. Second, the small campaigns are the jumpy ones: broad_prospecting and flash_sale_burst have the widest day-to-day CTR spread, which is exactly the noisiness credibility is designed to discipline.

The data you need, and how to build it

The export above is genuinely useful: delivery (reach, impressions), attention (clicks and CTR), and cost (CPM). It is also all most teams ever look at. For our question, who deserves the channel’s credit, it has one blind spot that matters: its reach is duplicated across days. Summing the column gives impressions-flavoured volume, not people, so it cannot tell you when a campaign has saturated its audience, the point where more spend just re-hits the same faces. And “how many people has this campaign reached” is precisely the evidence measure the credibility dial needs.

Unduplicated audience growth is not a column you can export directly, but you can build it from the same platform or API with a particular query pattern and no extra tracking. Call it the expanding window. Do not pull one report of daily reach. Pull a series of reports where you fix the start date and push the end date out one day at a time: [start to day 0], then [start to day 1], then [start to day 2], through [start to day N]. The platform de-duplicates within each report, so each one is the true unique reach so far, and the difference between consecutive reports is the count of genuinely new people reached that day. Stitch them together and you get a clean audience-growth curve instead of a duplicated daily count. If all you need is the flight total, one report over the whole flight window is enough, because the platform de-duplicates within a single report’s date range; the expanding-window series adds the growth curve on top.

fig, ax = plt.subplots(figsize=(12, 7), layout="constrained")
for campaign in campaigns:
    ax.plot(
        data.cum_reach.index,
        data.cum_reach[campaign] / 1e3,
        lw=1.9,
        color=colors[campaign],
        label=campaign,
    )
ax.set(
    title="Unduplicated audience growth (expanding-window reach)",
    xlabel="date",
    ylabel="cumulative unique reach (thousands)",
)
ax.xaxis.set_major_locator(mdates.MonthLocator())
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b"))
fig.legend(loc="outside center right");

Fig. The expanding-window result: cumulative unique reach per campaign. Most curves are C-shaped, steep while the campaign keeps finding new people and then flattening as it exhausts its audience, and that flattening point is the audience saturation the raw daily table cannot reveal. Scales differ enormously, from brand_always_on saturating near 62k unique users to niche_interest near 18k. The exception is broad_prospecting: it launched on April 1st against a very large audience and is still climbing when the window closes, so it has no saturation knee, just a truncated growth curve.

Put the two reach numbers side by side and the difference is stark.

volume = pd.DataFrame(
    {
        "reach_days": campaign_totals(data.metrics, campaigns)["reach_days"],
        "unique_reach": final_unique_reach(data.cum_reach, campaigns),
    }
).sort_values("unique_reach")
thousands = volume / 1e3
positions = np.arange(len(volume))

fig, ax = plt.subplots(figsize=(10, 6), layout="constrained")
ax.barh(
    positions + 0.2,
    thousands["reach_days"],
    height=0.38,
    color="C0",
    label="reach-days (exposure the credit is spent over)",
)
ax.barh(
    positions - 0.2,
    thousands["unique_reach"],
    height=0.38,
    color="C1",
    label="unique people reached (evidence that sets trust)",
)
for position, (reach_days, unique_reach) in enumerate(volume.to_numpy()):
    ax.annotate(
        f"{reach_days / unique_reach:.1f}x",
        (reach_days / 1e3, position + 0.2),
        xytext=(6, 0),
        textcoords="offset points",
        va="center",
        fontsize=12,
    )
# leave room on the right for the ratio annotations, and never clip a bar
ax.set(
    xscale="log",
    xlim=(thousands.to_numpy().min() / 2, thousands.to_numpy().max() * 2),
    yticks=positions,
    yticklabels=volume.index,
    xlabel="thousands (log scale)",
    title="Daily reach summed over days is not a count of people",
)
fig.legend(loc="outside lower center", ncol=2);

Fig. Reach-days against unique people reached, with their ratio annotated. The gap runs from 1.4x for broad_prospecting, which barely revisits anyone, to 24.7x for brand_always_on, which spends most of its delivery re-reaching people it already has. Summing the daily reach column therefore measures delivery volume, not audience size.

Two quantities, kept apart

The rest of the notebook uses these two words in exactly one sense each, because conflating them is the easiest way to get the arithmetic wrong.

• Unique reach \(n_c\), from the de-duplicated reach report, is the evidence. It enters the trust dial \(Z_c\) and nothing else.

• Reach-days \(\text{reach}_{ct}\), summing to \(R_c\), are the exposure base. They set the denominator of the rate we model and they spend that rate back out, day by day.

A person reached again tomorrow adds one reach-day but no new evidence about the campaign. The robustness checks at the end put this choice to the test by swapping the two.

The number we trust

The last input is the MMM’s answer for the channel as a whole.

social_daily = data.social.stack(sample=("chain", "draw"))
social_hdi = az.hdi(social_daily, dim="sample")

fig, ax = plt.subplots(figsize=(12, 7), layout="constrained")
ax.fill_between(
    data.mmm["date"],
    social_hdi.sel(ci_bound="lower").to_numpy(),
    social_hdi.sel(ci_bound="upper").to_numpy(),
    color="C0",
    alpha=0.25,
    label=f"{CI_PROB:.0%} HDI",
)
ax.plot(
    data.mmm["date"],
    social_daily.median("sample").to_numpy(),
    color="C0",
    lw=1.6,
    label="posterior median",
)
ax.set(
    title=f"Paid Social contribution from the MMM (total {channel_total / 1e3:,.0f}k)",
    xlabel="date",
    ylabel="daily contribution (outcome units)",
)
ax.xaxis.set_major_locator(mdates.MonthLocator())
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b"))
ax.legend(loc="upper left", fontsize=12);

Fig. What the MMM hands us: Paid Social’s daily incremental contribution with its 94% credible interval. That interval is a highest density interval (HDI), the narrowest range that holds 94% of the plausible values, and it is the uncertainty format used for every estimate in this notebook. The area under the curve, about 282k in the MMM’s outcome units, is the number we split across campaigns. Everything else in the notebook decides the proportions.

One more report worth pulling: campaign overlap

The same multi-report trick, asked of campaign groups, tells you how much your campaigns compete for the same people. For each pair, pull reach for A alone, for B alone, and for the combined A∪B group. The shared audience is then

\[\text{overlap}(A, B) = \text{Reach}(A) + \text{Reach}(B) - \text{Reach}(A \cup B).\]

Repeat over the pairs (and larger groups) you care about to fill an N by N matrix of shared people. That matrix is symmetric: the people in both A and B are the same people whichever way round you ask. Its diagonal is each campaign’s own audience size and the off-diagonals are shared users, zero where targeting is disjoint. In this synthetic example the matrix describes targeting universes, so a diagonal is a potential audience size rather than the realized cumulative unique reach used below.

Raw counts are hard to compare across campaigns of very different size, so the usual display divides each row by that row’s own audience. The cell in row A, column B then reads: “this share of A’s audience is also targeted by B”. That normalisation makes the matrix asymmetric, and the asymmetry is information. The same 22,990 shared people are 55% of spring_promo’s 42k audience but only 37% of brand_always_on’s 62k audience. The higher percentage always belongs to the smaller campaign, and it says that this campaign’s audience is largely nested inside the other’s, so it is the one more exposed to the competition. The figure shows both views side by side.

overlap = data.overlap.reindex(index=campaigns, columns=campaigns)
own_audience = pd.Series(np.diag(overlap), index=campaigns)
overlap_pct = 100 * overlap.div(own_audience, axis=0)

fig, (ax_people, ax_pct) = plt.subplots(
    ncols=2, figsize=(16, 7.5), layout="constrained"
)
sns.heatmap(
    overlap / 1e3,
    annot=True,
    fmt=".0f",
    cmap="Greens",
    linewidths=0.5,
    cbar_kws={"label": "thousands of people in both audiences"},
    ax=ax_people,
)
ax_people.set(
    title="Shared people (thousands): symmetric,\ndiagonal = own audience size",
    xlabel="",
    ylabel="",
)
sns.heatmap(
    overlap_pct,
    mask=np.eye(len(campaigns), dtype=bool),
    annot=True,
    fmt=".0f",
    cmap="Blues",
    linewidths=0.5,
    cbar_kws={"label": "% of the row's audience also targeted by the column"},
    ax=ax_pct,
)
ax_pct.set(
    title="Share of the row's audience the column also targets (%):\n"
    "directional, diagonal masked (always 100%)",
    xlabel="",
    ylabel="",
    yticklabels=[],
)
for ax in (ax_people, ax_pct):
    for tick_label in ax.get_xticklabels():
        tick_label.set(rotation=35, horizontalalignment="right")
fig.suptitle(
    "Audience overlap: how much campaigns compete for the same people", fontsize=15
);

Fig. Audience overlap from the pairwise reach queries. Left: shared people in thousands, one symmetric matrix, with each campaign’s own audience size on the diagonal. Right: the same matrix with each row divided by the row campaign’s own audience, so a cell is the share of the row’s audience that the column campaign also targets; the diagonal is masked because it is always 100%. Read the right panel by rows. 55% of spring_promo’s audience is also targeted by brand_always_on, but only 37% of brand_always_on’s audience is also targeted by spring_promo, because the same 23k shared people are a larger slice of the smaller audience. niche_interest is the most contested campaign: half of its 18k audience is also targeted by flash_sale_burst and 46% by video_product_launch. broad_prospecting is the umbrella: with 82k people it overlaps 40 to 44% of everyone else’s audience (read down its column), while any single other campaign covers at most 30% of its audience (read along its row). Pale cells are close to disjoint targeting. This is the competition structure inside your own account.

This notebook only needs the unique reach signal. The overlap matrix is not used below, but it comes from the same query pattern, it shows how much your campaigns are fighting over the same people, and it is the natural next input if you extend the method, so it is worth building at the same time.

Be clear about what that means: we measure the overlap here, we do not correct for it. The credibility model below treats each campaign’s unique reach as independent evidence, so two campaigns contesting the same 55% of an audience are trusted as if their reach were disjoint. Measured overlap is a diagnostic, not an adjustment; the closing section repeats this so it is not mistaken for one.

Step 1: each campaign’s own noisy story

To turn the dial we need two numbers per campaign: how much evidence it carries, its unique reach, and its own noisy story, a rough guess at how much revenue it drove per unit of exposure. Before we build them, here is every symbol the rest of the notebook uses. Each one is defined again where it first appears, so this table is a reference to come back to.

Notation at a glance

Symbol	Plain meaning	Where it comes from
\(c\), \(t\)	A campaign, a day	Indices
\(T\)	The MMM’s total social contribution over the window (posterior mean, about 282k)	The trusted number we divide up
\(\text{reach}_{ct}\)	Daily reach of campaign \(c\) on day \(t\)	Platform export
\(R_c = \sum_t \text{reach}_{ct}\)	Reach-days: the exposure the credit is spent over	Platform export, summed
\(n_c\)	Unique people reached over the flight: the evidence that sets trust	De-duplicated reach report
\(n^{\text{rel}}_c = n_c / \bar n\)	Unique reach relative to the average campaign (\(\bar n\) is the mean of the seven \(n_c\))	Derived; 1.0 means an average-sized campaign
\(\text{clickshare}_c\)	Campaign \(c\)’s share of the channel’s clicks	Platform export
\(y_c\)	The direct estimate: log credit per reach-day implied by the click split	Step 1
\(\theta_c\)	Campaign \(c\)’s true log rate per reach-day, the quantity we want	Step 2, estimated
\(\mu\)	The collective mean: the channel-average log rate	Step 2, estimated
\(\tau\)	How much true rates really differ across campaigns (a standard deviation on the log scale)	Step 2, estimated
\(s^2\)	Noise variance of a direct estimate for a campaign of average reach	Step 2, estimated
\(\psi_c = s^2 / n^{\text{rel}}_c\)	Noise variance of campaign \(c\)’s direct estimate: more reach, less noise	Assumption, Step 1
\(k = s^2 / \tau^2\)	The evidence threshold, in units of average reach; multiply by \(\bar n\) to get people	Step 2, derived
\(Z_c = \dfrac{n^{\text{rel}}_c}{n^{\text{rel}}_c + k} = \dfrac{\tau^2}{\tau^2 + \psi_c}\)	The trust dial: the share of campaign \(c\)’s own story we believe	Step 2, derived

For the own story we do something deliberately crude, exactly the kind of back-of-envelope split a platform dashboard invites. We take the MMM’s total social credit \(T\), hand it out in proportion to each campaign’s clicks, and express the result per reach-day:

\[y_c=\log\left(\frac{T \cdot \text{clickshare}_c}{R_c}\right),\qquad R_c=\sum_t \text{reach}_{ct}.\]

In words: \(y_c\) is the credit per reach-day that a click-share split would imply, on the log scale. We call it the direct estimate, because it comes straight from the campaign’s own numbers with no help from the other campaigns. It is noisy and non-incremental on purpose: this is the shaky signal we want to discipline, not trust. Two details are worth stating plainly.

• The denominator is reach-days, not people. That is the unit the daily allocation later spends, and matching the two is what makes the credit add up. Divide by unique reach instead and full trust in the click signal no longer recovers the click split even at the flight-total level, because a campaign’s frequency leaks into its share. (Step 3 spells out the exact sense in which full trust corresponds to the click split.)

• We work with the log of the rate, so “twice as effective” is the same size step whether a campaign is weak or strong. All the “log units” quoted below are steps on this scale; a step of 0.40 means a factor of 1.5.

One more thing to be precise about, since the framing is “one number you trust, seven shaky signals”. The MMM total \(T\) enters \(y_c\) as \(\log T + \log \text{clickshare}_c - \log R_c\), an identical shift for every campaign, which the collective mean \(\mu\) absorbs. So the MMM sets the level of the whole channel, and the ranking of the direct estimates comes from click share and exposure alone.

Clicks are a stand-in, and every click is priced the same

Two assumptions ride along with this construction, and neither is innocent.

• Clicks are a proxy, not the thing itself. Nothing in the machinery is attached to clicks: conversions, last-touch revenue, or a lift-informed signal slot into the same \(y_c\), and if you have a metric you believe tracks incrementality better, use it. What the model disciplines is the noise of whichever proxy you feed it; the proxy’s bias passes straight through (“it cures noise, not bias”, as the closing section puts it).

• A click is priced the same wherever it came from. The split hands out credit per click regardless of the campaign that earned it, yet the EDA showed median daily CTR differing by a factor of nearly four across campaigns. Two campaigns can deliver identical impressions with very different click counts simply because they were optimised for different things, say an awareness push against a click-optimised conversion campaign, and that difference is structural rather than noisy. Evidence turns the dial between a campaign’s own story and the collective average; it cannot repair a story that is told in the wrong currency.

The same precision applies to the MMM itself: it enters this construction only as the trusted total \(T\), not through its media variables or response curves. One structural piece of the MMM can be inherited cheaply: its adstock, the carryover that lets today’s exposure keep paying off for days afterwards. The adstock extension after the grading section does exactly that.

def naive_click_share(clicks: pd.Series) -> pd.Series:
    """Platform-style split of the channel total, proportional to clicks."""
    return clicks / clicks.sum()


def naive_log_rate(
    channel_total: float, click_share: pd.Series, exposure: pd.Series
) -> pd.Series:
    """Log contribution per unit of exposure implied by the naive click split."""
    return np.log(channel_total * click_share / exposure)


def direct_estimates(data: CampaignData) -> pd.DataFrame:
    """Per-campaign evidence, exposure, and the noisy direct estimate y_c."""
    campaigns = data.campaigns
    totals = campaign_totals(data.metrics, campaigns)
    channel_total = float(data.social.sum("date").mean())
    click_share = naive_click_share(totals["clicks"])
    return pd.DataFrame(
        {
            "unique_reach": final_unique_reach(data.cum_reach, campaigns),
            "reach_days": totals["reach_days"].astype(float),
            "clicks": totals["clicks"],
            "click_share": click_share,
            "naive_rate": channel_total * click_share / totals["reach_days"],
            "y_direct": naive_log_rate(
                channel_total, click_share, totals["reach_days"]
            ),
        }
    )


direct = direct_estimates(data)
y_direct = direct["y_direct"]
unique_reach = direct["unique_reach"]
direct.sort_values("unique_reach", ascending=False).round(4)

	unique_reach	reach_days	clicks	click_share	naive_rate	y_direct
brand_always_on	62386.0	1541137.0	40073	0.3903	0.0714	-2.6399
lookalike_conversions	42733.0	597676.0	16891	0.1645	0.0776	-2.5566
spring_promo	41719.0	438462.0	15985	0.1557	0.1001	-2.3020
flash_sale_burst	34802.0	139620.0	10084	0.0982	0.1982	-1.6183
video_product_launch	33427.0	519340.0	10628	0.1035	0.0562	-2.8794
broad_prospecting	24140.0	32921.0	1294	0.0126	0.1079	-2.2267
niche_interest	17987.0	254475.0	7709	0.0751	0.0831	-2.4872

The table gives seven plain numbers, one \(y_c\) per campaign. Nothing in the table says how reliable each of them is, and that is exactly what the trust dial needs to know. So we now make the one modelling assumption at the heart of the method, and the picture below draws it.

The assumption. Treat each direct estimate as a noisy reading of the campaign’s true log rate \(\theta_c\):

\[y_c = \theta_c + \varepsilon_c,\qquad \varepsilon_c \sim \text{Normal}\left(0, \sqrt{\psi_c}\right),\qquad \psi_c = \frac{s^2}{n^{\text{rel}}_c}.\]

The noise variance \(\psi_c\) is not measured from the data. It is assumed to shrink in proportion to how many people the campaign reached: a campaign of average size (\(n^{\text{rel}}_c = 1\)) has noise variance \(s^2\), a campaign that reached twice the average has half that variance, and so on. The noise standard deviation therefore falls like \(1/\sqrt{n_c}\). This is the same square-root rule that governs an opinion poll: an average over four times as many people has half the margin of error. Here the “poll” is the campaign’s click signal, and the people polled are the people it reached.

Where the bands come from. The single unknown in \(\psi_c\) is \(s\), the noise scale of an average campaign, and the model in Step 2 will estimate it. To draw a picture before fitting anything we need a stand-in, so we use a robust measure of how far the seven \(y_c\) spread around their median (the median absolute deviation, scaled so it estimates a standard deviation), which comes out at 0.275 log units (printed above the table). Each 94% band is then \(y_c \pm 1.88\sqrt{\psi_c}\), where 1.88 is the Normal quantile that leaves 3% in each tail. The table below lists the pieces for every campaign. Two examples: brand_always_on reached 1.70 times the average, so its noise standard deviation is \(0.275/\sqrt{1.70} \approx 0.21\) and its band half-width is \(1.88 \times 0.21 \approx 0.40\); niche_interest reached 0.49 times the average, so its standard deviation is \(0.275/\sqrt{0.49} \approx 0.39\) and its half-width about 0.74. The bands are the same recipe with the evidence dialled up or down.

def robust_scale(values: pd.Series) -> float:
    """Outlier-resistant estimate of the dispersion of the direct estimates."""
    return float(1.4826 * np.median(np.abs(values - np.median(values))))


# n_rel_c: unique reach relative to the average campaign (1.0 = average size)
relative_reach = unique_reach / unique_reach.mean()
# temporary stand-in for s, replaced by the fitted posterior in Step 2
s_plug_in = robust_scale(y_direct)
# psi_c = s^2 / n_rel_c, the assumed noise variance of each direct estimate
psi_plug_in = s_plug_in**2 / relative_reach
# 1.88: the Normal quantile that turns a standard deviation into a 94% half-width
z_score = float(norm.ppf(0.5 + CI_PROB / 2))
band_half_width = z_score * np.sqrt(psi_plug_in)
collective_mean = float(np.average(y_direct, weights=unique_reach))

print(f"plug-in s = {s_plug_in:.3f} log units, band multiplier = {z_score:.2f}")
band_table = pd.DataFrame(
    {
        "unique_reach": unique_reach.astype(int),
        "n_rel": relative_reach,
        "noise_sd": np.sqrt(psi_plug_in),
        "band_half_width": band_half_width,
        "y_direct": y_direct,
    }
).sort_values("unique_reach", ascending=False)
band_table.round(3)

plug-in s = 0.275 log units, band multiplier = 1.88

	unique_reach	n_rel	noise_sd	band_half_width	y_direct
brand_always_on	62386	1.698	0.211	0.396	-2.640
lookalike_conversions	42733	1.163	0.255	0.479	-2.557
spring_promo	41719	1.135	0.258	0.485	-2.302
flash_sale_burst	34802	0.947	0.282	0.531	-1.618
video_product_launch	33427	0.910	0.288	0.541	-2.879
broad_prospecting	24140	0.657	0.339	0.637	-2.227
niche_interest	17987	0.490	0.392	0.738	-2.487

ordered = unique_reach.sort_values().index
fig, ax = plt.subplots(figsize=(12, 7), layout="constrained")
ax.errorbar(
    unique_reach[ordered],
    y_direct[ordered],
    yerr=band_half_width[ordered],
    fmt="o",
    capsize=4,
    color="C0",
    label=f"direct estimate $y_c$ with assumed {CI_PROB:.0%} band",
)
ax.axhline(
    collective_mean,
    ls="--",
    color="black",
    label="collective mean (reach-weighted average of $y_c$)",
)
label_to_the_left = {"broad_prospecting", "flash_sale_burst"}
for campaign in ordered:
    to_the_right = campaign not in label_to_the_left
    ax.annotate(
        campaign,
        (unique_reach[campaign], y_direct[campaign]),
        xytext=(12, -3) if to_the_right else (-9, -3),
        textcoords="offset points",
        ha="left" if to_the_right else "right",
        fontsize=10,
    )
ax.text(
    0.02,
    0.04,
    "band half-width $= 1.88\\,\\sqrt{\\psi_c}$, with $\\psi_c = s^2 / n^{\\mathrm{rel}}_c$\n"
    f"plug-in $s$ = {s_plug_in:.3f} (robust spread of the seven $y_c$)",
    transform=ax.transAxes,
    fontsize=10,
    bbox={"boxstyle": "round", "facecolor": "white", "alpha": 0.9},
)
ax.margins(x=0.22, y=0.20)
ax.legend(loc="upper left")
ax.set(
    xscale="log",
    title=(
        "Direct estimates $y_c$ with assumed 94% bands "
        "(noise sd falls like $1/\\sqrt{n_c}$)"
    ),
    xlabel="unique reach $n_c$ (log scale)",
    ylabel="$y_c=\\log$ (naive credit per reach-day)",
);

Fig. Each campaign’s direct estimate \(y_c\) (dot) against its evidence \(n_c\) on a log axis, with the assumed 94% band. The band widths are not measured; they follow the assumption that noise falls like \(1/\sqrt{n_c}\), scaled by the plug-in \(s = 0.275\) from the table above. niche_interest carries the least evidence and gets the widest band; brand_always_on carries the most and gets the narrowest. Step 2 replaces the plug-in \(s\) with a value learned from the data.

Two campaigns anchor the range. flash_sale_burst sits highest: it ran 15 days and re-reached each person only four times, so its clicks are divided by very little exposure. video_product_launch sits lowest, with heavy delivery and the worst CTR in the account. niche_interest carries the least evidence of all, so its band is the widest on the chart even though its rate is close to the middle.

The dashed line is the collective mean: the reach-weighted average of the seven direct estimates, and a plug-in stand-in for the \(\mu\) the model will estimate. It is the level each campaign gets pulled toward, in proportion to how little evidence it carries.

Step 2: let the data set the dial

Now we let the data set the dial. Instead of decreeing “a campaign is only half-trusted until it reaches 50k people”, or handing everyone a flat Z = 0.5 by gut, we write the blend down as a small model and let it learn, alongside the channel average itself, the two numbers that set the dial, estimated from the seven campaigns together.

The model has two lines. The first is the measurement assumption from Step 1. The second says where the true rates come from:

\[y_c \sim \text{Normal}\left(\theta_c,\ \sqrt{\psi_c}\right),\qquad \psi_c = \frac{s^2}{n^{\text{rel}}_c}, \qquad\qquad \theta_c \sim \text{Normal}(\mu,\ \tau).\]

• \(\tau\) is how much campaigns genuinely differ from one another. It is the standard deviation of the true log rates \(\theta_c\) around the collective mean \(\mu\). If \(\tau\) is small, campaigns are basically alike, the group average is very informative and everyone gets pulled toward it hard.

• \(s^2\) is how noisy one campaign’s own signal is for a given amount of evidence. It is the noise variance of the direct estimate for a campaign of average size, and \(\psi_c = s^2 / n^{\text{rel}}_c\) scales it to each campaign. The noisier the signal, the more reach it takes to earn trust.

Their ratio is the evidence threshold, \(k = s^2/\tau^2\), measured in units of “average campaign reach” (multiply by the mean unique reach to get people), and each campaign’s dial follows automatically: \(Z_c = n^{\text{rel}}_c/(n^{\text{rel}}_c+k)\). Because we fit this as a Bayesian model we do not get a single \(Z\) per campaign, we get a distribution. Even the trust carries an error bar.

Statisticians call this partial pooling: every campaign borrows strength from the others, in proportion to how little evidence it has of its own, and the resulting pull toward the group average is called shrinkage. Both words appear in the figures below and mean exactly this.

You are already setting Z, you just do not call it that

Every day-to-day call about campaign performance is a Z in disguise, set by gut.

• Reading the dashboard at face value is Z = 1. “This campaign shows 8x ROAS, that one 2x, move the budget” trusts each number fully, no matter how few people it reached. It is why small campaigns swing between hero and flop week to week.

• Splitting the credit by delivery is Z = 0. You have thrown away each campaign’s own signal and handed out credit in proportion to reach alone. Splitting by spend, or evenly, is Z = 0 in spirit: the average rules and the campaign’s own story is ignored.

• “Wait for enough conversions before you trust the CPA” is choosing k, the evidence threshold where belief kicks in. Meta and Google do exactly this with their campaign learning phase.

• “That 12x looks too good, tiny audience, call it maybe 4x” is a blend. That is credibility done in your head, at an unspoken Z somewhere between 0 and 1.

You cannot avoid picking a Z. The only question is whether it is an arbitrary gut call that shifts by analyst and by mood, or one value estimated from how your campaigns actually behave, with its uncertainty attached. That is all the model below does.

The model is short, and it is the entire method. Everything before it was preparing its two inputs, and everything after it is presentation.

def build_credibility_model(
    y_direct: pd.Series,
    volume: pd.Series,
    *,
    sigma_s_scale: float = 0.5,
    log_k_sd: float = 0.75,
) -> pm.Model:
    """Buhlmann credibility model in which volume sets each campaign's trust dial."""
    campaigns = list(y_direct.index)
    volume_relative = (volume / volume.mean()).reindex(campaigns).to_numpy()
    with pm.Model(coords={"campaign": campaigns}) as model:
        volume_data = pm.Data("volume_relative", volume_relative, dims="campaign")
        mu = pm.Normal("mu", mu=float(y_direct.mean()), sigma=1.0)
        log_k = pm.Normal("log_k", 0.0, log_k_sd)
        k = pm.Deterministic("k", pm.math.exp(log_k))
        sigma_s = pm.HalfNormal("sigma_s", sigma_s_scale)
        s2 = pm.Deterministic("s2", sigma_s**2)
        tau = pm.Deterministic("tau", pm.math.sqrt(s2 / k))
        psi = pm.Deterministic("psi", s2 / volume_data, dims="campaign")
        pm.Deterministic("Z", volume_data / (volume_data + k), dims="campaign")
        theta_raw = pm.Normal("theta_raw", 0.0, 1.0, dims="campaign")
        theta = pm.Deterministic("theta", mu + tau * theta_raw, dims="campaign")
        pm.Normal(
            "y_obs",
            mu=theta,
            sigma=pm.math.sqrt(psi),
            observed=y_direct.to_numpy(),
            dims="campaign",
        )
    return model


credibility_model = build_credibility_model(y_direct, unique_reach)
pm.model_to_graphviz(credibility_model)

Fig. The model graph. The plate holds the seven campaigns, volume_relative is the data node carrying each campaign’s unique reach divided by the average unique reach (\(n^{\text{rel}}_c\)), and everything inside the plate is per-campaign. Only mu, log_k and sigma_s are global.

Reading it line by line, with the code names next to the symbols from the notation table:

• mu is \(\mu\), the channel-average log rate every campaign is judged against.

• log_k and k are \(\log k\) and \(k\), the evidence threshold in units of average reach. We sample its log so that the prior is easy to state.

• sigma_s is \(s\), the noise scale of one campaign’s direct estimate at average reach, and s2 is \(s^2\).

• tau is \(\tau = \sqrt{s^2/k}\), how much the true rates differ across campaigns. It is derived from \(s\) and \(k\) rather than sampled on its own (the dropdown explains why).

• psi is \(\psi_c = s^2 / n^{\text{rel}}_c\), the noise variance of each campaign’s direct estimate.

• Z is the trust dial \(Z_c = n^{\text{rel}}_c / (n^{\text{rel}}_c + k)\) that falls out of \(k\).

• theta is \(\theta_c\), each campaign’s trust-adjusted log rate: the blend of its own story and the group average. theta_raw is a standard Normal that is scaled by \(\tau\) and shifted by \(\mu\) to give \(\theta_c\), a standard reparametrisation that helps the sampler.

• y_obs is the observed direct estimate \(y_c\), with standard deviation \(\sqrt{\psi_c}\).

The model samples \(\log k\) directly, so the prior on the evidence threshold is a statement a media analyst can argue with: “a campaign of average size starts out about half believed, give or take”. The prior predictive check below shows what that implies before any data are used.

For the statistically minded: parametrisation, priors and a weak identification

We sample \(\log k\) and derive \(\tau=\sqrt{s^2/k}\), rather than sampling \(\tau\) and \(s^2\) separately, for two reasons. First, \(k\) is the actuarial quantity of interest, the amount of reach at which a campaign earns half its own trust, so putting a prior on it is meaningful. Second, the \((\tau, s^2)\) pair is only weakly identified from seven campaigns: sampling those two directly produced divergences at target_accept=0.95 while this version produces none.

Marginalising out \(\theta_c\), the model says

\[y_c \sim \text{Normal}\left(\mu,\; \sqrt{\tau^2+\psi_c}\right), \qquad \tau^2+\psi_c = s^2\left(\frac{1}{k}+\frac{1}{n^{\text{rel}}_c}\right),\]

so \(s^2\) and \(k\) enter only through that combination. What separates them is the spread of relative reach across campaigns, which here runs from 0.49 to 1.70. With seven campaigns and a 3.5-fold spread, that separation is weak, and the posterior for \(k\) stays wide. This is the honest reason the priors deserve a check rather than a shrug, and why the robustness checks sweep them.

Could more data sharpen it? The tempting route is time series: cut each flight into weeks and let dozens of weekly direct estimates \(y_{c,w}\) inform \(s^2\) and \(k\) instead of seven flight totals. We deliberately do not. The extra rows are far from independent evidence: unique reach is cumulative, so consecutive weekly evidence increments are strongly autocorrelated; the weekly click shares are still one composition per week; and within a campaign, whatever moves its share in one week (a creative, an audience, a promotion) tends to move it in the next. Modelling that dependence honestly costs more structure than seven campaigns can support, so this notebook keeps the flight-level model and spends its scepticism on the robustness checks instead.

On the sigma_s prior: a HalfNormal(0.5) implies simulated flights whose spread across campaigns has a median of 0.41, against the observed 0.40, with an interquartile range of 0.20 to 0.73. A tighter HalfNormal(0.4) would imply a median spread of 0.34 and a wider HalfNormal(0.7) 0.59, so the chosen prior covers the data comfortably from both sides without insisting on the observed value; the alternatives bracket it, and the robustness checks show none of them changes the conclusions.

One more disclosure: \(\mu\)’s prior is centered on the mean of the observed \(y_c\). With a prior standard deviation of 1.0 against a data spread of 0.40 the effect is mild, but it is the data informing a prior, not prior information.

with credibility_model:
    prior = pm.sample_prior_predictive(draws=1_000, random_seed=rng)

prior_z = prior.prior["Z"].stack(sample=("chain", "draw"))
prior_k_people = (
    prior.prior["k"].stack(sample=("chain", "draw")).to_numpy() * unique_reach.mean()
)
prior_spread = (
    prior.prior_predictive["y_obs"].stack(sample=("chain", "draw")).std("campaign")
)


def prior_spread_median(sigma_s_scale: float) -> float:
    """Median spread of the simulated flights implied by one `sigma_s` prior scale."""
    with build_credibility_model(y_direct, unique_reach, sigma_s_scale=sigma_s_scale):
        drawn = pm.sample_prior_predictive(draws=1_000, random_seed=seed)
    simulated = drawn.prior_predictive["y_obs"].stack(sample=("chain", "draw"))
    return float(simulated.std("campaign").median())


prior_spread_by_scale = pd.Series(
    {
        scale: prior_spread_median(scale)
        if scale != 0.5
        else float(prior_spread.median())
        for scale in (0.4, 0.5, 0.7)
    },
    name="prior median spread of y_c across campaigns",
).rename_axis("sigma_s prior scale")
print(f"observed spread of y_c across campaigns: {y_direct.std(ddof=1):.2f}")
display(prior_spread_by_scale.round(2))

prior_z_table = pd.DataFrame(
    {
        "unique_reach": unique_reach,
        "prior_median_Z": prior_z.median("sample").to_pandas(),
        "prior_Z_3%": prior_z.quantile(0.03, "sample").to_pandas(),
        "prior_Z_97%": prior_z.quantile(0.97, "sample").to_pandas(),
    }
).sort_values("unique_reach", ascending=False)
prior_z_table.round(3)

Sampling: [log_k, mu, sigma_s, theta_raw, y_obs]

Sampling: [log_k, mu, sigma_s, theta_raw, y_obs]

Sampling: [log_k, mu, sigma_s, theta_raw, y_obs]

observed spread of y_c across campaigns: 0.40

sigma_s prior scale
0.4    0.34
0.5    0.41
0.7    0.59
Name: prior median spread of y_c across campaigns, dtype: float64

	unique_reach	prior_median_Z	prior_Z_3%	prior_Z_97%
brand_always_on	62386.0	0.634	0.295	0.874
lookalike_conversions	42733.0	0.542	0.223	0.826
spring_promo	41719.0	0.537	0.219	0.823
flash_sale_burst	34802.0	0.491	0.189	0.795
video_product_launch	33427.0	0.481	0.183	0.788
broad_prospecting	24140.0	0.401	0.139	0.728
niche_interest	17987.0	0.333	0.108	0.667

fig, (ax_k, ax_z, ax_y) = plt.subplots(nrows=3, figsize=(10, 12), layout="constrained")

ax_k.hist(
    prior_k_people / 1e3,
    bins=np.geomspace(prior_k_people.min() / 1e3, prior_k_people.max() / 1e3, 70),
    color="C0",
    alpha=0.85,
)
# the prior median of k is exactly 1 in relative units, so it is the average reach
ax_k.axvline(
    unique_reach.mean() / 1e3,
    color="C2",
    ls="-",
    label=f"prior median = {unique_reach.mean() / 1e3:.0f}k",
)
ax_k.axvline(
    unique_reach.min() / 1e3, color="C1", ls="--", label="smallest campaign reach"
)
ax_k.axvline(
    unique_reach.max() / 1e3, color="C3", ls="--", label="largest campaign reach"
)
ax_k.set(
    xscale="log",
    xlabel="$k$ (thousand people, log scale)",
    title="Prior on the evidence threshold $k$, in people",
)
ax_k.legend(loc="upper left", fontsize=11)

z_order = unique_reach.sort_values().index.tolist()
prior_z_draws = pd.DataFrame(
    {campaign: prior_z.sel(campaign=campaign).to_numpy() for campaign in z_order}
).melt(var_name="campaign", value_name="Z")
sns.boxplot(
    data=prior_z_draws,
    x="Z",
    y="campaign",
    order=z_order,
    orient="h",
    width=0.6,
    showfliers=False,
    color="C0",
    ax=ax_z,
)
ax_z.scatter(
    [relative_reach[campaign] / (relative_reach[campaign] + 1) for campaign in z_order],
    np.arange(len(z_order)),
    color="C1",
    zorder=3,
    label="prior median",
)
ax_z.set(
    xlim=(0, 1),
    xlabel="$Z_c$",
    ylabel="",
    title="Prior on the trust dial $Z_c$",
)
# above the axes, so the legend cannot sit on top of a whisker
ax_z.legend(loc="lower right", bbox_to_anchor=(1.0, 1.01))

# bins spanning the full sampled range, so no part of the prior is clipped away
ax_y.hist(
    prior_spread.to_numpy(),
    bins=np.linspace(0, float(prior_spread.max()), 70),
    color="C0",
    alpha=0.85,
)
ax_y.axvline(
    float(prior_spread.median()),
    color="C2",
    ls="-",
    label=f"prior median = {float(prior_spread.median()):.2f}",
)
ax_y.axvline(
    float(y_direct.std(ddof=1)),
    color="C1",
    ls="--",
    label=f"observed spread of $y_c$ = {y_direct.std(ddof=1):.2f}",
)
ax_y.legend(loc="upper right")
ax_y.set(
    xlabel="spread of the seven simulated $y_c$ within a draw",
    title="Prior predictive spread of $y_c$ across campaigns",
);

Fig. Prior predictive checks, top to bottom. The top panel puts the evidence threshold \(k\) on a scale a media analyst can argue with. Its prior median is about 37k people, which is the average campaign’s unique reach in a field spanning 18k to 62k. In plain language: before seeing any data, a campaign of average size starts out about half believed, and the prior leaves plenty of room on both sides.

The middle panel shows what that implies for each campaign’s dial. The prior on \(Z_c\) deliberately does not span 0 to 1 for every campaign: because the prior median of \(k\) equals the average reach, the prior median of \(Z_c\) is \(n^{\text{rel}}_c/(n^{\text{rel}}_c+1)\), marked by the orange dots, which already encodes “bigger campaigns start out more trusted”. brand_always_on starts at a median of 0.63, niche_interest at 0.33, each with a wide range around it (see the table above). That is wide enough for the data to move the dial and narrow enough to rule out the two extremes a priori.

The bottom panel is the check on the noise prior: the simulated flights spread across campaigns by about as much as the real ones do (prior median 0.41 against an observed 0.40), with a wide range around it, so the prior covers the data without insisting on it.

Now we fit the model.

def fit_credibility(
    model: pm.Model, *, sample_seed: int = seed, progressbar: bool = True
) -> xr.DataTree:
    """Sample the credibility model with the settings used throughout the notebook."""
    with model:
        return pm.sample(
            draws=1_000,
            tune=1_000,
            chains=4,
            cores=4,
            target_accept=0.95,
            random_seed=sample_seed,
            progressbar=progressbar,
        )


idata = fit_credibility(credibility_model)

posterior = idata.posterior
divergences = int(idata.sample_stats["diverging"].sum())
max_rhat = float(az.summary(idata, var_names=["mu", "log_k", "sigma_s"])["r_hat"].max())
# k in units of the average campaign's unique reach, then converted to people
k_relative = float(posterior["k"].median())
k_people = k_relative * unique_reach.mean()
trust = posterior["Z"].mean(("chain", "draw")).to_pandas().reindex(campaigns)

print(
    f"{divergences} divergences, max R-hat {max_rhat:.3f}: the sampler is healthy.\n"
    f"Posterior median evidence threshold k = {k_relative:.2f} times the average "
    f"campaign's unique reach ({unique_reach.mean():,.0f} people), "
    f"that is {k_people:,.0f} people.\n"
    f"The fitted trust dial runs from {trust.min():.2f} ({trust.idxmin()}) "
    f"to {trust.max():.2f} ({trust.idxmax()})."
)
az.summary(idata, var_names=["mu", "k", "sigma_s", "tau"]).round(3)

NUTS[nutpie]: [mu, log_k, sigma_s, theta_raw]

0 divergences, max R-hat 1.000: the sampler is healthy.
Posterior median evidence threshold k = 1.01 times the average campaign's unique reach (36,742 people), that is 37,061 people.
The fitted trust dial runs from 0.34 (niche_interest) to 0.61 (brand_always_on).

	mean	sd	hdi94_lb	hdi94_ub	ess_bulk	ess_tail	r_hat	mcse_mean	mcse_sd
mu	-2.401	0.194	-2.8	-2	1785	1687	1.00	0.0049	0.0048
k	1.32	1.1	0.24	4.1	1719	1722	1.00	0.023	0.056
sigma_s	0.33	0.118	0.16	0.6	1513	1466	1.00	0.0029	0.0028
tau	0.332	0.133	0.16	0.64	3141	2350	1.00	0.0026	0.0035

Two things to read off. First, the sampler diagnostics: no divergences and \(\hat R\) of 1.00 mean the four chains agree, so the numbers below can be trusted as a faithful picture of the posterior. Second, the threshold. In the summary table k is in units of the average campaign’s reach, so a value of 1 means “the average campaign is believed exactly half-way”; the print statement converts it to people. It is reported as the posterior median rather than the mean, because \(k\) is sampled on the log scale and its posterior is right-skewed: the mean of 1.32 in the table is pulled up by a long right tail. Note what the fitted range of \(Z\) does and does not say: even the best-measured campaign in this account is only slightly more than half believed, and the smallest sits near a third.

azp.plot_trace_dist(idata, var_names=["mu", "k", "sigma_s"]);

Fig. Marginal posteriors and traces for the three sampled quantities. The chains mix, and k is visibly skewed and wide, which is the weak identification described in the dropdown above showing up in the output rather than hiding in it.

The fitted dial

This is the picture from the introduction, now with the data in it. The curve is the fitted trust dial, the band around it is the uncertainty in the threshold \(k\), and each dot is a real campaign placed at its unique reach, with a whisker for the uncertainty in its own \(Z_c\).

k_draws_people = posterior["k"].to_numpy().ravel() * unique_reach.mean()
z_curves = xr.DataArray(
    people_grid[None, :] / (people_grid[None, :] + k_draws_people[:, None]),
    dims=("sample", "n"),
)
z_curve_hdi = az.hdi(z_curves, dim="sample")
z_posterior = posterior["Z"].stack(sample=("chain", "draw"))
z_bounds = az.hdi(z_posterior, dim="sample")

fig, ax = plt.subplots(figsize=(10, 6.5), layout="constrained")
ax.fill_between(
    people_grid,
    z_curve_hdi.sel(ci_bound="lower").to_numpy(),
    z_curve_hdi.sel(ci_bound="upper").to_numpy(),
    color="C0",
    alpha=0.18,
    label=f"{CI_PROB:.0%} HDI of the curve (uncertainty in $k$)",
)
ax.plot(
    people_grid,
    z_curves.median("sample").to_numpy(),
    color="C0",
    lw=2,
    label="posterior median dial",
)

for campaign in campaigns:
    ax.errorbar(
        unique_reach[campaign],
        trust[campaign],
        yerr=[
            [
                trust[campaign]
                - float(z_bounds.sel(campaign=campaign, ci_bound="lower"))
            ],
            [
                float(z_bounds.sel(campaign=campaign, ci_bound="upper"))
                - trust[campaign]
            ],
        ],
        fmt="o",
        ms=8,
        capsize=3,
        elinewidth=1.2,
        alpha=0.85,
        color=colors[campaign],
        label=campaign,
        zorder=4,
    )
ax.axvline(k_people, ls=":", color="gray")
ax.annotate(
    f"fitted $k$ = {k_people / 1e3:.0f}k people:\na campaign this big is\nbelieved 50/50",
    (k_people, 0.5),
    xytext=(110, -150),
    textcoords="offset points",
    ha="left",
    fontsize=15,
    arrowprops={"arrowstyle": "->", "color": "gray"},
)
ax.legend(loc="upper left", fontsize=11)
ax.set(
    xscale="log",
    xlim=(5e3, 5e5),
    ylim=(0, 1),
    xlabel="unique people reached, $n_c$ (log scale)",
    ylabel="trust dial $Z_c = n_c / (n_c + k)$",
    title="The fitted dial: how much each campaign's own story is believed",
);

Fig. The fitted trust dial with the seven campaigns on it. Reading it left to right: niche_interest at 18k people is believed about a third, brand_always_on at 62k about 61%, and the whole account sits on the steep part of the curve, so extra reach still buys trust quickly. The band is wide, because seven campaigns pin down the threshold \(k\) only loosely; that is honest, and the robustness checks measure how much it matters. Compared with the introduction’s picture, the data have chosen the middle curve: a threshold of about 37k people.

Watch the rates move toward the collective centre

This is the payoff picture. Each orange dot is a campaign’s naive self-report and each arrow shows where the model moved it after weighing the evidence. The number under each campaign is its trust dial Z, and campaigns are sorted by unique reach with the smallest on the left.

The size of each move is not arbitrary. Credibility says

\[\theta_c - y_c = (1-Z_c)\,(\mu - y_c),\]

so a campaign moves by its distance from the collective mean, discounted by how much trust it has earned. That predicts the two things you see. flash_sale_burst starts furthest above the collective rate, because its clicks are divided by only four reach-days per person, and takes the largest step down. video_product_launch starts furthest below, with the account’s worst CTR against heavy delivery, and takes the largest step up. niche_interest has the lowest trust dial of all but barely moves, because it started close to the middle.

That is the discipline this method buys: a loud dashboard number is not a licence to scale a campaign, and a quiet one is not proof of failure. The adjustment is earned by evidence, and it protects a big steady campaign from being punished for one unlucky week just as much as it reins in a small loud one.

collective_rate = float(np.exp(posterior["mu"].mean()))
theta_mean = posterior["theta"].mean(("chain", "draw")).to_pandas().reindex(campaigns)
blend = trust * y_direct + (1 - trust) * float(posterior["mu"].mean())
blend_gap = float(np.abs(theta_mean - blend).max())

shrinkage = pd.DataFrame(
    {
        "unique_reach": unique_reach.astype(int),
        "reach_days": direct["reach_days"].astype(int),
        "Z": trust,
        "naive_rate": direct["naive_rate"],
        "credible_rate": np.exp(theta_mean),
        "collective_rate": collective_rate,
    }
).sort_values("unique_reach")

print(
    f"Collective rate exp(mu) = {collective_rate:.4f} per reach-day. The posterior mean "
    f"theta reproduces the blend Z * y + (1 - Z) * mu to within {blend_gap:.3f} log units."
)
shrinkage.round(4)

Collective rate exp(mu) = 0.0906 per reach-day. The posterior mean theta reproduces the blend Z * y + (1 - Z) * mu to within 0.004 log units.

	unique_reach	reach_days	Z	naive_rate	credible_rate	collective_rate
niche_interest	17987	254475	0.3434	0.0831	0.0882	0.0906
broad_prospecting	24140	32921	0.4051	0.1079	0.0969	0.0906
video_product_launch	33427	519340	0.4766	0.0562	0.0719	0.0906
flash_sale_burst	34802	139620	0.4855	0.1982	0.1323	0.0906
spring_promo	41719	438462	0.5259	0.1001	0.0958	0.0906
lookalike_conversions	42733	597676	0.5312	0.0776	0.0834	0.0906
brand_always_on	62386	1541137	0.6138	0.0714	0.0784	0.0906

positions = np.arange(len(shrinkage))
fig, ax = plt.subplots(figsize=(10, 6), layout="constrained")
ax.axhline(
    np.log(collective_rate), color="black", ls="--", label=r"collective mean $\mu$"
)
for position, (_, row) in enumerate(shrinkage.iterrows()):
    ax.annotate(
        "",
        xy=(position, np.log(row["credible_rate"])),
        xytext=(position, np.log(row["naive_rate"])),
        arrowprops={"arrowstyle": "->", "color": "C0", "lw": 1.5},
    )
ax.scatter(
    positions,
    np.log(shrinkage["naive_rate"]),
    color="C1",
    zorder=3,
    label="naive (own click signal)",
)
ax.scatter(
    positions,
    np.log(shrinkage["credible_rate"]),
    color="C0",
    zorder=3,
    label="credible (shrunk)",
)
ax.legend(loc="upper left")
ax.set(
    xticks=positions,
    xticklabels=[
        f"{campaign}\n$n$={reach:,}\n$Z$={z:.2f}"
        for campaign, reach, z in zip(
            shrinkage.index,
            shrinkage["unique_reach"],
            shrinkage["Z"],
            strict=True,
        )
    ],
    ylabel="log rate per reach-day",
    title="Credibility shrinkage, campaigns sorted by unique reach (left = smallest)",
)
ax.tick_params(axis="x", labelsize=9)
ax.margins(y=0.16);

Fig. Naive self-report (orange) to trust-adjusted estimate (blue), as log rates per reach-day. Each campaign moves toward the collective mean by a fraction \(1-Z_c\) of its distance from it, so the largest steps belong to the campaigns that combine a low trust dial with an extreme starting point. The collective rate is 0.0906 per reach-day, and the table above confirms that the posterior means really are the blend \(Z_c y_c + (1-Z_c)\mu\), to within 0.004 log units.

Step 3: hand the credit back, day by day

The trust-adjusted rates say how much each campaign earns per reach-day. To turn that into daily credit we score each campaign on each day by its credible rate times its reach that day, and hand out the MMM’s channel curve in those proportions. Concretely: if campaign A is twice as effective per reach-day, but campaign B delivered three times as many reach-days on Tuesday, B gets 1.5 times A’s share of Tuesday’s credit. In symbols, the weight of campaign \(c\) on day \(t\) and the credit it receives are

\[w_{ct} = \frac{e^{\theta_c}\,\text{reach}_{ct}}{\sum_{c'} e^{\theta_{c'}}\,\text{reach}_{c't}},\qquad \text{credit}_{ct} = w_{ct}\cdot \text{social}_t,\]

where \(\text{social}_t\) is the MMM’s social contribution on day \(t\) and \(e^{\theta_c}\) is campaign \(c\)’s trust-adjusted rate per reach-day.

The weights add to one every single day, so the campaign pieces always sum back exactly to the MMM’s channel total. We are not inventing a second model of the world, we are dividing the one number we trust.

Both sources of uncertainty travel through. The MMM posterior holds 4 chains of 100 draws, so we pair all 400 of them with 400 draws from the credibility posterior and carry the pairs through the allocation.

Because the rate is per reach-day, the two ends of the dial are quantities you already recognise:

• at \(Z=0\) every campaign is judged by the collective rate alone, so each day’s credit follows delivery, giving the daily-reach split;

• at \(Z=1\) every campaign is believed completely, so each day’s credit follows the click split spread over each campaign’s own daily delivery profile: the day-\(t\) score is \(\text{clickshare}_c \cdot \text{reach}_{ct}/R_c\), the click share times the fraction of the campaign’s reach-days that fell on day \(t\).

One subtlety at the \(Z=1\) end is worth being precise about, because the sweep section quotes both numbers. The reach-days cancel in the flight-total score, \(\sum_t \text{score}_{ct} \propto \text{clickshare}_c\), but the credit is handed out after a per-day renormalisation: each day’s scores are rescaled to sum to one against whichever campaigns delivered that day, and then divide that day’s MMM credit. So \(Z=1\) reproduces the click split exactly only in the special case where every campaign has the same normalised daily delivery profile, so that the profiles cancel day by day; otherwise it is a daily-renormalised cousin of it. That is why, in the grading section, the sweep’s \(Z=1\) setting and the naive click split score differently against the truth (an MAE of 0.0373 against 0.0432).

So the dial genuinely interpolates between (daily versions of) the two splits an analyst would reach for first, and the fitted answer sits in between. The common-dial sweep in the next section makes this visible.

One structural caveat: the weights are undefined on a day when the whole channel delivered nothing, since the scores would all be zero. Every date in this window has delivery somewhere. On a dark day you would need a convention, and the honest one is that a day with no delivery earns no credit.

def daily_matrix(
    metrics: pd.DataFrame, campaigns: list[str], column: str
) -> xr.DataArray:
    """One daily column of the platform export as a (date, campaign) array."""
    wide = metrics.pivot(index="date", columns="campaign", values=column).reindex(
        columns=campaigns
    )
    return xr.DataArray(
        wide.to_numpy(dtype=float),
        dims=("date", "campaign"),
        coords={"date": wide.index, "campaign": campaigns},
    )


def daily_reach(metrics: pd.DataFrame, campaigns: list[str]) -> xr.DataArray:
    """Daily duplicated reach per campaign, the exposure the credit is spent over."""
    return daily_matrix(metrics, campaigns, "reach")


def paired_draws(
    social: xr.DataArray,
    theta: xr.DataArray,
    *,
    pair_seed: int = seed,
) -> tuple[xr.DataArray, xr.DataArray]:
    """Pair MMM and credibility draws so both uncertainties propagate."""
    flat_social = social.stack(sample=("chain", "draw")).transpose("sample", "date")
    flat_theta = theta.stack(sample=("chain", "draw")).transpose("sample", "campaign")
    n_draws = flat_social.sizes["sample"]
    available = flat_theta.sizes["sample"]
    picked = np.random.default_rng(pair_seed).choice(
        available, size=n_draws, replace=available < n_draws
    )
    sample = np.arange(n_draws)
    return (
        flat_social.drop_vars(["sample", "chain", "draw"]).assign_coords(sample=sample),
        flat_theta.isel(sample=picked)
        .drop_vars(["sample", "chain", "draw"])
        .assign_coords(sample=sample),
    )


def allocation_weights(
    theta_draws: xr.DataArray, exposure: xr.DataArray
) -> xr.DataArray:
    """Daily weights proportional to the credible rate times daily exposure."""
    score = np.exp(theta_draws - theta_draws.max("campaign")) * exposure
    return score / score.sum("campaign")


def allocate(social_draws: xr.DataArray, weights: xr.DataArray) -> xr.DataArray:
    """Spend each MMM channel draw over campaigns using the daily weights."""
    return (weights * social_draws).rename("campaign_contribution")


def flight_share(contribution: xr.DataArray) -> xr.DataArray:
    """Flight-total share of the channel per campaign, one value per draw."""
    flight = contribution.sum("date")
    return flight / flight.sum("campaign")


reach = daily_reach(data.metrics, campaigns)
social_draws, theta_draws = paired_draws(data.social, posterior["theta"])
contribution = allocate(social_draws, allocation_weights(theta_draws, reach))
share = flight_share(contribution).transpose("sample", "campaign")

delivery_split = (reach.sum("date").to_pandas() / float(reach.sum())).reindex(campaigns)
click_split = direct["click_share"]

reconciliation = float(np.abs(contribution.sum("campaign") - social_draws).max())
print(
    f"Largest per-day gap between the campaign pieces and the MMM draw, across all "
    f"{share.sizes['sample']} paired draws: {reconciliation:.1e}. The split is exact."
)

Largest per-day gap between the campaign pieces and the MMM draw, across all 400 paired draws: 1.8e-12. The split is exact.

Here is the split itself: every day’s MMM credit divided among the campaigns, stacking back to the channel curve.

daily_mean = contribution.mean("sample")

fig, ax = plt.subplots(figsize=(12, 6.5), layout="constrained")
ax.stackplot(
    reach.date,
    [daily_mean.sel(campaign=campaign).to_numpy() for campaign in campaigns],
    labels=campaigns,
    colors=[colors[campaign] for campaign in campaigns],
    alpha=0.85,
)
ax.plot(
    reach.date,
    social_draws.mean("sample").to_numpy(),
    color="black",
    lw=1.4,
    label="social total (MMM)",
)
ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.1), ncol=4)
ax.set(
    title="Posterior-mean daily credit per campaign, stacking to the MMM social curve",
    xlabel="date",
    ylabel="contribution (outcome units)",
)
ax.xaxis.set_major_locator(mdates.MonthLocator())
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b"));

Fig. Daily credit per campaign (stacked areas, posterior mean) against the MMM’s social curve (black line). The pieces stack to the curve by construction. Within a campaign the day-to-day shape is inherited from its daily reach, which is why the broad_prospecting slice appears from nothing in April and the flash_sale_burst slice is a 15-day wedge in February. The grading section asks whether that inherited shape is right, and the adstock extension improves it.

How much of the click signal survives?

The split started from clicks, so a natural check is how the credited share compares with the click share. The scatter below puts each campaign’s share of the channel’s clicks on the x axis and its share of the credited contribution on the y axis (the posterior-mean daily credit summed over the window, as a share of the channel). A campaign on the dashed diagonal is credited exactly its click share. The hollow diamond above or below each dot is the same campaign’s share of reach-days, the delivery split, which is what the dial pulls toward.

credited_share = (
    (daily_mean.sum("date") / daily_mean.sum()).to_pandas().reindex(campaigns)
)
proxy_check = pd.DataFrame(
    {
        "click_share": click_split,
        "credited_share": credited_share,
        "reach_days_share": delivery_split,
        "credited_over_click": credited_share / click_split,
        "ctr_%": summary["ctr_%"].reindex(campaigns),
        "Z": trust,
    }
).sort_values("click_share", ascending=False)

# keep each label clear of the neighbouring campaign's markers
label_offsets = {
    "video_product_launch": (-12, -4),
    "spring_promo": (12, -14),
    "lookalike_conversions": (12, 4),
    "flash_sale_burst": (-12, 6),
}
fig, ax = plt.subplots(figsize=(9, 8), layout="constrained")
lower, upper = 0.007, 0.6
ax.plot(
    (lower, upper),
    (lower, upper),
    color="black",
    ls="--",
    lw=1,
    label="credited = clicks",
)
for campaign in campaigns:
    ax.plot(
        [click_split[campaign]] * 2,
        [delivery_split[campaign], credited_share[campaign]],
        color=colors[campaign],
        lw=1,
        alpha=0.6,
    )
    ax.scatter(
        click_split[campaign],
        delivery_split[campaign],
        marker="D",
        facecolors="none",
        edgecolors=colors[campaign],
        s=60,
        zorder=3,
    )
    ax.scatter(
        click_split[campaign],
        credited_share[campaign],
        color=colors[campaign],
        s=90,
        zorder=4,
    )
    offset = label_offsets.get(campaign, (10, -12))
    ax.annotate(
        campaign,
        (click_split[campaign], credited_share[campaign]),
        xytext=offset,
        textcoords="offset points",
        ha="right" if offset[0] < 0 else "left",
        fontsize=9,
    )
ax.scatter([], [], color="gray", s=90, label="credited share (credibility split)")
ax.scatter(
    [],
    [],
    marker="D",
    facecolors="none",
    edgecolors="gray",
    s=60,
    label="reach-days share (delivery)",
)
ax.legend(loc="lower right")
ax.set(
    xscale="log",
    yscale="log",
    xlim=(lower, upper),
    ylim=(lower, upper),
    xlabel="share of the channel's clicks (log scale)",
    ylabel="share of the channel's credited contribution (log scale)",
    title="Clicks are only partly honoured: credited share against click share",
)
proxy_check.round(3)

	click_share	credited_share	reach_days_share	credited_over_click	ctr_%	Z
brand_always_on	0.390	0.414	0.437	1.062	1.357	0.614
lookalike_conversions	0.165	0.170	0.170	1.032	2.046	0.531
spring_promo	0.156	0.140	0.124	0.897	2.284	0.526
video_product_launch	0.104	0.133	0.147	1.281	0.998	0.477
flash_sale_burst	0.098	0.056	0.040	0.569	2.723	0.486
niche_interest	0.075	0.076	0.072	1.016	2.173	0.343
broad_prospecting	0.013	0.011	0.009	0.910	3.807	0.405

Fig. Credited share against click share, both on log axes so the small campaigns stay visible. Three campaigns sit on the diagonal: lookalike_conversions, niche_interest and broad_prospecting have a click share that already matches their delivery share, so the dial has nothing to arbitrate and they are credited what their clicks suggest. The other four are where credibility does its work. flash_sale_burst earned 9.8% of the clicks but is credited 5.6% of the channel: its clicks are far above its delivery (4.0% of reach-days), and with a trust dial near 0.49 the model keeps only about half of the excess. video_product_launch moves the other way, from 10.4% of clicks to 13.3% of credit: it has the account’s lowest CTR against heavy delivery (14.7% of reach-days), so the group average tops it up. brand_always_on (39% to 41%) and spring_promo (15.6% to 14.0%) are milder versions of the same two stories.

Two things are worth noticing. First, the credited share lands between the click share and the delivery share, or very close to it: the dial interpolates, it does not invent. Exactly between is not guaranteed, because the weights are renormalised every day rather than once per flight, and niche_interest sits a hair above both. Second, the only change in ranking relative to the click split is that niche_interest moves above flash_sale_burst. Whether these moves were the right ones is a question for the grading section, which is allowed to look at the truth.

And here is what you bring to the meeting: each campaign’s credited contribution over the window, with its uncertainty, its share of the channel, and the trust dial that produced it.

flight_total = contribution.sum("date")
flight_bounds = az.hdi(flight_total, dim="sample")
meeting_table = (
    pd.DataFrame(
        {
            "credited_contribution": flight_total.mean("sample").to_pandas(),
            "hdi_low": flight_bounds.sel(ci_bound="lower").to_pandas(),
            "hdi_high": flight_bounds.sel(ci_bound="upper").to_pandas(),
            "share_%": 100 * share.mean("sample").to_pandas(),
            "Z": trust,
        }
    )
    .reindex(campaigns)
    .sort_values("credited_contribution", ascending=False)
)
meeting_table.round(
    {"credited_contribution": 0, "hdi_low": 0, "hdi_high": 0, "share_%": 1, "Z": 2}
)

	credited_contribution	hdi_low	hdi_high	share_%	Z
campaign
brand_always_on	116765.0	79341.0	158825.0	41.4	0.61
lookalike_conversions	47835.0	28134.0	69718.0	17.0	0.53
spring_promo	39345.0	22720.0	55161.0	14.0	0.53
video_product_launch	37360.0	22575.0	58672.0	13.3	0.48
niche_interest	21505.0	11559.0	34371.0	7.6	0.34
flash_sale_burst	15735.0	9528.0	22948.0	5.6	0.49
broad_prospecting	3233.0	1522.0	5063.0	1.1	0.41

The credited_contribution column sums to the MMM’s channel total, and its interval carries both the MMM’s uncertainty about the channel and the credibility model’s uncertainty about the split. On real data this table is the deliverable. Here we can go one step further and check it against the truth.

Did it work?

Because this is simulated data we can look at the truth we hid earlier and grade the result. One thing first: the hidden truth totals 540,350 in outcome units against the MMM’s 281,777 for this channel, a factor of 1.92, because in the simulation the MMM under-credits the channel as a whole. That gap is the MMM’s business, not ours: we divide what the MMM hands us, so we grade shares of the channel, not absolute levels.

We ask four questions of every method, and here is what each score means in plain terms:

• Ranking (Spearman correlation). Did we get the order of campaigns right? 1 is a perfect order, 0 is no relationship.

• Share error (MAE). How far off is each campaign’s share of the channel, on average? An MAE of 0.03 means the typical campaign is off by 3 percentage points of the channel’s credit.

• Probabilistic quality (CRPS). Our method returns a whole distribution while every simple heuristic returns one number. The continuous ranked probability score grades both on the same footing, in share units: for a single number it reduces exactly to the absolute error, and for a distribution it rewards being centred and being calibrated. It is a proper score, so padding the interval to look safe is punished.

• Coverage. How many of the seven true shares fall inside our 94% interval? Only the credibility split has an interval, so only it gets a coverage count. With seven campaigns this count is a coarse instrument, moving by a campaign or so between runs, so we read it as a sanity check rather than a target.

true_contribution = data.truth.sum(axis=0).reindex(campaigns)
true_share = true_contribution / true_contribution.sum()
share_bounds = az.hdi(share, dim="sample")

results = pd.DataFrame(
    {
        "true_share": true_share,
        "credibility_share": share.mean("sample").to_pandas(),
        "hdi_low": share_bounds.sel(ci_bound="lower").to_pandas(),
        "hdi_high": share_bounds.sel(ci_bound="upper").to_pandas(),
        "naive_clicks_share": click_split,
        "crps": per_observation_crps(true_share.to_numpy(), share.to_numpy()),
    }
).sort_values("true_share", ascending=False)
results.round(4)

	true_share	credibility_share	hdi_low	hdi_high	naive_clicks_share	crps
brand_always_on	0.4310	0.4140	0.3054	0.5078	0.3903	0.0130
lookalike_conversions	0.2406	0.1696	0.1070	0.2241	0.1645	0.0546
niche_interest	0.1096	0.0763	0.0467	0.1239	0.0751	0.0245
spring_promo	0.0878	0.1398	0.0904	0.1837	0.1557	0.0382
video_product_launch	0.0782	0.1329	0.0775	0.1920	0.1035	0.0382
flash_sale_burst	0.0419	0.0560	0.0338	0.0784	0.0982	0.0086
broad_prospecting	0.0110	0.0115	0.0054	0.0168	0.0126	0.0007

def score_shares(true_share: pd.Series, draws: np.ndarray) -> dict[str, float]:
    """Grade a share forecast given draws with shape (sample, campaign)."""
    truth = true_share.to_numpy()
    point = draws.mean(axis=0)
    scores = {
        "spearman": float(spearmanr(point, truth).correlation),
        "mae": float(np.abs(point - truth).mean()),
        "crps": float(crps(truth, draws)),
        "coverage": float("nan"),
    }
    if draws.shape[0] > 1:
        bounds = az.hdi(xr.DataArray(draws, dims=("sample", "campaign")), dim="sample")
        inside = (truth >= bounds.sel(ci_bound="lower").to_numpy()) & (
            truth <= bounds.sel(ci_bound="upper").to_numpy()
        )
        scores["coverage"] = float(inside.sum())
    return scores


totals = campaign_totals(data.metrics, campaigns)
baselines = {
    "daily-reach share": delivery_split,
    "spend share": totals["spend"] / totals["spend"].sum(),
    "naive clicks share": click_split,
    "impressions share": totals["impressions"] / totals["impressions"].sum(),
    "unique-reach share": unique_reach / unique_reach.sum(),
}

scoreboard = pd.DataFrame(
    [
        {"method": "credibility split", **score_shares(true_share, share.to_numpy())},
        *(
            {
                "method": name,
                **score_shares(
                    true_share, baseline.reindex(campaigns).to_numpy()[None, :]
                ),
            }
            for name, baseline in baselines.items()
        ),
    ]
).set_index("method")
scoreboard.assign(coverage=lambda df: df["coverage"].map("{:.0f} of 7".format)).replace(
    "nan of 7", ""
).round(4)

	spearman	mae	crps	coverage
method
credibility split	0.8929	0.0347	0.0254	5 of 7
daily-reach share	0.8571	0.0321	0.0321
spend share	0.8571	0.0322	0.0322
naive clicks share	0.7857	0.0432	0.0432
impressions share	0.6429	0.0481	0.0481
unique-reach share	0.6071	0.0864	0.0864

Five statements, each of which follows from the tables above.

1. Against the signal it disciplines, it improves everything. The naive clicks split scores 0.786 on Spearman and 0.0432 on both MAE and CRPS, which coincide for a point forecast. The credibility split scores 0.893, 0.0347 and 0.0254. In plain terms, the typical campaign is now off by 3.5 percentage points of the channel’s credit instead of 4.3, and the order is closer to right.

2. As a probabilistic forecast it is the best of everything tested, at a CRPS of 0.0254 against 0.0321 for the closest heuristic, the daily-reach share, because CRPS rewards calibration as well as centring.

3. On point accuracy alone the daily-reach share is still slightly better, 0.0321 against 0.0347. That is worth stating plainly rather than burying. Clicks are the weakest signal in this account, and a method that disciplines a signal inherits its slant.

4. 5 of the 7 true shares fall inside the 94% HDI. The interval carries parameter uncertainty only, so the real uncertainty is at least this wide: nothing in the model represents the gap between clicks and incrementality.

5. The dial moved six of the seven campaigns toward the truth. Compare the naive_clicks_share and credibility_share columns with true_share in the first table: every move goes the right way except video_product_launch, where the click share (0.104) was closer to the truth (0.078) than the credited share (0.133). The pull toward delivery helped where clicks over-claimed (flash_sale_burst, spring_promo) and hurt where a low click-through rate was honest. The model cannot tell those two cases apart, because nothing in it knows what a click is worth. That is the “cures noise, not bias” limit made concrete.

Keep the sample size in view. With 7 campaigns, a Spearman of 0.786 corresponds to a sum of squared rank differences of 12 and 0.893 to 6, so the improvement is 3 adjacent swaps’ worth of reordering, and MAE differences of 0.0026 are not evidence of much in either direction. Coverage is the most fragile score of the four: with seven campaigns the count moves by a campaign across seeds, as the robustness checks confirm, so it is a sanity check and not a quantity to tune against.

fig, ax = plt.subplots(figsize=(9, 8), layout="constrained")
# one square window wide enough for the widest HDI, so no interval is clipped
upper = 1.06 * max(
    results[["true_share", "credibility_share", "hdi_high"]].to_numpy().max(),
    click_split.max(),
    baselines["spend share"].max(),
)
limits = (-0.02, upper)
ax.plot((0, upper), (0, upper), color="black", ls="--", lw=1.0, label="perfect")
ax.errorbar(
    results["true_share"],
    results["credibility_share"],
    yerr=[
        results["credibility_share"] - results["hdi_low"],
        results["hdi_high"] - results["credibility_share"],
    ],
    fmt="o",
    capsize=4,
    color="C0",
    label=f"credibility ({CI_PROB:.0%} HDI)",
)
ax.scatter(
    true_share,
    click_split,
    marker="s",
    color="C1",
    zorder=3,
    label="naive clicks split",
)
ax.scatter(
    true_share,
    baselines["spend share"],
    marker="^",
    color="C2",
    zorder=3,
    label="spend share",
)
label_offsets = {"video_product_launch": (-9, -6), "brand_always_on": (-9, 6)}
for campaign in campaigns:
    offset = label_offsets.get(campaign, (9, -3))
    ax.annotate(
        campaign,
        (true_share[campaign], results.loc[campaign, "credibility_share"]),
        xytext=offset,
        textcoords="offset points",
        ha="left" if offset[0] > 0 else "right",
        fontsize=7,
    )
credibility_scores = score_shares(true_share, share.to_numpy())
ax.annotate(
    f"credibility: $\\rho$ = {credibility_scores['spearman']:.3f}, "
    f"CRPS = {credibility_scores['crps']:.4f}",
    (0.03, 0.95),
    xycoords="axes fraction",
)
ax.legend(loc="lower right")
ax.set(
    xlim=limits,
    ylim=limits,
    xlabel="true share of the channel",
    ylabel="predicted share of the channel",
    title="Predicted against true campaign shares",
);

Fig. Predicted against true flight shares, with the 94% HDI on the credibility estimate. Points above the dashed line are over-credited and points below it are under-credited. lookalike_conversions is the largest miss, under-credited by 0.07 in share, followed by spring_promo and video_product_launch over-credited by about 0.05 each. broad_prospecting lands almost exactly and brand_always_on within 0.02. The naive clicks split (squares) is further from the line on six of the seven campaigns.

The same comparison, campaign by campaign, with the interval laid flat so the coverage count can be read off directly:

# rows sorted by true share, largest first
forest_order = list(results.index)
pc = azp.plot_forest(
    xr.Dataset({"share": share.sel(campaign=forest_order)}),
    sample_dims=["sample"],
    combined=True,
    point_estimate="mean",
    ci_kind="hdi",
    ci_probs=(0.5, CI_PROB),
    labels=["campaign"],
    figure_kwargs={"figsize": (11, 6.5), "layout": "constrained"},
)
ax = pc.viz["/"]["plot"].sel(column="forest").item()
# y position of each campaign's row, read back from the point-estimate markers
row_y = {
    campaign: pc.viz["point_estimate"]["share"]
    .sel(campaign=campaign)
    .item()
    .get_offsets()[0, 1]
    for campaign in forest_order
}
for campaign in forest_order:
    y = row_y[campaign]
    ax.scatter(
        true_share[campaign],
        y,
        marker="|",
        s=500,
        linewidths=3,
        color="black",
        zorder=5,
    )
    ax.scatter(click_split[campaign], y, marker="s", color="C1", s=55, zorder=5)
    ax.scatter(
        baselines["spend share"][campaign], y, marker="^", color="C2", s=60, zorder=5
    )
handles = [
    Line2D([0], [0], color="C0", lw=4, label="credibility posterior: 50% HDI (thick)"),
    Line2D(
        [0],
        [0],
        color="C0",
        lw=1.5,
        marker="o",
        label=f"{CI_PROB:.0%} HDI (thin) and posterior mean",
    ),
    Line2D(
        [0],
        [0],
        color="black",
        marker="|",
        ms=14,
        mew=3,
        ls="none",
        label="hidden truth",
    ),
    Line2D([0], [0], color="C1", marker="s", ls="none", label="naive clicks split"),
    Line2D([0], [0], color="C2", marker="^", ls="none", label="spend share"),
]
ax.legend(handles=handles, loc="lower right", fontsize=10)
ax.set(xlabel="share of the channel's credit")
fig = pc.viz["/"]["figure"].values.item()
fig.suptitle(
    "Who gets the credit: credibility posterior against the truth and two heuristics",
    fontsize=14,
);

Fig. A forest plot of the credibility posterior for each campaign’s share of the channel, sorted by the true share. The thick blue segment is the 50% HDI, the thin segment the 94% HDI, and the open circle the posterior mean. The black tick is the hidden truth, the orange square the naive clicks split and the green triangle the spend share. The truth tick sits inside the thin segment for five of the seven campaigns; lookalike_conversions and spring_promo are the two misses, and both marginal cases are close: spring_promo misses its interval by 0.003 in share while video_product_launch makes it by 0.001. Where the square sits far from the tick and the circle sits between them, credibility has moved the answer in the right direction without going far enough.

Sweeping the dial by hand

The most useful validation in the notebook grades the method rather than one point estimate. The fitted model gave each campaign its own trust dial. Here we ask a simpler question: if every campaign were given the same trust level \(Z\), how would each setting score against the truth?

Two things to hold in mind while reading the figure.

• The two ends of the sweep are splits you already know. At \(Z = 0\) every campaign is judged by the collective rate alone, so the credit follows delivery: this is the daily-reach share from the scoreboard. At \(Z = 1\) every campaign is believed completely, so the credit follows the click split spread over each campaign’s daily delivery, renormalised day by day (Step 3 explains why this is not quite the plain click split from the scoreboard).

• The fitted model is not a point on this curve, because its dials differ by campaign, from 0.34 for niche_interest to 0.61 for brand_always_on. The shaded band marks that range, and the dashed horizontal lines mark the fitted model’s own scores, so you can see where campaign-specific dials land relative to one common dial. The fitted model never saw the truth, so this is a fair test.

def share_at_common_z(z: float) -> pd.Series:
    """Flight shares if every campaign were trusted by exactly the same amount z."""
    # The blend's (1 - z) * mu term is deliberately dropped: it is the same constant on
    # every theta, so it cancels in the normalisation that turns scores into shares.
    theta_z = xr.DataArray(
        z * y_direct.to_numpy()[None, :],
        dims=("sample", "campaign"),
        coords={"sample": [0], "campaign": campaigns},
    )
    mean_social = social_draws.mean("sample")
    shares = flight_share(allocate(mean_social, allocation_weights(theta_z, reach)))
    return shares.isel(sample=0).to_pandas().reindex(campaigns)


dial_grid = np.linspace(0.0, 1.0, 21)
sweep = pd.DataFrame(
    [
        {
            "Z": z,
            **{
                key: value
                for key, value in score_shares(
                    true_share, share_at_common_z(z).to_numpy()[None, :]
                ).items()
                if key in {"mae", "spearman"}
            },
        }
        for z in dial_grid
    ]
).set_index("Z")

# where the ranking changes, and which two campaigns swap places there
step_index = int(
    np.flatnonzero(np.abs(np.diff(sweep["spearman"].to_numpy())) > 1e-9)[0]
)
z_below, z_above = sweep.index[step_index], sweep.index[step_index + 1]
rank_below = share_at_common_z(z_below).rank(ascending=False)
rank_above = share_at_common_z(z_above).rank(ascending=False)
swapped = rank_below.index[rank_below != rank_above].tolist()
best_index = int(sweep["mae"].to_numpy().argmin())
print(
    f"Spearman steps between Z = {z_below:.2f} and Z = {z_above:.2f}; "
    f"the campaigns that change rank are {swapped}.\n"
    f"Fitted model (campaign-specific dials): MAE {credibility_scores['mae']:.4f}, "
    f"Spearman {credibility_scores['spearman']:.3f}."
)
sweep.iloc[[0, best_index, step_index, step_index + 1, len(sweep) - 1]].round(4)

Spearman steps between Z = 0.35 and Z = 0.40; the campaigns that change rank are ['spring_promo', 'video_product_launch'].
Fitted model (campaign-specific dials): MAE 0.0347, Spearman 0.893.

	spearman	mae
Z
0.00	0.8571	0.0322
0.20	0.8571	0.0310
0.35	0.8571	0.0318
0.40	0.8929	0.0321
1.00	0.8929	0.0373

fig, (ax_mae, ax_rho) = plt.subplots(
    nrows=2, sharex=True, figsize=(12, 8.5), layout="constrained"
)
for ax in (ax_mae, ax_rho):
    ax.axvspan(
        trust.min(),
        trust.max(),
        color="C0",
        alpha=0.15,
        label="fitted range of the campaign dials $Z_c$",
    )
ax_mae.plot(
    sweep.index,
    sweep["mae"],
    color="C0",
    marker="o",
    ms=4,
    label="one common $Z$ for every campaign",
)
ax_mae.axhline(
    credibility_scores["mae"],
    color="C3",
    ls="--",
    label=f"fitted model, campaign-specific dials (MAE {credibility_scores['mae']:.4f})",
)
ax_mae.axvline(
    sweep.index[best_index],
    color="gray",
    ls=":",
    label=f"best common dial, $Z$ = {sweep.index[best_index]:.2f} (MAE {sweep['mae'].min():.4f})",
)
ax_rho.plot(
    sweep.index,
    sweep["spearman"],
    color="C1",
    marker="s",
    ms=4,
    label="one common $Z$ for every campaign",
)
ax_rho.axhline(
    credibility_scores["spearman"],
    color="C3",
    ls="--",
    label=(
        "fitted model, campaign-specific dials "
        f"(Spearman {credibility_scores['spearman']:.3f})"
    ),
)
ax_mae.set(
    ylabel="MAE on shares (lower = better)",
    title="One common trust dial for every campaign: no single $Z$ wins on both scores",
)
ax_rho.set(
    ylabel="Spearman (higher = better)",
    xlabel="common trust dial $Z$ applied to every campaign",
    xticks=np.linspace(0, 1, 11),
)
ax_rho.text(
    0.0,
    -0.16,
    "$Z=0$: trust the average\n(delivery split)",
    transform=ax_rho.get_xaxis_transform(),
    ha="left",
    va="top",
    fontsize=10,
    style="italic",
)
ax_rho.text(
    1.0,
    -0.16,
    "$Z=1$: trust the dashboard\n(click split, renormalised daily)",
    transform=ax_rho.get_xaxis_transform(),
    ha="right",
    va="top",
    fontsize=10,
    style="italic",
)
ax_mae.margins(y=0.25)
ax_rho.margins(y=0.25)
ax_mae.legend(loc="upper left", fontsize=9)
ax_rho.legend(loc="lower right", fontsize=9);

Fig. One common trust dial swept from 0 to 1 and graded against the hidden truth. Top: the share error (MAE). Bottom: the ranking quality (Spearman). Three takeaways.

1. Trusting nobody and trusting everybody both lose. The MAE curve is U-shaped: 0.0322 at \(Z=0\) (the delivery split), 0.0373 at \(Z=1\) (the daily click split), and a minimum of 0.0310 near \(Z=0.20\). Some pooling is better than either extreme, which is the case for partial pooling in one picture.

2. The ranking is a step, not a curve. Below \(Z=0.40\), video_product_launch is ranked above spring_promo; from \(Z=0.40\) upward they swap, which matches the truth, and Spearman jumps from 0.857 to 0.893. No common dial gets the full order right, because niche_interest is under-credited at every setting.

3. The two scores disagree, and the fitted model sides with the ranking. The fitted dials (shaded, 0.34 to 0.61) sit past the step, so the model recovers the better order. It pays for that in MAE: its own score of 0.0347 (dashed line) is above the best common dial (0.0310) and slightly above the common-dial curve inside its own range. A common dial near 0.20 would have scored better on MAE while getting the order wrong.

The main takeaway: an interior amount of trust is the right kind of answer, the amount the data chose is defensible without ever looking at the truth, and no single number summarises “better” here, which is why the scoreboard reports ranking, error and calibration side by side.

Is the daily shape right?

The stacked figure in Step 3 takes for granted that the day-to-day shape of each campaign’s credit is right. Now that we may look at the truth, we can check: how well does each campaign’s daily credit track its true daily contribution?

daily_correlation = pd.Series(
    {
        campaign: float(
            np.corrcoef(
                daily_mean.sel(campaign=campaign).to_numpy(),
                data.truth[campaign].to_numpy(),
            )[0, 1]
        )
        for campaign in campaigns
    }
).sort_values(ascending=False)

focus = "flash_sale_burst"
focus_daily = contribution.sel(campaign=focus)
focus_bounds = az.hdi(focus_daily, dim="sample")

fig, ax = plt.subplots(figsize=(12, 6), layout="constrained")
ax.plot(
    reach.date,
    focus_daily.mean("sample").to_numpy(),
    color="C0",
    label="posterior mean",
)
ax.fill_between(
    reach.date,
    focus_bounds.sel(ci_bound="lower").to_numpy(),
    focus_bounds.sel(ci_bound="upper").to_numpy(),
    color="C0",
    alpha=0.3,
    label=f"{CI_PROB:.0%} HDI",
)
ax.plot(
    data.truth.index, data.truth[focus].to_numpy(), color="C1", ls="--", label="true"
)
ax.xaxis.set_major_locator(mdates.MonthLocator())
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b"))
ax.margins(y=0.20)
ax.legend(loc="upper left")
ax.set(
    title=f"Daily contribution posterior for {focus} "
    f"(daily correlation {daily_correlation[focus]:.2f})",
    xlabel="date",
    ylabel="contribution (outcome units)",
);

daily_correlation.rename("daily correlation with the truth").round(2)

flash_sale_burst         0.94
spring_promo             0.91
video_product_launch     0.89
broad_prospecting        0.86
brand_always_on          0.77
lookalike_conversions    0.56
niche_interest           0.28
Name: daily correlation with the truth, dtype: float64

Fig. The daily posterior for flash_sale_burst against the hidden truth. Notice the tail: the truth decays for about a week after the burst stops, while our allocation drops to exactly zero the day delivery ends, because the daily shape is inherited from daily reach.

The correlations measure the same thing for every campaign. The inherited shape holds well for five campaigns, from flash_sale_burst at 0.94 down to brand_always_on at 0.77, and poorly for lookalike_conversions at 0.56 and niche_interest at 0.28. That is the expected failure mode: the daily shape is inherited entirely from daily reach, so a campaign whose true impact does not track its delivery cannot be fixed by any reweighting of the level. The missing tail, however, is repairable, and the MMM already tells us how.

Extension: letting the credit outlive the flight

The daily allocation spends each day’s MMM credit in proportion to that day’s raw reach, so a campaign’s credit stops dead the moment its delivery does. The MMM itself disagrees: it was fit with a geometric adstock on the channel’s media variable, so the channel curve it hands us already contains carryover, revenue arriving days after the exposure that caused it. Splitting an adstocked curve by un-adstocked reach hands the tail of every burst to whichever campaigns happen to still be delivering.

This is also the answer to a fair objection raised about the method: the split should be aligned with the MMM’s own structure, not just anchored to its total. The repair is one line: push each campaign’s daily reach through the same transform the MMM used, geometric adstock, before computing the daily weights. Because the weights are renormalised every day, the split still reconciles to the MMM total exactly; only the shape of each campaign’s credit changes. And crucially, the decay rate is not a new free dial to agonise over. Your fitted MMM has a posterior for the channel’s adstock decay (adstock_alpha in PyMC-Marketing), and inheriting it is precisely what “aligned with the main model’s structure” means.

The data package here ships only the contribution curve, not the fitted decay, so this notebook does the next best thing: sweep the decay rate and grade each setting against the hidden truth, the same oracle-grading already used for the \(Z\) sweep. Two caveats to read alongside the table. This grading is a validation your own data cannot repeat, which is why the production recommendation is to inherit the MMM’s decay rather than tune it. And in this simulation each campaign’s true carryover decays at its own rate, so a single channel-level decay, which is all a channel-level MMM can offer, is itself an approximation.

def adstocked_reach(reach: xr.DataArray, alpha: float, l_max: int = 12) -> xr.DataArray:
    """Daily reach pushed through the MMM's geometric adstock along the date axis."""
    values = (
        geometric_adstock(
            ptx.as_xtensor(pt.as_tensor_variable(reach.to_numpy()), dims=reach.dims),
            alpha=alpha,
            l_max=l_max,
            normalize=True,
            dim="date",
        )
        .transpose(*reach.dims)
        .eval()
    )
    return xr.DataArray(values, dims=reach.dims, coords=reach.coords)


def daily_truth_correlation(contribution: xr.DataArray) -> pd.Series:
    """Correlation of each campaign's posterior-mean daily credit with the truth."""
    daily = contribution.mean("sample")
    return pd.Series(
        {
            campaign: float(
                np.corrcoef(
                    daily.sel(campaign=campaign).to_numpy(),
                    data.truth[campaign].to_numpy(),
                )[0, 1]
            )
            for campaign in campaigns
        }
    )


def evaluate_decay(alpha: float) -> dict[str, float]:
    """Allocate with adstocked reach at one decay rate and grade the result."""
    contribution_alpha = allocate(
        social_draws, allocation_weights(theta_draws, adstocked_reach(reach, alpha))
    )
    correlations = daily_truth_correlation(contribution_alpha)
    shares = flight_share(contribution_alpha).transpose("sample", "campaign")
    scores = score_shares(true_share, shares.to_numpy())
    return {
        "mean_daily_corr": float(correlations.mean()),
        **{f"corr_{campaign}": correlations[campaign] for campaign in campaigns},
        "spearman": scores["spearman"],
        "mae": scores["mae"],
        "crps": scores["crps"],
    }


decay_grid = (0.0, 0.2, 0.4, 0.6, 0.8)
corr_columns = [f"corr_{campaign}" for campaign in campaigns]
adstock_sweep = pd.DataFrame(
    [{"alpha": alpha, **evaluate_decay(alpha)} for alpha in decay_grid]
).set_index("alpha")
adstock_sweep.round(3)

	mean_daily_corr	corr_brand_always_on	corr_spring_promo	corr_lookalike_conversions	corr_video_product_launch	corr_flash_sale_burst	corr_niche_interest	corr_broad_prospecting	spearman	mae	crps
alpha
0.0	0.744	0.774	0.912	0.559	0.885	0.939	0.285	0.856	0.893	0.035	0.025
0.2	0.761	0.796	0.914	0.590	0.895	0.944	0.317	0.870	0.893	0.035	0.026
0.4	0.779	0.816	0.912	0.634	0.903	0.942	0.359	0.887	0.893	0.035	0.026
0.6	0.797	0.834	0.900	0.695	0.905	0.923	0.414	0.908	0.893	0.036	0.026
0.8	0.787	0.824	0.858	0.729	0.876	0.856	0.434	0.929	0.893	0.036	0.027

def mean_share_at(alpha: float) -> pd.Series:
    """Posterior-mean flight share of each campaign at one decay rate."""
    weights = allocation_weights(theta_draws, adstocked_reach(reach, alpha))
    return (
        flight_share(allocate(social_draws, weights))
        .mean("sample")
        .to_pandas()
        .reindex(campaigns)
    )


best_alpha = float(adstock_sweep["mean_daily_corr"].idxmax())
display_alpha = best_alpha if best_alpha > 0 else 0.4
corr_change = pd.Series(
    (
        adstock_sweep.loc[display_alpha, corr_columns]
        - adstock_sweep.loc[0.0, corr_columns]
    ).to_numpy(),
    index=campaigns,
).sort_values(ascending=False)
share_shift = float((mean_share_at(display_alpha) - mean_share_at(0.0)).abs().max())

print(
    f"Best decay rate by mean daily correlation: {best_alpha:.1f}. Largest change in any "
    f"campaign's mean flight share between no carryover and that decay: {share_shift:.4f}."
)
corr_change.rename(
    f"change in daily correlation, decay {display_alpha:.1f} vs 0"
).round(2)

Best decay rate by mean daily correlation: 0.6. Largest change in any campaign's mean flight share between no carryover and that decay: 0.0026.

lookalike_conversions    0.14
niche_interest           0.13
brand_always_on          0.06
broad_prospecting        0.05
video_product_launch     0.02
spring_promo            -0.01
flash_sale_burst        -0.02
Name: change in daily correlation, decay 0.6 vs 0, dtype: float64

The sweep favours carrying credit forward: the mean daily correlation with the truth rises from 0.744 with no carryover to 0.797 at a decay rate of 0.6 (row alpha = 0.0 reproduces the raw-reach allocation exactly, since a zero decay is the identity). The campaigns that gain are the ones whose credit had been stopping dead: lookalike_conversions +0.14, niche_interest +0.13, brand_always_on +0.06, broad_prospecting +0.05, video_product_launch +0.02. The correction is not free everywhere: flash_sale_burst loses 0.02 and spring_promo 0.01, because at a decay of 0.6 their tails are stretched slightly further than the truth’s. Window-total shares barely move: the largest change in any campaign’s mean flight share between no carryover and a decay of 0.6 is 0.0026, so this is a repair to when the credit lands, not to who gets it. The scoreboard columns confirm it: Spearman, MAE and CRPS are essentially flat across the decay grid.

adstock_contribution = allocate(
    social_draws, allocation_weights(theta_draws, adstocked_reach(reach, display_alpha))
)
focus = "flash_sale_burst"
focus_adstock = adstock_contribution.sel(campaign=focus)
focus_adstock_bounds = az.hdi(focus_adstock, dim="sample")

fig, ax = plt.subplots(figsize=(12, 7), layout="constrained")
ax.fill_between(
    reach.date,
    focus_adstock_bounds.sel(ci_bound="lower").to_numpy(),
    focus_adstock_bounds.sel(ci_bound="upper").to_numpy(),
    color="C0",
    alpha=0.3,
    label=f"adstocked, {CI_PROB:.0%} HDI",
)
ax.plot(
    reach.date,
    focus_adstock.mean("sample").to_numpy(),
    color="C0",
    label=f"adstocked allocation (decay {display_alpha:.1f})",
)
ax.plot(
    reach.date,
    contribution.sel(campaign=focus).mean("sample").to_numpy(),
    color="C2",
    ls=":",
    label="raw-reach allocation",
)
ax.plot(
    data.truth.index, data.truth[focus].to_numpy(), color="C1", ls="--", label="true"
)
ax.xaxis.set_major_locator(mdates.MonthLocator())
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b"))
ax.margins(y=0.20)
ax.legend(loc="upper left")
ax.set(
    title=f"Adstocked against raw-reach daily credit for {focus}",
    xlabel="date",
    ylabel="contribution (outcome units)",
);

Fig. The flash_sale_burst daily credit under the adstocked allocation, against the raw-reach allocation and the hidden truth. By construction the adstocked credit decays after the flight ends instead of stopping dead; how well that tail matches the truth is what the sweep table above measures. One boundary artefact is worth knowing about: adstocked reach ramps up from an empty history at the left edge of the window, while the true carryover of the always-on campaigns is already at steady state there, so the first stretch of the window is where this correction is least trustworthy.

Robustness checks

A sceptical reader has three fair questions at this point, and each deserves a test rather than a shrug.

1. The fitted dial’s band was wide. How much of the answer is the prior?

2. We asserted that unique reach is the evidence and reach-days are the exposure. Is that choice doing real work, or would the delivery measure do just as well?

3. The seven shares add to one, so the model’s assumption that the seven campaign signals are independent is not literally true. Does ignoring that hurt?

All three checks re-fit the model and re-grade the split with the Step 3 (raw reach) allocation, because the adstock extension showed that adstock changes when the credit lands and not who gets it. Questions 2 and 3 are, in different ways, both stress tests of the notebook’s thesis that reach is evidence. A fourth question, what to do when unique reach is not available at all, gets its own section after these three.

How much of the answer is the prior?

The dropdown in Step 2 explained why \(k\) and \(s^2\) are only weakly separated by seven campaigns, which is exactly the situation in which priors can quietly steer the answer. So we vary both priors that touch the dial, the scale of sigma_s and the width of the prior on \(\log k\), and grade every combination.

def evaluate_variant(
    volume: pd.Series,
    *,
    y: pd.Series | None = None,
    daily_exposure: xr.DataArray | None = None,
    sigma_s_scale: float = 0.5,
    log_k_sd: float = 0.75,
    run_seed: int = seed,
) -> dict[str, float]:
    """Fit one configuration end to end and grade the split it produces.

    By default the direct estimate is ``y_direct`` and the credit is spent over daily
    reach; the fallback section near the end swaps both.
    """
    y = y_direct if y is None else y
    daily_exposure = reach if daily_exposure is None else daily_exposure
    model = build_credibility_model(
        y, volume, sigma_s_scale=sigma_s_scale, log_k_sd=log_k_sd
    )
    variant = fit_credibility(model, sample_seed=run_seed, progressbar=False)
    variant_social, variant_theta = paired_draws(
        data.social, variant.posterior["theta"], pair_seed=run_seed
    )
    variant_share = flight_share(
        allocate(variant_social, allocation_weights(variant_theta, daily_exposure))
    ).transpose("sample", "campaign")
    dial = variant.posterior["Z"].mean(("chain", "draw")).to_numpy()
    return {
        "divergences": int(variant.sample_stats["diverging"].sum()),
        "k_median": float(variant.posterior["k"].median()) * float(volume.mean()),
        "Z_min": float(dial.min()),
        "Z_max": float(dial.max()),
        **score_shares(true_share, variant_share.to_numpy()),
    }


prior_sensitivity = pd.DataFrame(
    [
        {
            "sigma_s": sigma_s_scale,
            "sd_log_k": log_k_sd,
            **evaluate_variant(
                unique_reach, sigma_s_scale=sigma_s_scale, log_k_sd=log_k_sd
            ),
        }
        for sigma_s_scale in (0.4, 0.5, 0.7)
        for log_k_sd in (0.5, 0.75)
    ]
).set_index(["sigma_s", "sd_log_k"])
prior_sensitivity.round(4)

NUTS[nutpie]: [mu, log_k, sigma_s, theta_raw]

NUTS[nutpie]: [mu, log_k, sigma_s, theta_raw]

NUTS[nutpie]: [mu, log_k, sigma_s, theta_raw]

NUTS[nutpie]: [mu, log_k, sigma_s, theta_raw]

NUTS[nutpie]: [mu, log_k, sigma_s, theta_raw]

NUTS[nutpie]: [mu, log_k, sigma_s, theta_raw]

		divergences	k_median	Z_min	Z_max	spearman	mae	crps	coverage
sigma_s	sd_log_k
0.4	0.50	0	36681.0065	0.3367	0.6223	0.8929	0.0341	0.0251	5.0
	0.75	0	37336.0913	0.3405	0.6105	0.8929	0.0338	0.0245	6.0
0.5	0.50	0	37399.0742	0.3328	0.6177	0.8929	0.0348	0.0249	5.0
	0.75	0	37060.9792	0.3434	0.6138	0.8929	0.0347	0.0254	5.0
0.7	0.50	0	37831.6077	0.3297	0.6145	0.8929	0.0343	0.0247	6.0
	0.75	0	39093.8979	0.3335	0.6019	0.8929	0.0346	0.0252	5.0

The ranking is identical across all six settings, MAE stays between 0.0338 and 0.0348 and CRPS between 0.0245 and 0.0254, so the conclusions are not a prior artefact. The evidence threshold drifts from 36,681 to 39,094 people, about 7%, upward with the prior scale on sigma_s, which is the weak identification showing up as mild prior sensitivity exactly where the algebra said it would. No configuration produced a single divergence. Coverage sits at 5 of 7 in four settings and 6 of 7 in two, with no pattern, which is the fragility flagged in the grading section.

What we learn

The prior nudges the threshold \(k\) by a few thousand people and leaves the ranking, the share errors and the probabilistic score essentially untouched. The dial is set by the data, not by us.

Is unique reach really the right evidence?

The “Two quantities” callout in the EDA asserted that unique people, not reach-days, are the evidence. That is directly testable: swap the volume measure that sets \(Z\) from unique reach \(n_c\) to reach-days \(R_c\), at three seeds so that sampling noise is not mistaken for a result. The k_median column is in the units of whichever measure sets the dial, so it reads as people in the first block and as reach-days in the second.

evidence_choice = pd.DataFrame(
    [
        {
            "volume": name,
            "seed": run_seed,
            **evaluate_variant(volume, run_seed=run_seed),
        }
        for name, volume in [
            ("unique reach (evidence)", unique_reach),
            ("reach-days (delivery)", direct["reach_days"]),
        ]
        for run_seed in (42, 1, 7)
    ]
).set_index(["volume", "seed"])
evidence_choice.round(4)

NUTS[nutpie]: [mu, log_k, sigma_s, theta_raw]

NUTS[nutpie]: [mu, log_k, sigma_s, theta_raw]

NUTS[nutpie]: [mu, log_k, sigma_s, theta_raw]

NUTS[nutpie]: [mu, log_k, sigma_s, theta_raw]

NUTS[nutpie]: [mu, log_k, sigma_s, theta_raw]

NUTS[nutpie]: [mu, log_k, sigma_s, theta_raw]

		divergences	k_median	Z_min	Z_max	spearman	mae	crps	coverage
volume	seed
unique reach (evidence)	42	0	37060.9792	0.3434	0.6138	0.8929	0.0347	0.0254	5.0
	1	0	38103.5874	0.3360	0.6070	0.8929	0.0345	0.0254	4.0
	7	0	40297.8043	0.3279	0.5975	0.8929	0.0350	0.0253	4.0
reach-days (delivery)	42	1	587708.0941	0.0667	0.7045	0.8929	0.0327	0.0251	3.0
	1	0	573633.5694	0.0677	0.7096	0.8929	0.0319	0.0246	3.0
	7	0	582217.3447	0.0656	0.7026	0.8929	0.0324	0.0248	4.0

Read that honestly: this dataset does not settle the question, and on the graded scores it mildly favours the delivery measure. Driving the dial by reach-days gives the same ranking at 0.893, a CRPS of 0.0246 to 0.0251 against 0.0253 to 0.0254, slightly better MAE at 0.0319 to 0.0327 against 0.0345 to 0.0350, and coverage of 3 to 4 of 7 against 4 to 5 of 7. With only seven campaigns we cannot separate two measures that are this strongly correlated with each other. The reach-days fits are also the harder ones to sample, with 1 divergence across the three seeds against 0 for unique reach. What the table does establish is stability: three independently seeded fits move MAE and CRPS by at most 0.0008, and the coverage count by at most one campaign.

Look at the Z_min and Z_max columns, though. With reach-days as the evidence the dial spreads from 0.07 to 0.71, because reach-days span a factor of 47 across campaigns where unique reach spans 3.5. In practice that means brand_always_on would be believed almost fully, not because it reached the most people but because it re-hit the same people 25 times each, while broad_prospecting would be nearly ignored despite reaching 24k distinct people in 30 days. That is the delivery measure rewarding how hard we pushed, not how much we learned.

What we learn

We are overriding a small score difference on conceptual grounds, and it is worth being explicit about that. A person re-reached tomorrow adds delivery but no new information about the campaign, so unique reach is the quantity that answers “how much do I know about this campaign?”, which is what a credibility volume measure is supposed to be. An MAE gap of 0.0024 on seven campaigns is not the kind of evidence that should redefine what evidence means.

Does ignoring the composition hurt?

The seven naive rates come from seven click shares that add to one, so they are slices of one pie: if one slice is too big, the others must be too small. The model treats the seven signals as independent and ignores that built-in push-and-pull. Here we fit two variants that respect the pie and grade them on the same scoreboard.

• A compositional twin with the same trust dial. Same model, same hierarchy, same evidence-based \(Z_c\), but the likelihood is written on the differences between each campaign’s signal and a reference campaign, which is what a set of shares actually pins down. If ignoring the composition mattered, this twin would move the split.

• A model of the click counts themselves. The textbook way to model a pie of counts is a Dirichlet-multinomial on the raw clicks, with the same hierarchy on effectiveness. The catch is what plays the role of evidence: in this likelihood a share’s noise shrinks with the number of clicks, so a campaign with thousands of clicks looks precisely measured whether or not its clicks mean anything. That is the strongest possible test of “reach is evidence”, because it swaps reach out for clicks.

For the statistically minded: the two likelihoods

ALR logistic-normal twin. Suppose each campaign’s unnormalised log signal is \(z_c=\theta_c+\varepsilon_c\) with the familiar reach-scaled noise \(\varepsilon_c\sim\text{Normal}(0,\psi_c)\), and the observed click shares are the softmax of the \(z\)’s. Then the six observable contrasts \(y_c-y_{\text{ref}}\) are multivariate normal with mean \(\theta_c-\theta_{\text{ref}}\) and covariance \(\operatorname{diag}(\psi_{-\text{ref}})+\psi_{\text{ref}}\mathbf{1}\mathbf{1}^{\top}\), and that rank-one term is exactly the induced negative dependence. The implementation never writes that matrix down: conditional on the reference campaign’s own noise \(\varepsilon_{\text{ref}}\), the contrasts are independent Normals sharing the \(-\varepsilon_{\text{ref}}\) shift, so one extra scalar latent reproduces the covariance exactly. Nothing else changes: same hierarchy, same \(\psi_c=s^2/n^{\text{rel}}_c\), same dial \(Z_c\). The likelihood now only identifies contrasts of \(\theta\), so \(\mu\) stays close to its prior, but the daily allocation is invariant to a common shift of \(\theta\) (the softmax eats it), so nothing downstream is lost. It also makes explicit that the seven \(y_c\) carry only six degrees of freedom, since \(\sum_c R_c\,e^{y_c}=T\) holds by construction.

Dirichlet-multinomial. \(\text{clicks}\sim\text{DirichletMultinomial}(N,\kappa\,p)\) with \(p_c\propto R_c\,e^{\theta_c}\) and the same hierarchy on \(\theta\). The concentration \(\kappa\) is asked to learn the overdispersion from a single observed composition, and whether that leaves any room for reach to set the trust is exactly what the comparison measures.

def build_compositional_model(
    y_direct: pd.Series,
    volume: pd.Series,
    *,
    sigma_s_scale: float = 0.5,
    log_k_sd: float = 0.75,
) -> pm.Model:
    """ALR logistic-normal twin of the credibility model, same trust dial."""
    campaigns = list(y_direct.index)
    volume_relative = (volume / volume.mean()).reindex(campaigns).to_numpy()
    y_alr = (y_direct.iloc[:-1] - y_direct.iloc[-1]).to_numpy()
    with pm.Model(coords={"campaign": campaigns, "contrast": campaigns[:-1]}) as model:
        mu = pm.Normal("mu", mu=float(y_direct.mean()), sigma=1.0)
        log_k = pm.Normal("log_k", 0.0, log_k_sd)
        k = pm.Deterministic("k", pm.math.exp(log_k))
        sigma_s = pm.HalfNormal("sigma_s", sigma_s_scale)
        s2 = pm.Deterministic("s2", sigma_s**2)
        pm.Deterministic("tau", pm.math.sqrt(s2 / k))
        psi = s2 / volume_relative
        pm.Deterministic("Z", volume_relative / (volume_relative + k), dims="campaign")
        theta_raw = pm.Normal("theta_raw", 0.0, 1.0, dims="campaign")
        theta = pm.Deterministic(
            "theta", mu + pm.math.sqrt(s2 / k) * theta_raw, dims="campaign"
        )
        # the reference campaign's own noise, shared by every contrast: this scalar
        # latent IS the rank-one negative dependence, without writing the matrix down
        eps_ref = pm.Normal("eps_ref", 0.0, pm.math.sqrt(psi[-1]))
        pm.Normal(
            "y_alr_obs",
            mu=theta[:-1] - theta[-1] - eps_ref,
            sigma=pm.math.sqrt(psi[:-1]),
            observed=y_alr,
            dims="contrast",
        )
    return model


def build_dirichlet_multinomial_model(
    clicks: pd.Series,
    exposure: pd.Series,
    *,
    sigma_s_scale: float = 0.5,
    log_k_sd: float = 0.75,
) -> pm.Model:
    """Dirichlet-multinomial on click counts, same hierarchy on theta."""
    campaigns = list(clicks.index)
    log_exposure = np.log((exposure / exposure.mean()).reindex(campaigns).to_numpy())
    with pm.Model(coords={"campaign": campaigns}) as model:
        mu = pm.Normal("mu", 0.0, 1.0)
        log_k = pm.Normal("log_k", 0.0, log_k_sd)
        k = pm.Deterministic("k", pm.math.exp(log_k))
        sigma_s = pm.HalfNormal("sigma_s", sigma_s_scale)
        s2 = pm.Deterministic("s2", sigma_s**2)
        pm.Deterministic("tau", pm.math.sqrt(s2 / k))
        theta_raw = pm.Normal("theta_raw", 0.0, 1.0, dims="campaign")
        theta = pm.Deterministic(
            "theta", mu + pm.math.sqrt(s2 / k) * theta_raw, dims="campaign"
        )
        kappa = pm.LogNormal("kappa", mu=np.log(1_000.0), sigma=2.0)
        p = pm.math.softmax(theta + log_exposure)
        pm.DirichletMultinomial(
            "clicks_obs",
            n=int(clicks.sum()),
            a=kappa * p,
            observed=clicks.to_numpy(),
            dims="campaign",
        )
    return model


def variant_share(idata: xr.DataTree, *, pair_seed: int) -> xr.DataArray:
    """Flight-share draws implied by one fitted variant's theta posterior."""
    variant_social, variant_theta = paired_draws(
        data.social, idata.posterior["theta"], pair_seed=pair_seed
    )
    return flight_share(
        allocate(variant_social, allocation_weights(variant_theta, reach))
    ).transpose("sample", "campaign")


compositional_idata = fit_credibility(
    build_compositional_model(y_direct, unique_reach), progressbar=False
)
compositional_share = variant_share(compositional_idata, pair_seed=seed)
dm_idata = fit_credibility(
    build_dirichlet_multinomial_model(direct["clicks"], direct["reach_days"]),
    progressbar=False,
)
dm_share = variant_share(dm_idata, pair_seed=seed)

variant_scoreboard = pd.DataFrame(
    [
        {
            "method": name,
            "divergences": int(fit.sample_stats["diverging"].sum()),
            **score_shares(true_share, draws.to_numpy()),
        }
        for name, fit, draws in [
            ("independent Normal (main)", idata, share),
            ("ALR logistic-normal", compositional_idata, compositional_share),
            ("Dirichlet-multinomial", dm_idata, dm_share),
        ]
    ]
).set_index("method")
variant_scoreboard.assign(
    coverage=lambda df: df["coverage"].map("{:.0f} of 7".format)
).round(4)

NUTS[nutpie]: [mu, log_k, sigma_s, theta_raw, eps_ref]

NUTS[nutpie]: [mu, log_k, sigma_s, theta_raw, kappa]

	divergences	spearman	mae	crps	coverage
method
independent Normal (main)	0	0.8929	0.0347	0.0254	5 of 7
ALR logistic-normal	0	0.8929	0.0343	0.0251	4 of 7
Dirichlet-multinomial	10	0.8929	0.0360	0.0299	3 of 7

alr_share_gap = float(
    (compositional_share.mean("sample").to_pandas() - share.mean("sample").to_pandas())
    .abs()
    .max()
)
dm_click_gap = float((dm_share.mean("sample").to_pandas() - click_split).abs().max())
main_click_gap = float((share.mean("sample").to_pandas() - click_split).abs().max())
kappa_posterior = dm_idata.posterior["kappa"].to_numpy().ravel()
kappa_prior_iqr = float(
    np.exp(np.log(1_000.0) + 2.0 * 0.6745) - np.exp(np.log(1_000.0) - 2.0 * 0.6745)
)
kappa_posterior_iqr = float(
    np.quantile(kappa_posterior, 0.75) - np.quantile(kappa_posterior, 0.25)
)

pd.Series(
    {
        "largest mean-share gap, ALR twin vs main": alr_share_gap,
        "largest mean-share gap, main vs naive click split": main_click_gap,
        "largest mean-share gap, Dirichlet-multinomial vs naive click split": dm_click_gap,
        "kappa prior interquartile range": kappa_prior_iqr,
        "kappa posterior interquartile range": kappa_posterior_iqr,
    },
    name="value",
).round(4)

largest mean-share gap, ALR twin vs main                                 0.0009
largest mean-share gap, main vs naive click split                        0.0422
largest mean-share gap, Dirichlet-multinomial vs naive click split       0.0304
kappa prior interquartile range                                       3594.0704
kappa posterior interquartile range                                   3642.4506
Name: value, dtype: float64

Reading the two tables in both directions.

The compositional twin changes remarkably little. The two models are practically the same split: the largest difference in any campaign’s mean share is 0.0009. On the graded scores the twin lands at MAE 0.0343, CRPS 0.0251 and Spearman 0.893 against 0.0347, 0.0254 and 0.893 for the independent likelihood, with no divergences in either. The push-and-pull between the slices is real, but with seven campaigns and noise this large it barely moves the split. That is a useful negative result: the independent likelihood is not what stands between this method and better accuracy.

The Dirichlet-multinomial is the cautionary tale. Its mean shares end up within 0.0304 of the naive click split, against 0.0422 for the credibility split, so it drifts back toward trusting the dashboard. It scores MAE 0.0360, CRPS 0.0299 and Spearman 0.893, with 10 divergences, so it is both less accurate and harder to sample. And its concentration parameter, the one number that could have added scepticism, is left almost where the prior put it: a posterior interquartile range of 3,642 against 3,594 in the prior. One observed composition cannot teach it how much to distrust clicks, so the amount of trust ends up set by a prior nobody can argue with, not by reach.

What we learn

Respecting the pie is not what limits accuracy here. Swapping reach for clicks as the evidence is what breaks the method: thousands of clicks look precise, but “precise” is not “right”. The noise that matters is the gap between clicks and incrementality, and it does not shrink with click volume. Reach measures how much we know; clicks measure how loud the signal is. That is why reach does the disciplining in this notebook, and why it can, precisely because it is not the quantity being divided up.

What if you cannot get unique reach?

Everything above rests on one de-duplicated number per campaign: the unique people it reached over its flight. Before treating that as a hard requirement, be precise about what the model needs. It does not need the expanding-window series. That series gives the audience growth curve and the saturation diagnostic, which are useful but optional. The model needs only the flight total \(n_c\), and any platform that reports reach for a date range de-duplicates within that range, so a single report per campaign over its full flight is enough. The fallback question is therefore narrower than it sounds: what if all you have is the daily export, with its daily reach, impressions and clicks, and no de-duplicated number at all?

The daily export still brackets the answer. Unique reach can never be smaller than the largest single day’s reach (those people were all distinct on that day) and never larger than the reach-days total (that is the count if nobody was ever reached twice):

\[\max_t \text{reach}_{ct} \;\le\; n_c \;\le\; R_c = \sum_t \text{reach}_{ct}.\]

Here is where the seven campaigns sit inside their brackets. In practice the unique_reach column is the number you would not have.

reach_bounds = pd.DataFrame(
    {
        "peak_daily_reach": totals["peak_daily_reach"],
        "unique_reach": unique_reach.astype(int),
        "reach_days": totals["reach_days"],
        "unique_over_peak": unique_reach / totals["peak_daily_reach"],
        "reach_days_over_unique": totals["reach_days"] / unique_reach,
    }
).sort_values("unique_reach", ascending=False)
spread = (
    reach_bounds[["peak_daily_reach", "unique_reach", "reach_days"]]
    .agg(lambda column: column.max() / column.min())
    .rename("largest campaign / smallest campaign")
)
display(spread.round(1))
reach_bounds.round(2)

peak_daily_reach     9.3
unique_reach         3.5
reach_days          46.8
Name: largest campaign / smallest campaign, dtype: float64

	peak_daily_reach	unique_reach	reach_days	unique_over_peak	reach_days_over_unique
brand_always_on	17748	62386	1541137	3.52	24.70
lookalike_conversions	8024	42733	597676	5.33	13.99
spring_promo	10423	41719	438462	4.00	10.51
flash_sale_burst	11227	34802	139620	3.10	4.01
video_product_launch	9365	33427	519340	3.57	15.54
broad_prospecting	1906	24140	32921	12.67	1.36
niche_interest	3426	17987	254475	5.25	14.15

The brackets are wide. The true unique reach is 3 to 13 times the peak day and 1.4 to 25 times smaller than the reach-days total, and the ratio differs by campaign in a way that tracks frequency: brand_always_on re-reaches the same people for months, broad_prospecting barely revisits anyone. That is exactly why the two bounds make different evidence measures. Reach-days reward frequency, so the always-on campaign looks like it carries 47 times the evidence of the smallest one; peak daily reach ignores frequency altogether and compresses the spread to 9 times, closer to the 3.5 times of unique reach.

We can test each fallback the same way as everything else: swap it into the model as the evidence measure, re-fit, re-allocate and grade against the hidden truth. The evaluate_variant helper from the robustness checks already does this, and two of its arguments let us also swap the direct estimate and the daily exposure, which the impressions-only recipe needs.

• Reach-days as evidence. Reach only, no de-duplication: the delivery measure already tested in the robustness checks, repeated here for comparison.

• Peak daily reach as evidence. Reach only, using the lower end of the bracket as a conservative stand-in for people.

• Impressions only. No reach column at all. The direct estimate becomes the click split per impression, \(y_c = \log(T\cdot\text{clickshare}_c/\text{impressions}_c)\), which is a log click-through rate up to a constant; impressions set the trust dial, and the daily credit is spent over daily impressions instead of daily reach.

• Clicks only is the one case with no experiment to run. If clicks are both the signal and the exposure, the rate \(T \cdot \text{clickshare}_c / \text{clicks}_c\) is the same number for every campaign, so there is no own story left to shrink and the only split available is the click split itself (the naive clicks share row of the scoreboard). Making clicks the evidence instead is the Dirichlet-multinomial trap of the previous section: thousands of clicks look precise, whether or not they mean anything.

impressions = totals["impressions"].astype(float)
y_impressions = naive_log_rate(channel_total, direct["click_share"], impressions)
daily_impressions = daily_matrix(data.metrics, campaigns, "impressions")

fallback_recipes = {
    "unique reach (full recipe)": {"volume": unique_reach},
    "reach-days as evidence": {"volume": direct["reach_days"]},
    "peak daily reach as evidence": {
        "volume": totals["peak_daily_reach"].astype(float)
    },
    "impressions only (CTR signal)": {
        "volume": impressions,
        "y": y_impressions,
        "daily_exposure": daily_impressions,
    },
}
fallback_scoreboard = pd.DataFrame(
    [
        {"recipe": name, **evaluate_variant(**kwargs)}
        for name, kwargs in fallback_recipes.items()
    ]
).set_index("recipe")
fallback_scoreboard.assign(
    coverage=lambda df: df["coverage"].map("{:.0f} of 7".format)
).round(4)

NUTS[nutpie]: [mu, log_k, sigma_s, theta_raw]

NUTS[nutpie]: [mu, log_k, sigma_s, theta_raw]

NUTS[nutpie]: [mu, log_k, sigma_s, theta_raw]

NUTS[nutpie]: [mu, log_k, sigma_s, theta_raw]

	divergences	k_median	Z_min	Z_max	spearman	mae	crps	coverage
recipe
unique reach (full recipe)	0	37060.9792	0.3434	0.6138	0.8929	0.0347	0.0254	5 of 7
reach-days as evidence	1	587708.0941	0.0667	0.7045	0.8929	0.0327	0.0251	3 of 7
peak daily reach as evidence	0	8117.5777	0.2121	0.6681	0.8929	0.0350	0.0249	5 of 7
impressions only (CTR signal)	2	875652.8889	0.0463	0.7474	0.8571	0.0360	0.0287	2 of 7

Read the table against the full recipe in the first row. The k_median column is in the units of whichever measure sets the dial: people for the first and third rows, reach-days for the second, impressions for the fourth.

Peak daily reach is a serviceable stand-in. It reproduces the ranking (0.893), lands within 0.0005 of the full recipe on MAE (0.0350 against 0.0347) and CRPS (0.0249 against 0.0254), keeps 5 of 7 truths inside the interval, and samples without divergences. Its dials run from 0.21 to 0.67, a wider spread than unique reach gives (0.34 to 0.61), because a burst that hits 11k people in a day looks better evidenced than a slow campaign that never exceeds 3.5k. But it does not reward frequency, so it does not hand the always-on campaign the near-total trust that reach-days would.

Reach-days buys a little MAE with over-confidence. As in the robustness checks, it scores slightly better on MAE (0.0327) and CRPS (0.0251) but covers only 3 of 7 truths, spreads the dial from 0.07 to 0.70, and produced a divergence. The always-on campaign is believed almost fully because it re-hit the same people 25 times each, not because it reached more of them.

Impressions only is the weakest recipe. Ranking drops to 0.857, CRPS worsens to 0.0287, only 2 of 7 truths land inside the interval, and the sampler reports divergences. Two things go wrong at once: the signal is now a click-through rate, which is the currency most distorted by what a campaign was optimised for, and the credit is spent over impressions, so the daily shape follows how often a campaign was shown rather than how many people saw it.

What we learn

If you cannot get a de-duplicated reach report, prefer peak daily reach as the evidence measure: it is reach-only, conservative, and on this data it performs about as well as unique reach. Reach-days and impressions are usable, but they make trust follow delivery rather than people, which is the very thing the method is meant to avoid, and impressions-only also swaps the signal for a click-through rate. Clicks alone leave nothing to discipline. Whatever the fallback, the structural price is the same: the dial no longer measures how many people you learned from, only how hard you pushed. Get the de-duplicated number if you possibly can; it is one report per campaign.

What it can and cannot do

What it is genuinely good at

• It refuses to be fooled by small, loud campaigns. The less evidence a campaign has, the harder its wild self-report is reined in, so bursty campaigns stop hijacking the credit. flash_sale_burst is the case in point: its naive rate per reach-day is the highest in the account, and the model cuts it by about a third (0.40 log units, a factor of 1.5).

• It never breaks the budget. The campaign pieces sum back to the MMM’s channel total exactly, every day, to numerical precision. It divides the number you trust rather than inventing a new one. Note this is a property of the construction, any renormalised split has it, so treat it as a floor the method clears rather than evidence that the split is decision-grade.

• It is the best probabilistic forecast of everything tested here, at a CRPS of 0.0254 in share units against 0.0321 for the best simple heuristic.

• The dial it picks is defensible without the oracle. Sweeping a common \(Z\) shows partial pooling beats both extremes, and the fitted dials land in the interior and recover the campaign ranking without ever seeing the truth.

• It is built on the safest idea in the toolbox, partial pooling, so it behaves well and is easy to explain to whoever asks.

What it cannot do (read before trusting)

• It cures noise, not bias. It fixes the “small campaigns are jumpy” problem. It does not fix the fact that clicks are not incrementality. lookalike_conversions stays under-credited by 0.07 in share (7 percentage points of the channel’s credit) and niche_interest by 0.03, because no method that only sees delivery data can recover what delivery data does not contain.

• A campaign ranking is not reallocation advice. The split hands out average credit for spend already delivered, not the marginal return on the next dollar. No per-campaign response curve is fitted anywhere in the method, and saturation is estimated at channel grain, so reading “move budget to the top-ranked campaign” out of the split silently assumes that a campaign’s contribution scales with its own spend and that the channel’s saturation curve is indifferent to which campaign the dollar went to. Neither assumption is tested here, and this limit is separate from the clicks-are-not-incrementality one: it would still hold if the proxy measured incrementality perfectly. For the forward-looking question, take the split’s output to the budget allocation workflow, which models marginal response explicitly.

• On point accuracy it does not beat a delivery split here. MAE is 0.0347 against 0.0321 for the daily-reach share, and a common dial near 0.20 would have done better still on MAE while getting the ranking wrong. Its measurable gains are ranking, calibration and a defensible trust weight, not a jump in point accuracy.

• The daily shape is inherited, not estimated. Within a campaign the day-to-day pattern is its daily reach, which tracks the truth well for five campaigns and badly for lookalike_conversions at 0.56 and niche_interest at 0.28. The adstock extension shows how much of that is repairable by inheriting the MMM’s carryover, and a campaign whose true impact simply does not track its delivery stays beyond any reweighting.

• The uncertainty is a lower bound. Five of seven true shares land inside the 94% HDI, but the interval carries parameter uncertainty only, and the clicks-to-incrementality gap is nowhere in the model.

• The likelihood ignores that the target is compositional, and we tested the fix. The compositional twin barely moves the split, so this shortcut is not what limits accuracy here; the Dirichlet-multinomial alternative trades the reach-based dial for a prior-driven one and drifts toward the click split.

• It measures audience overlap, it does not correct for it. The overlap matrix in the data section is a diagnostic, nothing downstream consumes it. The credibility weights treat each campaign’s unique reach as independent evidence, so campaigns contesting the same audience are trusted as if their reach were disjoint. Accounting for correlated reach is a possible extension, not something this model already does.

• Evidence can be confounded with a truncated flight. broad_prospecting is the only campaign still acquiring new people when the window closes, so its unique reach measures 30 days of growth rather than a saturated audience. It is discounted for being young as much as for being small, and nothing in the model knows the difference.

• It needs a people count, and the fallbacks cost something. Without a de-duplicated reach number, peak daily reach is the safest stand-in and performs about as well here; reach-days and impressions make trust follow delivery instead of people, and clicks alone leave nothing to discipline.

• Only seven campaigns. The dial’s tuning is itself uncertain, the posterior for \(k\) says so, and no single synthetic dataset can settle a comparison against any baseline.

Shown daily, applies weekly

Everything above runs on a daily MMM with daily campaign exports. On a weekly MMM the same machinery applies one level coarser: the exposure base becomes reach-weeks, the allocation weights renormalise per week, and the expanding-window unique-reach query is unchanged. What you lose is resolution, not validity. With far fewer periods per flight there is no within-week shape to inherit, short bursts blur into their week, and campaign flights overlap more heavily inside each period, so the split leans harder on the trust dial and less on timing. The evidence measure, unique people reached, does not change with the grain, which is one more reason to prefer it over volume.

Run it on your own data

1. Pull de-duplicated reach per campaign over its flight; one report per campaign is enough for the model. For the audience growth curve, pull the expanding-window series instead: fix the report’s start date and extend the end date one day at a time, so the day-over-day differences are genuinely new people (the cum_reach table above). If no de-duplicated number is available at all, use each campaign’s peak daily reach as the evidence, as the fallback section shows.

2. Build any naive split of the MMM channel total from what your dashboard already gives you, whether clicks, conversions or last-touch revenue. It is allowed to be bad, that is the point. Express it per reach-day, using the same exposure you will allocate over. In the notebook this is direct_estimates, and the MMM total comes from mmm.idata.posterior["channel_contribution"] selected on the channel.

3. Fit the model with build_credibility_model and fit_credibility, and read off each campaign’s trust dial \(Z_c\). Even before the split, \(Z\) alone tells you which campaign numbers deserve to be taken seriously.

4. Allocate the MMM channel curve daily with daily_reach, allocation_weights and allocate, scoring each campaign by its trust-adjusted rate times its daily reach. Push that daily reach through adstocked_reach first, with the channel’s fitted decay (adstock_alpha in your MMM’s posterior), so a campaign’s credit can outlive its flight the way the MMM says the channel’s impact does. The pieces sum to the MMM total by construction either way.

5. Sweep the dial against any ground truth you have, even one lift test, with share_at_common_z. It is the cheapest check that the answer is not an artefact of one fit.

The bottom line. This is an uncertainty-aware way to split a channel’s MMM credit across campaigns, most valuable exactly where platform numbers are least reliable, on the small and bursty campaigns. Two structural refinements were tested above and are cheap to keep: spending the credit over adstocked reach aligns the split’s timing with the carryover the MMM already believes in, and the compositional likelihood check shows the independence shortcut is not what limits accuracy. Its one hard limit is that delivery data cannot see true incrementality. To break through it, pair the method with a lift test: feed a campaign’s measured incremental lift in as a prior, and the trust dial will do the rest.

Related notebooks. The Media Mix Modeling (MMM) in Python: End-to-End Example shows where the channel contribution posterior comes from; the adstock guide covers the geometric adstock inherited in the extension; the budget allocation notebook is the natural next step once campaign credit is settled; the lift test calibration notebook shows how experimental evidence enters an MMM; and the PIE notebook predicts campaign incrementality from experiments directly.

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

Last updated: Wed, 26 Aug 2026

Python implementation: CPython
Python version       : 3.14.2
IPython version      : 9.15.0

pymc_marketing: 1.1.0
pytensor      : 3.2.4

IPython       : 9.15.0
arviz         : 1.2.0
arviz_plots   : 1.2.0
matplotlib    : 3.10.9
numpy         : 2.4.6
pandas        : 2.3.3
pymc          : 6.2.0
pymc_marketing: 1.1.0
pytensor      : 3.2.4
scipy         : 1.18.0
seaborn       : 0.13.2
xarray        : 2026.4.0

Watermark: 2.6.0