Media Mix Models

Build MMMs of increasing complexity: adstock and saturation transforms, marketing funnel mediation, and hierarchical geo-varying effects.

A marketing analyst runs a regression of weekly sales on TV and digital spend. The model says TV is weak and digital is strong. Based on this, the team shifts budget from TV to digital — and six months later, sales are down. What went wrong?

Several things might have been missed:

• TV effects carry over across weeks and saturate at high spend levels — a linear regression captures neither
• TV spend also drove search traffic, which drove sales — the regression gave search credit for TV’s work
• TV was twice as effective in urban markets, but the model averaged over all regions

This notebook addresses all three problems with progressively richer models.

Model	What it adds	Key pathmc features
1. Transforms	Adstock carry-over + logistic saturation	adstock(), logistic_saturation(), panel mode
2. Marketing funnel	Direct vs indirect effects through a mediator	Labeled coefficients, := defined parameters, effect()
3. Hierarchical geo	Region-varying media coefficients	pooling={"slopes": [...]}, random slopes

Each model adds a new dimension to the analysis. You can stop after any section and apply what you’ve learned.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import arviz as az
import pymc as pm
import pathmc

FIG_WIDTH = 8
FIG_HEIGHT = 4

COLOR_TV = "#2171b5"
COLOR_DIGITAL = "#e6550d"
COLOR_BASELINE = "#636363"

COLOR_BRAND = "#2171b5"
COLOR_PERF = "#e6550d"
COLOR_SEARCH = "#31a354"
COLOR_INDIRECT = "#756bb1"

Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.

Model 1: Adstock and saturation transforms

Two empirical regularities shape how advertising affects sales:

1. Adstock: effects carry over across time periods — a TV ad this week still influences next week’s sales. Geometric carry-over: y_t = x_t + \text{decay} \cdot y_{t-1}.
2. Logistic saturation: each additional unit of spend produces less incremental effect. y = (1 - \exp(-\lambda \cdot x)) / (1 + \exp(-\lambda \cdot x)), which is \tanh(\lambda x / 2).

This model captures both transforms in a panel of 4 regions over 30 weeks.

The causal structure

Figure 1: MMM DAG with adstock and logistic saturation transforms applied to TV and digital spend before entering the sales regression.

Simulate panel data

rng = np.random.default_rng(42)

regions = ["North", "South", "East", "West"]
n_weeks = 30
true_intercepts = {"North": 50, "South": 60, "East": 45, "West": 55}
true_decay_tv = 0.7
true_decay_dig = 0.5
true_lam_tv = 0.3
true_lam_dig = 0.4
true_b_tv = 8.0
true_b_dig = 6.0

rows = []
for region in regions:
    adstocked_tv = 0.0
    adstocked_dig = 0.0
    for week in range(1, n_weeks + 1):
        tv = rng.uniform(10, 50)
        digital = rng.uniform(5, 30)
        adstocked_tv = tv + true_decay_tv * adstocked_tv
        adstocked_dig = digital + true_decay_dig * adstocked_dig
        sat_tv = (1 - np.exp(-true_lam_tv * adstocked_tv)) / (
            1 + np.exp(-true_lam_tv * adstocked_tv)
        )
        sat_dig = (1 - np.exp(-true_lam_dig * adstocked_dig)) / (
            1 + np.exp(-true_lam_dig * adstocked_dig)
        )
        sales = (
            true_intercepts[region]
            + true_b_tv * sat_tv
            + true_b_dig * sat_dig
            + 0.1 * week
            + rng.normal(scale=1.5)
        )
        rows.append({
            "region": region,
            "week": week,
            "tv": tv,
            "digital": digital,
            "trend": week,
            "sales": sales,
        })

df = pd.DataFrame(rows)
print(f"Panel: {len(regions)} regions × {n_weeks} weeks = {len(df)} rows")
df.head()

Panel: 4 regions × 30 weeks = 120 rows

	region	week	tv	digital	trend	sales
0	North	1	40.958242	15.971961	1	65.205473
1	North	2	37.894721	7.354434	2	62.220828
2	North	3	40.445588	24.651608	3	64.274769
3	North	4	28.015438	14.269951	4	65.566626
4	North	5	35.754605	25.569040	5	65.201263

Visualise sales by region

REGION_COLORS = {
    "North": "#2171b5",
    "South": "#e6550d",
    "East": "#31a354",
    "West": "#756bb1",
}

fig, ax = plt.subplots(figsize=(FIG_WIDTH, FIG_HEIGHT))
for region in regions:
    d = df[df["region"] == region]
    ax.plot(d["week"], d["sales"], label=region, color=REGION_COLORS[region], alpha=0.8)
ax.set_xlabel("Week")
ax.set_ylabel("Sales")
ax.legend()
plt.tight_layout()
plt.show()

Figure 2: Simulated weekly sales for four regions with adstock carry-over and saturation effects.

Specify and fit the model

The DSL supports transforms directly in the formula. Adstock and logistic saturation compose naturally: logistic_saturation(adstock(tv, decay=theta_tv), lam=lam_tv). pooling="partial" adds a random intercept per region and drops the fixed intercept. The formula starts with 0 + so equations() does not print a phantom 1 + intercept that is not in the design matrix.

spec = """
sales ~ 0 + b_tv*logistic_saturation(adstock(tv, decay=theta_tv), lam=lam_tv)
      + b_dig*logistic_saturation(adstock(digital, decay=theta_dig), lam=lam_dig)
      + trend
"""

model = pathmc.model(
    spec,
    data=df,
    panel={"unit": "region", "time": "week"},
    pooling="partial",
)
model.graph()

Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.

model.equations()

\begin{aligned} \beta_{sales} &\sim \text{Normal}(mu=0,\, sigma=10) \\ \sigma_{sales} &\sim \text{HalfNormal}(sigma=1) \\ \mu_{alpha,sales} &\sim \text{Normal}(mu=0,\, sigma=10) \\ \sigma_{alpha,sales} &\sim \text{HalfNormal}(sigma=1) \\ \alpha_{sales} &\sim \text{Normal}(mu\_alpha,\, sigma\_alpha) \\ \theta_{tv} &\sim \text{Beta}(alpha=2,\, beta=2) \\ \lambda_{tv} &\sim \text{HalfNormal}(sigma=1) \\ \theta_{dig} &\sim \text{Beta}(alpha=2,\, beta=2) \\ \lambda_{dig} &\sim \text{HalfNormal}(sigma=1) \\[6pt] \mu_{sales} &= b_{tv} \cdot \operatorname{logistic\_saturation}(\operatorname{adstock}(\mathrm{tv},\, \theta_{tv}),\, \lambda_{tv}) + b_{dig} \cdot \operatorname{logistic\_saturation}(\operatorname{adstock}(\mathrm{digital},\, \theta_{dig}),\, \lambda_{dig}) + \mathrm{trend} \\ \mathrm{sales} &\sim \text{Normal}(\mu_{sales},\, \sigma_{sales}) \end{aligned}

pm.model_to_graphviz(model.pymc_model)

/Users/benjamv/git/pathmc-worktrees/issue-538/.venv/lib/python3.12/site-packages/pytensor/link/c/cmodule.py:2053: UserWarning: `pytensor.config.cxx` is not an identifiable `g++` compiler. PyTensor will disable compiler optimizations specific to `g++`. At worst, this could cause slow downs.
Those parameters can be added manually via the `cxxflags` setting.
  warnings.warn(
/Users/benjamv/git/pathmc-worktrees/issue-538/.venv/lib/python3.12/site-packages/pytensor/link/c/cmodule.py:2053: UserWarning: `pytensor.config.cxx` is not an identifiable `g++` compiler. PyTensor will disable compiler optimizations specific to `g++`. At worst, this could cause slow downs.
Those parameters can be added manually via the `cxxflags` setting.
  warnings.warn(
/Users/benjamv/git/pathmc-worktrees/issue-538/.venv/lib/python3.12/site-packages/pytensor/link/c/cmodule.py:2053: UserWarning: `pytensor.config.cxx` is not an identifiable `g++` compiler. PyTensor will disable compiler optimizations specific to `g++`. At worst, this could cause slow downs.
Those parameters can be added manually via the `cxxflags` setting.
  warnings.warn(
/Users/benjamv/git/pathmc-worktrees/issue-538/.venv/lib/python3.12/site-packages/pytensor/link/c/cmodule.py:2053: UserWarning: `pytensor.config.cxx` is not an identifiable `g++` compiler. PyTensor will disable compiler optimizations specific to `g++`. At worst, this could cause slow downs.
Those parameters can be added manually via the `cxxflags` setting.
  warnings.warn(
/Users/benjamv/git/pathmc-worktrees/issue-538/.venv/lib/python3.12/site-packages/pytensor/link/c/cmodule.py:2053: UserWarning: `pytensor.config.cxx` is not an identifiable `g++` compiler. PyTensor will disable compiler optimizations specific to `g++`. At worst, this could cause slow downs.
Those parameters can be added manually via the `cxxflags` setting.
  warnings.warn(
/Users/benjamv/git/pathmc-worktrees/issue-538/.venv/lib/python3.12/site-packages/pytensor/link/c/cmodule.py:2053: UserWarning: `pytensor.config.cxx` is not an identifiable `g++` compiler. PyTensor will disable compiler optimizations specific to `g++`. At worst, this could cause slow downs.
Those parameters can be added manually via the `cxxflags` setting.
  warnings.warn(
/Users/benjamv/git/pathmc-worktrees/issue-538/.venv/lib/python3.12/site-packages/pytensor/link/c/cmodule.py:2053: UserWarning: `pytensor.config.cxx` is not an identifiable `g++` compiler. PyTensor will disable compiler optimizations specific to `g++`. At worst, this could cause slow downs.
Those parameters can be added manually via the `cxxflags` setting.
  warnings.warn(
/Users/benjamv/git/pathmc-worktrees/issue-538/.venv/lib/python3.12/site-packages/pytensor/link/c/cmodule.py:2053: UserWarning: `pytensor.config.cxx` is not an identifiable `g++` compiler. PyTensor will disable compiler optimizations specific to `g++`. At worst, this could cause slow downs.
Those parameters can be added manually via the `cxxflags` setting.
  warnings.warn(
/Users/benjamv/git/pathmc-worktrees/issue-538/.venv/lib/python3.12/site-packages/pytensor/link/c/cmodule.py:2053: UserWarning: `pytensor.config.cxx` is not an identifiable `g++` compiler. PyTensor will disable compiler optimizations specific to `g++`. At worst, this could cause slow downs.
Those parameters can be added manually via the `cxxflags` setting.
  warnings.warn(
/Users/benjamv/git/pathmc-worktrees/issue-538/.venv/lib/python3.12/site-packages/pytensor/link/c/cmodule.py:2053: UserWarning: `pytensor.config.cxx` is not an identifiable `g++` compiler. PyTensor will disable compiler optimizations specific to `g++`. At worst, this could cause slow downs.
Those parameters can be added manually via the `cxxflags` setting.
  warnings.warn(
/Users/benjamv/git/pathmc-worktrees/issue-538/.venv/lib/python3.12/site-packages/pytensor/link/c/cmodule.py:2053: UserWarning: `pytensor.config.cxx` is not an identifiable `g++` compiler. PyTensor will disable compiler optimizations specific to `g++`. At worst, this could cause slow downs.
Those parameters can be added manually via the `cxxflags` setting.
  warnings.warn(
/Users/benjamv/git/pathmc-worktrees/issue-538/.venv/lib/python3.12/site-packages/pytensor/link/c/cmodule.py:2053: UserWarning: `pytensor.config.cxx` is not an identifiable `g++` compiler. PyTensor will disable compiler optimizations specific to `g++`. At worst, this could cause slow downs.
Those parameters can be added manually via the `cxxflags` setting.
  warnings.warn(
/Users/benjamv/git/pathmc-worktrees/issue-538/.venv/lib/python3.12/site-packages/pytensor/link/c/cmodule.py:2053: UserWarning: `pytensor.config.cxx` is not an identifiable `g++` compiler. PyTensor will disable compiler optimizations specific to `g++`. At worst, this could cause slow downs.
Those parameters can be added manually via the `cxxflags` setting.
  warnings.warn(

PyMC plate diagram for the panel MMM with adstock and saturation transforms.

Sample

NoteSaturation–coefficient degeneracy

The b * saturation(lam * x) parameterization creates a funnel-shaped degeneracy between the coefficient and saturation steepness. This can cause divergences and max tree depth warnings. Consider tighter priors, reparameterization, or a non-centered saturation formulation if sampling is difficult.

model.fit(random_seed=42)

Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
NUTS[nutpie]: [alpha_sales_offset, sigma_alpha_sales, mu_alpha_sales, lam_dig, theta_dig, lam_tv, theta_tv, beta_sales, sigma_sales]
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.

Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.

<xarray.DataTree>
Group: /
├── Group: /posterior
│       Dimensions:             (chain: 4, draw: 1000, unit: 4, sales_predictors: 3,
│                                mu_sales_dim_0: 30, mu_sales_dim_1: 4)
│       Coordinates:
│         * chain               (chain) int64 32B 0 1 2 3
│         * draw                (draw) int64 8kB 0 1 2 3 4 5 ... 994 995 996 997 998 999
│         * unit                (unit) object 32B 'East' 'North' 'South' 'West'
│         * sales_predictors    (sales_predictors) object 24B 'tv' 'digital' 'trend'
│         * mu_sales_dim_0      (mu_sales_dim_0) int64 240B 0 1 2 3 4 ... 25 26 27 28 29
│         * mu_sales_dim_1      (mu_sales_dim_1) int64 32B 0 1 2 3
│       Data variables:
│           alpha_sales_offset  (chain, draw, unit) float64 128kB -2.29 ... 0.5133
│           mu_alpha_sales      (chain, draw) float64 32kB 20.78 34.41 ... 14.32 16.56
│           beta_sales          (chain, draw, sales_predictors) float64 96kB 24.63 .....
│           sigma_alpha_sales   (chain, draw) float64 32kB 2.784 4.107 ... 3.638 3.438
│           lam_dig             (chain, draw) float64 32kB 0.3309 1.832 ... 1.889 0.2791
│           theta_dig           (chain, draw) float64 32kB 0.5953 0.5657 ... 0.9706
│           lam_tv              (chain, draw) float64 32kB 0.8441 1.519 ... 1.793 0.4877
│           theta_tv            (chain, draw) float64 32kB 0.7207 0.2873 ... 0.3393
│           sigma_sales         (chain, draw) float64 32kB 1.326 1.356 ... 1.281 1.4
│           mu_sales            (chain, draw, mu_sales_dim_0, mu_sales_dim_1) float64 4MB ...
│           alpha_sales         (chain, draw, unit) float64 128kB 14.41 19.6 ... 18.33
│       Attributes:
│           created_at:                 2026-10-30T08:06:21.502587+00:00
│           creation_library:           ArviZ
│           creation_library_version:   1.2.0
│           creation_library_language:  Python
│           sample_dims:                ['chain', 'draw']
│           inference_library:          nutpie
│           inference_library_version:  0.16.10
│           sampling_time:              16.471026182174683
│           tuning_steps:               400
├── Group: /sample_stats
│       Dimensions:                   (chain: 4, draw: 1000)
│       Coordinates:
│         * chain                     (chain) int64 32B 0 1 2 3
│         * draw                      (draw) int64 8kB 0 1 2 3 4 ... 995 996 997 998 999
│       Data variables: (12/20)
│           depth                     (chain, draw) uint64 32kB 5 8 7 8 9 ... 7 7 7 8 6
│           maxdepth_reached          (chain, draw) bool 4kB False False ... False False
│           step_size                 (chain, draw) float64 32kB 0.09302 ... 0.08983
│           transformation_update_id  (chain, draw) int64 32kB 0 0 0 0 0 0 ... 0 0 0 0 0
│           step_size_bar             (chain, draw) float64 32kB 0.1007 ... 0.09101
│           mean_tree_accept          (chain, draw) float64 32kB 0.9958 ... 0.6168
│           ...                        ...
│           fisher_distance           (chain, draw) float64 32kB 208.4 ... 1.126e+03
│           transformation_index      (chain, draw) int64 32kB 339 339 339 ... 339 339
│           diverging                 (chain, draw) bool 4kB False False ... False False
│           divergence_draw           (chain, draw) uint64 32kB 0 0 0 0 0 ... 0 0 0 0 0
│           divergence_message        (chain, draw) object 32kB None None ... None None
│           divergence_energy_error   (chain, draw) float64 32kB nan nan nan ... nan nan
│       Attributes:
│           created_at:                  2026-10-30T08:06:21.495966+00:00
│           creation_library:            ArviZ
│           creation_library_version:    1.2.0
│           creation_library_language:   Python
│           sample_dims:                 ['chain', 'draw']
│           inference_library:           nutpie
│           inference_library_version:   0.16.10
│           inference_library_settings:  {"sampler": "nuts", "adaptation": "diag", "s...
├── Group: /constant_data
│       Dimensions:                    (_do_intervene_sales_dim_0: 30,
│                                       _do_intervene_sales_dim_1: 4, tv_dim_0: 30,
│                                       tv_dim_1: 4, trend_dim_0: 30, trend_dim_1: 4,
│                                       digital_dim_0: 30, digital_dim_1: 4)
│       Coordinates:
│         * _do_intervene_sales_dim_0  (_do_intervene_sales_dim_0) int64 240B 0 1 ... 29
│         * _do_intervene_sales_dim_1  (_do_intervene_sales_dim_1) int64 32B 0 1 2 3
│         * tv_dim_0                   (tv_dim_0) int64 240B 0 1 2 3 4 ... 26 27 28 29
│         * tv_dim_1                   (tv_dim_1) int64 32B 0 1 2 3
│         * trend_dim_0                (trend_dim_0) int64 240B 0 1 2 3 ... 26 27 28 29
│         * trend_dim_1                (trend_dim_1) int64 32B 0 1 2 3
│         * digital_dim_0              (digital_dim_0) int64 240B 0 1 2 3 ... 27 28 29
│         * digital_dim_1              (digital_dim_1) int64 32B 0 1 2 3
│       Data variables:
│           _do_intervene_sales        (_do_intervene_sales_dim_0, _do_intervene_sales_dim_1) float64 960B ...
│           _use_observed_carry        int32 4B 1
│           tv                         (tv_dim_0, tv_dim_1) float64 960B 30.75 ... 14.68
│           trend                      (trend_dim_0, trend_dim_1) float64 960B 1.0 .....
│           digital                    (digital_dim_0, digital_dim_1) float64 960B 12...
│       Attributes:
│           created_at:                 2026-10-30T08:06:21.499991+00:00
│           creation_library:           ArviZ
│           creation_library_version:   1.2.0
│           creation_library_language:  Python
│           inference_library:          pymc
│           inference_library_version:  6.2.0
│           sample_dims:                []
├── Group: /observed_data
│       Dimensions:      (sales_dim_0: 30, sales_dim_1: 4)
│       Coordinates:
│         * sales_dim_0  (sales_dim_0) int64 240B 0 1 2 3 4 5 6 ... 23 24 25 26 27 28 29
│         * sales_dim_1  (sales_dim_1) int64 32B 0 1 2 3
│       Data variables:
│           sales        (sales_dim_0, sales_dim_1) float64 960B 59.36 65.21 ... 72.55
│       Attributes:
│           created_at:                 2026-10-30T08:06:21.501748+00:00
│           creation_library:           ArviZ
│           creation_library_version:   1.2.0
│           creation_library_language:  Python
│           inference_library:          pymc
│           inference_library_version:  6.2.0
│           sample_dims:                []
└── Group: /log_likelihood
        Dimensions:      (chain: 4, draw: 1000, sales_dim_0: 30, sales_dim_1: 4)
        Coordinates:
          * chain        (chain) int64 32B 0 1 2 3
          * draw         (draw) int64 8kB 0 1 2 3 4 5 6 ... 993 994 995 996 997 998 999
          * sales_dim_0  (sales_dim_0) int64 240B 0 1 2 3 4 5 6 ... 23 24 25 26 27 28 29
          * sales_dim_1  (sales_dim_1) int64 32B 0 1 2 3
        Data variables:
            sales        (chain, draw, sales_dim_0, sales_dim_1) float64 4MB -1.546 ....
        Attributes:
            created_at:                 2026-10-30T08:06:23.227509+00:00
            creation_library:           ArviZ
            creation_library_version:   1.2.0
            creation_library_language:  Python
            inference_library:          pymc
            inference_library_version:  6.2.0
            sample_dims:                ['chain', 'draw']

xarray.DataTree

• /posterior(26)

  • Dimensions:
    • chain: 4
    • draw: 1000
    • unit: 4
    • sales_predictors: 3
    • mu_sales_dim_0: 30
    • mu_sales_dim_1: 4
  • Coordinates: (6)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999], shape=(1000,))

    • unit
      (unit)
      object
      'East' 'North' 'South' 'West'

      array(['East', 'North', 'South', 'West'], dtype=object)

    • sales_predictors
      (sales_predictors)
      object
      'tv' 'digital' 'trend'

      array(['tv', 'digital', 'trend'], dtype=object)

    • mu_sales_dim_0
      (mu_sales_dim_0)
      int64
      0 1 2 3 4 5 6 ... 24 25 26 27 28 29

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29])

    • mu_sales_dim_1
      (mu_sales_dim_1)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

  • Data variables: (11)
    • alpha_sales_offset
      (chain, draw, unit)
      float64
      -2.29 -0.425 2.98 ... 1.601 0.5133

      array([[[-2.29008954, -0.42496647,  2.97983217,  1.56183114],[-1.64039053, -0.44816836,  1.88616846,  0.94178383],[-1.42490628,  0.04020685,  2.67091357,  1.47507107],...,[-2.60143075, -0.91008996,  1.84789962,  0.52104548],[-2.42984143, -1.07893867,  1.81733898,  0.79390373],[-2.02538104, -0.57637411,  2.43659527,  1.4117505 ]],[[-0.67835557,  0.59707946,  2.73849828,  1.83744898],[-2.4837976 , -0.56142943,  3.33080617,  1.76334479],[-2.55058392, -0.55984794,  3.0634884 ,  1.47113387],...,[-2.05376612, -0.11630501,  3.51292236,  1.98854766],[-1.98341131, -0.25633961,  2.80663809,  1.51482682],[-1.82476736, -0.22285089,  2.80769952,  1.39874314]],[[-2.5824543 , -1.30098459,  1.38796789,  0.4624621 ],[-1.99490967, -0.31307972,  3.25782947,  1.98220738],[-2.34069645, -0.23269052,  3.84888945,  2.12021549],...,[-1.60017959,  0.16609716,  3.63488501,  2.28520378],[-1.59133595,  0.10988573,  3.58747849,  2.15649626],[-1.81960519,  0.14361998,  3.42810064,  1.93707794]],[[-3.15644302, -1.54377378,  1.74793644,  0.23963159],[-3.39195658, -1.62562879,  1.32965789,  0.18053084],[-2.98738656, -1.74286356,  0.85984612, -0.32204656],...,[-3.17334579, -1.27408926,  2.30780274,  0.85412621],[-2.60278578, -1.2117044 ,  1.48212787,  0.36955101],[-2.66703667, -1.03142198,  1.6011697 ,  0.51329829]]],shape=(4, 1000, 4))

    • mu_alpha_sales
      (chain, draw)
      float64
      20.78 34.41 40.56 ... 14.32 16.56

      array([[20.78083749, 34.40898935, 40.56352722, ..., 28.04056517,27.78427795, 28.15646447],[21.89141231, 24.56536531, 25.26597222, ..., 23.36732623,18.78062355, 19.31217556],[42.06538588, 21.7254304 , 22.82427724, ..., 25.70235044,25.5599965 , 25.20761597],[22.69376971, 30.00004012, 30.80695914, ..., 28.13287952,14.3200506 , 16.56073026]], shape=(4, 1000))

    • beta_sales
      (chain, draw, sales_predictors)
      float64
      24.63 19.65 ... 25.55 0.07918

      array([[[24.63351911, 19.64877388,  0.11280026],[16.17587203, 15.38195848,  0.1107729 ],[ 7.7061759 , 15.58983996,  0.11905205],...,[13.3863884 , 25.94635586,  0.10967116],[12.86323171, 26.20285164,  0.1255306 ],[10.99755898, 26.77109591,  0.08296415]],[[26.21257805, 13.77118328,  0.08575858],[16.50329329, 24.50150295,  0.08256477],[16.50802403, 24.120953  ,  0.08606829],...,[20.18544316, 20.49128384,  0.12210831],[26.47008726, 19.82891949,  0.09548955],[26.2421751 , 19.39863377,  0.09703768]],[[ 9.84603752, 16.18882768,  0.11677591],[13.16291252, 29.6757382 ,  0.09842344],[10.08829949, 31.55992447,  0.09894398],...,[22.9610933 , 14.94158677,  0.10466981],[23.01737297, 15.12477647,  0.10282143],[23.4045029 , 15.38437879,  0.10406987]],[[29.59656159, 16.44277115,  0.0981386 ],[22.44489067, 16.94032016,  0.096765  ],[22.05070726, 17.61868319,  0.10670145],...,[17.62944362, 21.67759845,  0.11523706],[27.15946475, 26.92915465,  0.10820139],[26.53076516, 25.55192152,  0.07917873]]], shape=(4, 1000, 3))

    • sigma_alpha_sales
      (chain, draw)
      float64
      2.784 4.107 3.563 ... 3.638 3.438

      array([[2.78382976, 4.10658528, 3.56293128, ..., 3.29160138, 3.37406924,3.22713083],[4.29392319, 2.48580886, 2.59487929, ..., 2.58044449, 3.09142488,3.16352845],[3.613106  , 2.70458513, 2.4462756 , ..., 2.75482449, 2.82427258,2.84222497],[3.03167915, 3.06836032, 3.81293674, ..., 2.71989762, 3.63760481,3.43800066]], shape=(4, 1000))

    • lam_dig
      (chain, draw)
      float64
      0.3309 1.832 ... 1.889 0.2791

      array([[0.33089651, 1.83169046, 0.75716335, ..., 1.76859841, 0.67979643,1.60348055],[0.36973544, 0.75954567, 1.69142263, ..., 1.83113964, 0.55971376,0.51428214],[0.43478478, 0.42416226, 2.65061823, ..., 0.53205138, 0.41210545,0.59972953],[0.69744851, 1.35234774, 0.54045545, ..., 2.51395045, 1.88850645,0.27909103]], shape=(4, 1000))

    • theta_dig
      (chain, draw)
      float64
      0.5953 0.5657 ... 0.03791 0.9706

      array([[0.59527858, 0.56571855, 0.4874882 , ..., 0.47499008, 0.46824888,0.47644364],[0.60615502, 0.63708566, 0.6990028 , ..., 0.82892722, 0.42284641,0.59381502],[0.7286007 , 0.4120705 , 0.4072078 , ..., 0.38659099, 0.43281879,0.53958799],[0.72265434, 0.71410729, 0.28984627, ..., 0.52072899, 0.0379137 ,0.97061575]], shape=(4, 1000))

    • lam_tv
      (chain, draw)
      float64
      0.8441 1.519 1.646 ... 1.793 0.4877

      array([[0.84414467, 1.51904908, 1.64617108, ..., 1.34657431, 0.66963427,1.25744293],[2.21576058, 0.29039246, 2.44260344, ..., 0.27530268, 0.55836903,2.04115136],[1.09775129, 1.42850834, 0.76125802, ..., 1.20736596, 0.94888772,0.49976095],[1.56253792, 1.22031723, 0.70305584, ..., 1.15717107, 1.79254957,0.48768832]], shape=(4, 1000))

    • theta_tv
      (chain, draw)
      float64
      0.7207 0.2873 ... 0.6487 0.3393

      array([[0.72065866, 0.28725197, 0.27119088, ..., 0.06378632, 0.5274306 ,0.39477409],[0.67067366, 0.57727616, 0.51704411, ..., 0.71367915, 0.13884074,0.09650477],[0.75721581, 0.37557579, 0.60533951, ..., 0.51539441, 0.58800586,0.60829306],[0.3889351 , 0.46350697, 0.63815698, ..., 0.29936229, 0.64874407,0.33930302]], shape=(4, 1000))

    • sigma_sales
      (chain, draw)
      float64
      1.326 1.356 1.334 ... 1.281 1.4

      array([[1.32606735, 1.35646222, 1.33394038, ..., 1.30799278, 1.31793789,1.32091782],[1.38232012, 1.25232265, 1.4680443 , ..., 1.44279146, 1.29946715,1.2736775 ],[1.35327651, 1.39238312, 1.25410788, ..., 1.28447323, 1.24459205,1.21304929],[1.29197034, 1.29119082, 1.31865933, ..., 1.43696629, 1.28139939,1.39953237]], shape=(4, 1000))

    • mu_sales
      (chain, draw, mu_sales_dim_0, mu_sales_dim_1)
      float64
      58.26 63.79 73.15 ... 76.52 72.78

      array([[[[58.25775633, 63.79479141, 73.15255193, 69.23641816],[58.88661329, 63.9579488 , 73.56254214, 69.63631282],[59.01511462, 64.21809025, 73.69684661, 69.74792766],...,[61.84129264, 67.03692296, 76.51600301, 72.5694069 ],[61.92877712, 67.1509058 , 76.62968009, 72.68220992],[62.00418629, 67.26183681, 76.74248355, 72.794969  ]],[[59.34118916, 64.23715115, 73.82330441, 69.9451084 ],[59.45196206, 64.34792406, 73.93407731, 70.0558813 ],[59.56273497, 64.45869696, 74.04485021, 70.1666542 ],...,[62.33205754, 67.22801954, 76.81417279, 72.93597678],[62.44283044, 67.33879244, 76.92494569, 73.04674968],[62.55360335, 67.44956534, 77.0357186 , 73.15752258]],[[58.89996342, 64.12167486, 73.49435496, 69.23376106],[59.02079888, 64.24057395, 73.61392517, 69.35322403],[59.13985364, 64.35995346, 73.73298075, 69.472276  ],...,...[62.0353836 , 67.20116691, 76.94354643, 72.9896951 ],[62.15062066, 67.31640397, 77.05878349, 73.10493216],[62.26585772, 67.43164103, 77.17402055, 73.22016922]],[[59.04896534, 64.10916963, 73.90826686, 69.86115193],[59.15716673, 64.21735507, 74.01646826, 69.96935333],[59.26536812, 64.32557242, 74.12466965, 70.07755461],...,[61.97040294, 67.03013904, 76.82970398, 72.78258954],[62.07795952, 67.13880863, 76.93790586, 72.89079093],[62.18658224, 67.24700905, 77.04610725, 72.99898396]],[[58.19380487, 64.59106924, 73.35563408, 69.68786459],[59.60415445, 65.16919544, 74.28678103, 70.56598722],[59.71027633, 65.33481618, 74.38574313, 70.64547823],...,[61.69112708, 67.31438905, 76.36524397, 72.62514134],[61.77004279, 67.39357075, 76.4444227 , 72.70432008],[61.84928714, 67.47274948, 76.52360143, 72.78349285]]]],shape=(4, 1000, 30, 4))

    • alpha_sales
      (chain, draw, unit)
      float64
      14.41 19.6 29.08 ... 22.07 18.33

      array([[[14.40561808, 19.59780318, 29.07618295, 25.12870948],[27.67258575, 32.56854774, 42.154701  , 38.27650499],[35.48668406, 40.70678145, 50.07980874, 45.81910407],...,[19.47769212, 25.0449118 , 34.12311413, 29.7556392 ],[19.58582472, 24.14386417, 33.91610552, 30.46296412],[21.62029487, 26.2964298 , 36.01967616, 32.71236804]],[[18.9786056 , 24.45522563, 33.65031357, 29.78127711],[18.39111923, 23.16975906, 32.8451128 , 28.9487034 ],[18.64751482, 23.81323438, 33.21535484, 29.08338704],...,[18.06769675, 23.0672076 , 32.4322274 , 28.49866308],[12.64905649, 17.98816891, 27.45713437, 23.46359688],[13.53947209, 18.60718042, 28.19441288, 23.73713928]],[[32.73470476, 37.36479064, 47.08026098, 43.73631048],[16.33002736, 20.87867965, 30.53650756, 27.08647901],[17.09828862, 22.25505209, 32.2397216 , 28.01090865],...,[21.29413651, 26.15991896, 35.7158207 , 31.99768578],[21.06563001, 25.87034377, 35.69201363, 31.65052976],[20.03588866, 25.61581626, 34.9510492 , 30.71322725]],[[13.1244472 , 18.01354294, 27.99295216, 23.4202558 ],[19.59229514, 25.01202524, 34.07990962, 30.55397378],[19.41624315, 24.16153065, 34.08549799, 29.57901598],...,[19.50170387, 24.66748718, 34.4098667 , 30.45601537],[ 4.85214455,  9.91234884, 19.71144608, 15.66433115],[ 7.39145642, 13.01470081, 22.06555276, 18.32545014]]],shape=(4, 1000, 4))

  • Attributes: (9)

    created_at :

    2026-10-30T08:06:21.502587+00:00

    creation_library :

    ArviZ

    creation_library_version :

    1.2.0

    creation_library_language :

    Python

    sample_dims :

    ['chain', 'draw']

    inference_library :

    nutpie

    inference_library_version :

    0.16.10

    sampling_time :

    16.471026182174683

    tuning_steps :

    400

  /sample_stats(30)

  • Dimensions:
    • chain: 4
    • draw: 1000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999], shape=(1000,))

  • Data variables: (20)
    • depth
      (chain, draw)
      uint64
      5 8 7 8 9 9 8 7 ... 5 6 5 7 7 7 8 6

      array([[5, 8, 7, ..., 7, 4, 6],[5, 7, 6, ..., 7, 7, 5],[5, 8, 5, ..., 8, 5, 6],[6, 7, 6, ..., 7, 8, 6]], shape=(4, 1000), dtype=uint64)

    • maxdepth_reached
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],[False, False, False, ..., False, False, False],[False, False, False, ..., False, False, False],[False, False, False, ..., False, False, False]], shape=(4, 1000))

    • step_size
      (chain, draw)
      float64
      0.09302 0.09667 ... 0.09357 0.08983

      array([[0.09302007, 0.09666563, 0.10418245, ..., 0.10276652, 0.09778782,0.09196321],[0.09531171, 0.08496788, 0.09745103, ..., 0.09966393, 0.09773026,0.09291642],[0.10200083, 0.10869376, 0.09932861, ..., 0.09676462, 0.10446596,0.09752979],[0.09125252, 0.09198268, 0.09611871, ..., 0.09094499, 0.09357339,0.0898342 ]], shape=(4, 1000))

    • transformation_update_id
      (chain, draw)
      int64
      0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0

      array([[0, 0, 0, ..., 0, 0, 0],[0, 0, 0, ..., 0, 0, 0],[0, 0, 0, ..., 0, 0, 0],[0, 0, 0, ..., 0, 0, 0]], shape=(4, 1000))

    • step_size_bar
      (chain, draw)
      float64
      0.1007 0.1007 ... 0.09101 0.09101

      array([[0.10069426, 0.10069426, 0.10069426, ..., 0.10069426, 0.10069426,0.10069426],[0.09349036, 0.09349036, 0.09349036, ..., 0.09349036, 0.09349036,0.09349036],[0.10025469, 0.10025469, 0.10025469, ..., 0.10025469, 0.10025469,0.10025469],[0.09100975, 0.09100975, 0.09100975, ..., 0.09100975, 0.09100975,0.09100975]], shape=(4, 1000))

    • mean_tree_accept
      (chain, draw)
      float64
      0.9958 0.7119 ... 0.836 0.6168

      array([[0.99578089, 0.71190701, 0.99232313, ..., 0.96811771, 0.86306515,0.98425555],[0.76705521, 0.97291857, 0.99881881, ..., 0.88222119, 0.7068945 ,0.92285033],[0.86976432, 0.85837611, 0.99179291, ..., 0.24494007, 0.95605042,0.93514992],[0.91163586, 0.98954209, 0.98578163, ..., 1.        , 0.83600028,0.61684962]], shape=(4, 1000))

    • mean_tree_accept_sym
      (chain, draw)
      float64
      0.9238 0.7964 ... 0.836 0.7274

      array([[0.92383447, 0.79641249, 0.78471684, ..., 0.98173314, 0.91771783,0.9302765 ],[0.8550203 , 0.91229257, 0.98872642, ..., 0.93059874, 0.75495708,0.95861866],[0.92418218, 0.88728456, 0.94380159, ..., 0.34995664, 0.92281933,0.91992611],[0.94508398, 0.99269873, 0.99139373, ..., 0.81293971, 0.83603451,0.72743316]], shape=(4, 1000))

    • n_steps
      (chain, draw)
      uint64
      31 319 191 511 ... 191 191 271 95

      array([[ 31, 319, 191, ..., 239,  31,  95],[ 31, 255,  63, ..., 127, 255,  31],[ 47, 255,  63, ..., 255,  63, 127],[ 63, 255,  63, ..., 191, 271,  95]], shape=(4, 1000), dtype=uint64)

    • max_energy_error
      (chain, draw)
      float64
      0.3172 -0.9743 ... -0.7199 -1.284

      array([[ 0.3172272 , -0.97432192,  1.04914684, ..., -0.10022436,-0.35870934,  0.33182402],[-0.63430244,  0.40470771,  0.06682697, ..., -0.32562633,-1.51641276, -0.16650896],[-0.35143074, -0.58889615,  0.3042587 , ..., -3.59809921,0.36445798, -0.33148085],[-0.24186276, -0.04228073, -0.04739058, ...,  0.78242255,-0.71986568, -1.28381087]], shape=(4, 1000))

    • tuning
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],[False, False, False, ..., False, False, False],[False, False, False, ..., False, False, False],[False, False, False, ..., False, False, False]], shape=(4, 1000))

    • index_in_trajectory
      (chain, draw)
      int64
      -16 -99 68 -186 ... 69 39 -100 -31

      array([[ -16,  -99,   68, ...,  108,   13,  -28],[ -13,  102,  -20, ...,   46,   38,   -8],[  15,  144,  -22, ...,  125,    4,    6],[ -43,  -69,   36, ...,   39, -100,  -31]], shape=(4, 1000))

    • logp
      (chain, draw)
      float64
      -243.0 -245.7 ... -247.6 -253.6

      array([[-243.01306021, -245.69345817, -245.34008302, ..., -247.32348102,-245.6849171 , -246.61165937],[-247.94687886, -246.02139871, -246.96981238, ..., -247.57188746,-243.09035141, -245.22254802],[-249.43790193, -246.86451266, -249.3542421 , ..., -244.67021641,-244.50724307, -244.7816169 ],[-244.83924552, -244.955936  , -246.37950787, ..., -246.76212296,-247.58006267, -253.57003818]], shape=(4, 1000))

    • energy
      (chain, draw)
      float64
      228.5 235.5 234.5 ... 241.2 238.6

      array([[228.47609509, 235.51126192, 234.47690317, ..., 236.41434546,233.12203127, 232.76369832],[238.40149395, 236.98470432, 233.79259681, ..., 235.33935941,238.65255704, 233.87461301],[243.04509124, 235.66117237, 236.91294697, ..., 236.3496099 ,232.27734182, 233.74768646],[235.70811244, 233.04533093, 231.47708691, ..., 235.11495152,241.16436076, 238.63831273]], shape=(4, 1000))

    • energy_error
      (chain, draw)
      float64
      -0.1138 0.6719 ... -0.3006 1.101

      array([[-0.11383958,  0.67191933, -0.91069441, ...,  0.04659589,0.27949516, -0.24676583],[ 0.32039483, -0.31456287, -0.03113501, ...,  0.13504815,-0.34645177, -0.0013487 ],[ 0.15657585,  0.04175885, -0.05831063, ...,  0.57400379,-0.07257658, -0.27149622],[-0.08719671, -0.01045752,  0.04451152, ..., -0.23531959,-0.30061949,  1.10072565]], shape=(4, 1000))

    • fisher_distance
      (chain, draw)
      float64
      208.4 1.041e+03 ... 434.3 1.126e+03

      array([[ 208.4482042 , 1041.02454315,  683.22480749, ...,  183.1063662 ,403.67998277,  211.2903906 ],[ 508.63786967,  156.61524029,   84.53952115, ...,  148.27975099,62.16450034,   68.07926862],[ 833.89019335,  431.48454621,  605.85343369, ...,  382.67548454,335.63641053,  146.0900792 ],[ 177.10669614,  196.55764679,  301.66772766, ...,  391.10611813,434.26101053, 1126.2685612 ]], shape=(4, 1000))

    • transformation_index
      (chain, draw)
      int64
      339 339 339 339 ... 339 339 339 339

      array([[339, 339, 339, ..., 339, 339, 339],[339, 339, 339, ..., 339, 339, 339],[337, 337, 337, ..., 337, 337, 337],[339, 339, 339, ..., 339, 339, 339]], shape=(4, 1000))

    • diverging
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],[False, False, False, ..., False, False, False],[False, False, False, ..., False, False, False],[False, False, False, ..., False, False, False]], shape=(4, 1000))

    • divergence_draw
      (chain, draw)
      uint64
      0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0

      array([[0, 0, 0, ..., 0, 0, 0],[0, 0, 0, ..., 0, 0, 0],[0, 0, 0, ..., 0, 0, 0],[0, 0, 0, ..., 0, 0, 0]], shape=(4, 1000), dtype=uint64)

    • divergence_message
      (chain, draw)
      object
      None None None ... None None None

      array([[None, None, None, ..., None, None, None],[None, None, None, ..., None, None, None],[None, None, None, ..., None, None, None],[None, None, None, ..., None, None, None]],shape=(4, 1000), dtype=object)

    • divergence_energy_error
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],[nan, nan, nan, ..., nan, nan, nan],[nan, nan, nan, ..., nan, nan, nan],[nan, nan, nan, ..., nan, nan, nan]], shape=(4, 1000))

  • Attributes: (8)

    created_at :

    2026-10-30T08:06:21.495966+00:00

    creation_library :

    ArviZ

    creation_library_version :

    1.2.0

    creation_library_language :

    Python

    sample_dims :

    ['chain', 'draw']

    inference_library :

    nutpie

    inference_library_version :

    0.16.10

    inference_library_settings :

    {"sampler": "nuts", "adaptation": "diag", "settings": {"num_tune": 400, "num_draws": 1000, "maxdepth": 10, "mindepth": 0, "store_gradient": false, "store_unconstrained": false, "store_transformed": false, "max_energy_error": 1000.0, "store_divergences": false, "adapt_options": {"step_size_settings": {"target_accept": 0.8, "initial_step": 0.1, "jitter": 0.1, "adapt_options": {"method": "DualAverage", "dual_average": {"k": 0.75, "t0": 10.0, "gamma": 0.05, "max_step_size": 3.141592653589793}, "adam": {"beta1": 0.9, "beta2": 0.999, "epsilon": 1e-08, "learning_rate": 0.05}}}, "mass_matrix_options": {"store_mass_matrix": false, "use_grad_based_estimate": true}, "early_window": 0.3, "step_size_window": 0.15, "mass_matrix_switch_freq": 80, "early_mass_matrix_switch_freq": 10, "mass_matrix_update_freq": 1, "mass_matrix_window_growth": 1.5}, "check_turning": true, "target_integration_time": null, "trajectory_kind": "Euclidean", "num_chains": 4, "seed": 534082709, "extra_doublings": 0}}

  /constant_data(20)

  • Dimensions:
    • _do_intervene_sales_dim_0: 30
    • _do_intervene_sales_dim_1: 4
    • tv_dim_0: 30
    • tv_dim_1: 4
    • trend_dim_0: 30
    • trend_dim_1: 4
    • digital_dim_0: 30
    • digital_dim_1: 4
  • Coordinates: (8)
    • _do_intervene_sales_dim_0
      (_do_intervene_sales_dim_0)
      int64
      0 1 2 3 4 5 6 ... 24 25 26 27 28 29

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29])

    • _do_intervene_sales_dim_1
      (_do_intervene_sales_dim_1)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • tv_dim_0
      (tv_dim_0)
      int64
      0 1 2 3 4 5 6 ... 24 25 26 27 28 29

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29])

    • tv_dim_1
      (tv_dim_1)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • trend_dim_0
      (trend_dim_0)
      int64
      0 1 2 3 4 5 6 ... 24 25 26 27 28 29

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29])

    • trend_dim_1
      (trend_dim_1)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • digital_dim_0
      (digital_dim_0)
      int64
      0 1 2 3 4 5 6 ... 24 25 26 27 28 29

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29])

    • digital_dim_1
      (digital_dim_1)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

  • Data variables: (5)
    • _do_intervene_sales
      (_do_intervene_sales_dim_0, _do_intervene_sales_dim_1)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan],[nan, nan, nan, nan]])

    • _use_observed_carry
      ()
      int32
      1

      array(1, dtype=int32)

    • tv
      (tv_dim_0, tv_dim_1)
      float64
      30.75 40.96 27.85 ... 39.94 14.68

      array([[30.75433442, 40.95824194, 27.84625102, 37.22858107],[36.46645053, 37.89472116, 35.21130372, 30.18104948],[39.87158445, 40.44558808, 14.72023608, 11.73900248],[19.638823  , 28.01543752, 37.98828535, 19.97550303],[16.13137266, 35.7546048 , 41.15003616, 11.9701648 ],[44.98248163, 19.08954887, 20.88966247, 14.06687612],[41.09619353, 43.10524688, 28.23105159, 46.99367476],[15.7559001 , 24.18103873, 33.16878276, 46.10612576],[15.27288871, 41.13533988, 40.34078119, 41.17910163],[30.25319726, 11.75215063, 35.09235363, 15.38208834],[17.99102607, 39.79048624, 13.37777285, 44.30288576],[39.55936016, 24.81838824, 29.75963277, 35.58253085],[22.03794767, 15.19686021, 14.13611871, 29.11509096],[48.22493028, 36.79255979, 47.00480474, 24.7004724 ],[10.99437966, 43.29039206, 33.63661966, 25.17856319],[14.23589615, 21.53312416, 29.29213748, 47.95431707],[48.64927649, 17.9963281 , 29.46633323, 23.1344482 ],[42.17443662, 36.59403426, 32.8691221 , 36.10349362],[10.71347136, 28.35663102, 23.26275989, 35.31670517],[41.87268353, 14.58120294, 10.8644832 , 18.09618373],[34.24063283, 32.60944426, 15.60996356, 34.57518962],[26.50286277, 32.14317603, 36.88960372, 13.37972875],[36.07724102, 11.23271338, 39.07978457, 48.35520984],[19.91358252, 26.34114575, 46.64047381, 41.30496419],[42.66275054, 12.33210967, 32.19409878, 20.13533415],[33.77734655, 36.47666059, 42.33005888, 34.646628  ],[22.41338696, 36.57254161, 21.63671353, 32.89722056],[16.62126891, 16.6789168 , 47.4417428 , 40.5569356 ],[14.61960255, 38.89437402, 27.4038824 , 33.98373766],[16.98565884, 30.041791  , 39.94432078, 14.68414264]])

    • trend
      (trend_dim_0, trend_dim_1)
      float64
      1.0 1.0 1.0 1.0 ... 30.0 30.0 30.0

      array([[ 1.,  1.,  1.,  1.],[ 2.,  2.,  2.,  2.],[ 3.,  3.,  3.,  3.],[ 4.,  4.,  4.,  4.],[ 5.,  5.,  5.,  5.],[ 6.,  6.,  6.,  6.],[ 7.,  7.,  7.,  7.],[ 8.,  8.,  8.,  8.],[ 9.,  9.,  9.,  9.],[10., 10., 10., 10.],[11., 11., 11., 11.],[12., 12., 12., 12.],[13., 13., 13., 13.],[14., 14., 14., 14.],[15., 15., 15., 15.],[16., 16., 16., 16.],[17., 17., 17., 17.],[18., 18., 18., 18.],[19., 19., 19., 19.],[20., 20., 20., 20.],[21., 21., 21., 21.],[22., 22., 22., 22.],[23., 23., 23., 23.],[24., 24., 24., 24.],[25., 25., 25., 25.],[26., 26., 26., 26.],[27., 27., 27., 27.],[28., 28., 28., 28.],[29., 29., 29., 29.],[30., 30., 30., 30.]])

    • digital
      (digital_dim_0, digital_dim_1)
      float64
      12.9 15.97 14.53 ... 27.27 7.193

      array([[12.8982263 , 15.97196099, 14.52553065, 14.84076142],[14.34144322,  7.3544337 , 14.04531526, 26.87512356],[11.5615129 , 24.65160763, 29.04744161,  9.53746024],[ 8.06894831, 14.26995061, 11.64674904, 19.28081629],[ 9.4817077 , 25.56904033, 22.92225473, 14.34035346],[ 9.91086664, 18.86461968,  7.40977405, 25.83646384],[29.29566065, 20.79160998, 10.05908412,  7.47782854],[ 5.34840719, 29.26745061,  9.41931957, 29.48926701],[21.94146684,  9.8659677 , 22.9865739 , 21.0620819 ],[22.35656089,  8.8572373 , 19.60244922, 18.4017009 ],[25.10311315, 29.18774331, 15.39518505, 16.56998414],[ 8.27644379, 16.73889528, 13.24653031, 11.66158294],[17.21460113, 16.89262316, 19.6911143 , 15.42223421],[12.16115567, 15.92879797, 19.52652849, 14.15981125],[18.87995106, 25.11910894,  5.57009678, 22.14358364],[ 8.50848993, 22.0623876 , 24.56838068, 27.90870049],[19.90106383,  5.18405674, 17.26767486, 18.38586975],[16.68454004, 22.62913447, 16.83723503, 25.1097957 ],[ 7.72859992, 19.2185299 , 18.01681006, 12.20389035],[10.81601855, 21.71007404, 25.6572981 , 22.36995322],[26.69347384, 24.12497144, 18.85090359,  7.37739371],[14.35460109, 18.98017902, 12.0308446 , 28.39849557],[26.68726579, 15.91793473, 24.2161873 , 25.0221044 ],[10.91655907, 26.33507683, 10.75534977, 24.87787097],[ 7.631952  , 12.0345973 , 14.2730571 , 19.75189738],[ 8.6543311 , 18.92580381, 12.92847232,  9.2822826 ],[ 8.59679832, 15.15967154, 17.87642823, 16.64962882],[12.11800206,  5.56780183,  9.11519544, 24.98111791],[ 5.52870041, 16.54693076, 29.8093891 , 28.28090589],[ 6.33454832,  8.80780257, 27.26981227,  7.1927253 ]])

  • Attributes: (7)

    created_at :

    2026-10-30T08:06:21.499991+00:00

    creation_library :

    ArviZ

    creation_library_version :

    1.2.0

    creation_library_language :

    Python

    inference_library :

    pymc

    inference_library_version :

    6.2.0

    sample_dims :

    []

  /observed_data(10)

  • Dimensions:
    • sales_dim_0: 30
    • sales_dim_1: 4
  • Coordinates: (2)
    • sales_dim_0
      (sales_dim_0)
      int64
      0 1 2 3 4 5 6 ... 24 25 26 27 28 29

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29])

    • sales_dim_1
      (sales_dim_1)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

  • Data variables: (1)
    • sales
      (sales_dim_0, sales_dim_1)
      float64
      59.36 65.21 73.43 ... 76.82 72.55

      array([[59.35894385, 65.2054733 , 73.42543162, 69.59579108],[59.84020072, 62.22082831, 72.02694502, 68.11764897],[57.4614564 , 64.27476914, 74.89965784, 73.52084465],[61.97960871, 65.56662585, 76.34880194, 69.97874314],[61.68617691, 65.20126303, 73.09957118, 69.89167614],[58.25282183, 63.16167428, 74.09117934, 69.99411774],[59.69231519, 64.42270531, 75.35061538, 71.09099747],[61.91064425, 64.56820575, 72.62835739, 69.09271246],[56.82472899, 65.69845209, 75.24899653, 69.46174278],[58.17177879, 65.64588891, 74.64096081, 70.22884231],[61.47385362, 64.33163536, 75.42997745, 69.74471907],[59.47369041, 66.89345528, 75.735234  , 69.64205138],[61.10714679, 64.06327101, 74.3403615 , 72.89096179],[59.35594884, 66.21471555, 73.33496608, 69.90711638],[60.76483873, 65.67502812, 73.29298144, 69.99831092],[60.94181067, 65.93539293, 76.85718454, 73.1984673 ],[59.42010709, 66.13340589, 76.32160616, 73.31590106],[59.91391637, 65.32048646, 76.20155555, 70.35514344],[59.8056176 , 65.48728251, 74.98184604, 72.38436649],[60.75474199, 67.45241634, 76.23557249, 70.71402203],[64.41656011, 66.34412935, 73.58797447, 73.83075631],[62.88577238, 67.39001999, 76.6210006 , 71.09956965],[59.09943179, 67.5869583 , 77.7741081 , 73.56684253],[59.75479575, 64.70006896, 74.95227413, 69.87942389],[59.7418178 , 66.7136255 , 75.26783819, 70.49328928],[64.18991166, 66.83780308, 77.99718168, 73.33292123],[62.85978178, 67.38514238, 76.63881997, 72.15651972],[61.59711986, 66.22570305, 76.11698068, 74.22431365],[61.48286605, 66.19570635, 77.15888134, 72.78635944],[62.67097155, 67.66838745, 76.82204246, 72.55249613]])

  • Attributes: (7)

    created_at :

    2026-10-30T08:06:21.501748+00:00

    creation_library :

    ArviZ

    creation_library_version :

    1.2.0

    creation_library_language :

    Python

    inference_library :

    pymc

    inference_library_version :

    6.2.0

    sample_dims :

    []

  /log_likelihood(12)

  • Dimensions:
    • chain: 4
    • draw: 1000
    • sales_dim_0: 30
    • sales_dim_1: 4
  • Coordinates: (4)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999], shape=(1000,))

    • sales_dim_0
      (sales_dim_0)
      int64
      0 1 2 3 4 5 6 ... 24 25 26 27 28 29

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29])

    • sales_dim_1
      (sales_dim_1)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

  • Data variables: (1)
    • sales
      (chain, draw, sales_dim_0, sales_dim_1)
      float64
      -1.546 -1.767 ... -1.278 -1.269

      array([[[[-1.54595157, -1.76700062, -1.22232916, -1.23787849],[-1.45971534, -2.05917881, -1.87164793, -1.85694223],[-1.88751281, -1.20206966, -1.61252732, -5.24871663],...,[-1.2181087 , -1.38827442, -1.24642859, -1.9798844 ],[-1.25769355, -1.46059026, -1.28078692, -1.20424049],[-1.32757477, -1.248153  , -1.20295598, -1.21787348]],[[-1.2239042 , -1.47861555, -1.26683583, -1.25697703],[-1.2647778 , -2.45331894, -2.21218099, -2.24467872],[-2.42365445, -1.23301137, -1.42237836, -4.28106206],...,[-1.37059473, -1.4968196 , -1.35590537, -1.67485689],[-1.47423614, -1.57888731, -1.23868978, -1.24224339],[-1.22756184, -1.23683032, -1.23622553, -1.3232911 ]],[[-1.26627085, -1.53713678, -1.20841063, -1.24390448],[-1.39574065, -2.35335528, -1.91476105, -1.63605419],[-1.99864096, -1.20911478, -1.58954678, -5.81282387],...,...[-1.32798285, -1.51188156, -1.4469093 , -1.65057131],[-1.38944455, -1.58559882, -1.28389888, -1.30604771],[-1.32121296, -1.29504467, -1.31147177, -1.38941819]],[[-1.1961506 , -1.53287581, -1.23788157, -1.18833375],[-1.3089558 , -2.3807028 , -2.37220198, -2.21099653],[-2.15779518, -1.16767722, -1.34978179, -4.77723203],...,[-1.20932173, -1.3639447 , -1.3215744 , -1.79983531],[-1.27472917, -1.43773488, -1.18176054, -1.17021225],[-1.23833929, -1.22095993, -1.1821792 , -1.2275957 ]],[[-1.60162164, -1.35144013, -1.2563203 , -1.25724077],[-1.26929993, -3.47412742, -2.55871753, -2.78527227],[-2.54603873, -1.54192638, -1.3224963 , -3.36560169],...,[-1.25733263, -1.55763521, -1.27081033, -1.90789846],[-1.27612914, -1.62136194, -1.38538086, -1.25679479],[-1.42742785, -1.26484705, -1.27781307, -1.26869791]]]],shape=(4, 1000, 30, 4))

  • Attributes: (7)

    created_at :

    2026-10-30T08:06:23.227509+00:00

    creation_library :

    ArviZ

    creation_library_version :

    1.2.0

    creation_library_language :

    Python

    inference_library :

    pymc

    inference_library_version :

    6.2.0

    sample_dims :

    ['chain', 'draw']

Results

model.summary()

	mean	sd	eti89_lb	eti89_ub	ess_bulk	ess_tail	r_hat	mcse_mean	mcse_sd
alpha_sales_offset[East]	-2.233316	0.624859	-3.234351	-1.258798	2049.373678	2089.746306	1.002678	0.013823	0.009554
alpha_sales_offset[North]	-0.619756	0.512710	-1.429654	0.190414	2046.396929	2007.077572	1.002609	0.011340	0.007999
alpha_sales_offset[South]	2.405143	0.589116	1.469266	3.339014	2230.468168	2099.418540	1.002027	0.012473	0.008967
alpha_sales_offset[West]	1.136707	0.511378	0.328093	1.950571	2151.456992	2060.438335	1.002584	0.011024	0.007901
mu_alpha_sales	21.865802	7.950126	9.357160	34.711160	929.508089	1069.596424	1.003145	0.261223	0.175237
...	...	...	...	...	...	...	...	...	...
mu_sales[29, 3]	72.767321	0.322305	72.256068	73.288583	4133.499727	3578.561962	0.999883	0.005010	0.003617
alpha_sales[East]	14.863741	8.016064	2.356212	27.871525	911.368172	1139.444051	1.003287	0.266227	0.178905
alpha_sales[North]	19.937892	8.017359	7.347588	32.900286	904.055903	1102.438498	1.003383	0.267227	0.179910
alpha_sales[South]	29.450793	8.021985	16.835363	42.330324	908.812177	1100.964943	1.003312	0.266895	0.179836
alpha_sales[West]	25.461902	8.016374	12.866825	38.438535	909.697773	1138.353640	1.003344	0.266466	0.179043

138 rows × 9 columns

model.effects_summary()

	mean	sd	hdi_3%	hdi_97%
name
b_tv	22.450390	8.070091	7.52803	37.758076
b_dig	21.644288	7.955049	6.12319	35.743989

The transform parameters should be close to the true values: theta_tv ≈ 0.7, theta_dig ≈ 0.5, lam_tv ≈ 0.3, lam_dig ≈ 0.4.

Posterior predictive check

The .predict() method runs posterior predictive sampling — a key diagnostic to verify that the model can reproduce the observed data patterns.

idata = model.predict()

Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Sampling: [sales]

pp_var = "sales_obs" if "sales_obs" in idata.posterior_predictive else "sales"
pp = idata.posterior_predictive[pp_var]
pp_mean_2d = pp.mean(dim=("chain", "draw")).values
sorted_regions = sorted(regions)

fig, axes = plt.subplots(
    2, 2, figsize=(FIG_WIDTH, FIG_HEIGHT * 1.5), sharex=True, sharey=True
)
for ax, region in zip(axes.flat, regions):
    obs = df.loc[df["region"] == region, "sales"].values
    weeks = df.loc[df["region"] == region, "week"].values
    r_idx = sorted_regions.index(region)
    pp_region = pp_mean_2d[:, r_idx]

    ax.scatter(
        weeks,
        obs,
        s=20,
        alpha=0.6,
        color=REGION_COLORS[region],
        zorder=3,
        label="Observed",
    )
    ax.plot(
        weeks,
        pp_region,
        "-",
        alpha=0.8,
        color=REGION_COLORS[region],
        lw=2,
        zorder=4,
        label="Predicted mean",
    )
    ax.set_title(region, fontweight="bold")
    ax.legend(fontsize=7, loc="lower right")

fig.supxlabel("Week")
fig.supylabel("Sales")
plt.tight_layout()
plt.show()

Figure 3: Posterior predictive check by region: observed sales (points) vs posterior predictive mean (line).

Counterfactual: what if we doubled TV spend?

Using time-forward do(), we simulate what sales would have been under a different spend scenario. The transforms (adstock + saturation) are automatically recomputed under the intervention.

r_baseline = model.do(
    set={"tv": 30.0, "digital": 15.0}, simulate_over="time", kind="mean"
)
r_double_tv = model.do(
    set={"tv": 60.0, "digital": 15.0}, simulate_over="time", kind="mean"
)
contrast = r_double_tv - r_baseline

/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/2237643227.py:4: UserWarning: Intervention value 60.00 for 'tv' is outside the observed data range [10.71, 48.65]. Results are extrapolations and should be interpreted with caution.
  r_double_tv = model.do(

fig, ax = plt.subplots(figsize=(FIG_WIDTH, FIG_HEIGHT * 1.2))

weeks = np.arange(1, n_weeks + 1)
avg_obs = df.groupby("week")["sales"].mean().values

pp_var = "sales_obs" if "sales_obs" in idata.posterior_predictive else "sales"
pp = idata.posterior_predictive[pp_var]
pp_mean_2d = pp.mean(dim=("chain", "draw")).values
pp_avg = pp_mean_2d.mean(axis=1) if pp_mean_2d.ndim == 2 else pp_mean_2d

ax.scatter(
    weeks,
    avg_obs,
    s=30,
    color="black",
    alpha=0.5,
    zorder=5,
    label="Observed (region avg)",
)
ax.plot(weeks, pp_avg, "-", color="black", lw=2, zorder=4, label="Model fit (PPC mean)")

bl_draws = r_baseline.draws("sales")
bl_mean = bl_draws.mean()
bl_hdi = r_baseline.hdi("sales", prob=0.94)
ax.axhspan(bl_hdi[0], bl_hdi[1], alpha=0.15, color=COLOR_BASELINE, zorder=1)
ax.axhline(
    bl_mean,
    color=COLOR_BASELINE,
    ls="--",
    lw=1.5,
    zorder=2,
    label=f"do(TV=30) = {bl_mean:.1f}",
)

dbl_draws = r_double_tv.draws("sales")
dbl_mean = dbl_draws.mean()
dbl_hdi = r_double_tv.hdi("sales", prob=0.94)
ax.axhspan(dbl_hdi[0], dbl_hdi[1], alpha=0.15, color=COLOR_TV, zorder=1)
ax.axhline(
    dbl_mean,
    color=COLOR_TV,
    ls="--",
    lw=1.5,
    zorder=2,
    label=f"do(TV=60) = {dbl_mean:.1f}",
)

mid_x = weeks[-1] + 1.5
ax.annotate(
    "",
    xy=(mid_x, dbl_mean),
    xytext=(mid_x, bl_mean),
    arrowprops=dict(arrowstyle="<->", color=COLOR_TV, lw=2),
)
lift = contrast.mean("sales")
ax.text(
    mid_x + 0.5,
    (bl_mean + dbl_mean) / 2,
    f"Δ = {lift:+.1f}",
    color=COLOR_TV,
    fontweight="bold",
    va="center",
    fontsize=10,
)

ax.set_xlabel("Week")
ax.set_ylabel("Sales (region average)")
ax.set_xlim(0, weeks[-1] + 4)
ax.legend(loc="lower right", fontsize=8)
plt.tight_layout()
plt.show()

Figure 4: Counterfactual analysis in data space. Observed sales (points) and model fit (solid line) alongside scenario predictions from do(). The vertical gap between bands is the causal effect of doubling TV spend.

The horizontal bands show the do() operator’s predictions: what the model expects average sales to be under each fixed-spend scenario.

fig, axes = plt.subplots(1, 2, figsize=(FIG_WIDTH, FIG_HEIGHT))

ax = axes[0]
for draws, color, label in [
    (r_baseline.draws("sales"), COLOR_BASELINE, "TV = 30 (baseline)"),
    (r_double_tv.draws("sales"), COLOR_TV, "TV = 60 (doubled)"),
]:
    x_kde, y_kde, _ = az.kde(draws)
    ax.plot(x_kde, y_kde, color=color, lw=2, label=label)
    ax.fill_between(x_kde, y_kde, alpha=0.3, color=color)
    ax.axvline(draws.mean(), color=color, ls="--", alpha=0.7, lw=1)
ax.set_xlabel("Mean sales")
ax.set_ylabel("Density")
ax.legend(fontsize=8)
ax.set_title("Scenario comparison")

ax = axes[1]
incr = contrast.draws("sales")
x_kde, y_kde, _ = az.kde(incr)
ax.plot(x_kde, y_kde, color=COLOR_TV, lw=2)
ax.fill_between(x_kde, y_kde, alpha=0.3, color=COLOR_TV)
hdi = contrast.hdi("sales", prob=0.94)
ax.axvspan(hdi[0], hdi[1], alpha=0.15, color=COLOR_TV, label="94% HDI")
ax.axvline(incr.mean(), color=COLOR_TV, ls="--", lw=1.5)
ax.axvline(0, color="black", ls=":", alpha=0.4)
ax.set_xlabel("Incremental sales (TV 60 − TV 30)")
ax.set_ylabel("Density")
ax.legend(fontsize=8)
ax.set_title("Causal effect of doubling TV")

plt.tight_layout()
plt.show()

Figure 5: Left: posterior distributions of mean sales under baseline and doubled TV spend. Right: posterior distribution of the causal effect (incremental sales from doubling TV) with 94% HDI.

With saturation transforms, doubling TV spend does not double the effect — diminishing returns kick in.

Channel contribution

How much does each channel contribute to sales at the same spend level? We set each channel to 30 while holding the other at zero, then compare the lift over a no-spend baseline.

r_no_spend = model.do(
    set={"tv": 0.0, "digital": 0.0}, simulate_over="time", kind="mean"
)
r_tv_only = model.do(
    set={"tv": 30.0, "digital": 0.0}, simulate_over="time", kind="mean"
)
r_dig_only = model.do(
    set={"tv": 0.0, "digital": 30.0}, simulate_over="time", kind="mean"
)

tv_lift = r_tv_only - r_no_spend
dig_lift = r_dig_only - r_no_spend

/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1376975462.py:1: UserWarning: Intervention value 0.00 for 'tv' is outside the observed data range [10.71, 48.65]. Results are extrapolations and should be interpreted with caution.
  r_no_spend = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1376975462.py:1: UserWarning: Intervention value 0.00 for 'digital' is outside the observed data range [5.18, 29.81]. Results are extrapolations and should be interpreted with caution.
  r_no_spend = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1376975462.py:4: UserWarning: Intervention value 0.00 for 'digital' is outside the observed data range [5.18, 29.81]. Results are extrapolations and should be interpreted with caution.
  r_tv_only = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1376975462.py:7: UserWarning: Intervention value 0.00 for 'tv' is outside the observed data range [10.71, 48.65]. Results are extrapolations and should be interpreted with caution.
  r_dig_only = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1376975462.py:7: UserWarning: Intervention value 30.00 for 'digital' is outside the observed data range [5.18, 29.81]. Results are extrapolations and should be interpreted with caution.
  r_dig_only = model.do(

fig, ax = plt.subplots(figsize=(FIG_WIDTH, FIG_HEIGHT))

for draws_obj, color, label in [
    (tv_lift, COLOR_TV, "TV (spend = 30)"),
    (dig_lift, COLOR_DIGITAL, "Digital (spend = 30)"),
]:
    draws = draws_obj.draws("sales")
    x_kde, y_kde, _ = az.kde(draws)
    ax.plot(x_kde, y_kde, color=color, lw=2, label=label)
    ax.fill_between(x_kde, y_kde, alpha=0.25, color=color)
    ax.axvline(draws.mean(), color=color, ls="--", alpha=0.7, lw=1)

ax.axvline(0, color="black", ls=":", alpha=0.4)
ax.set_xlabel("Incremental sales lift vs no spend")
ax.set_ylabel("Density")
ax.legend()
plt.tight_layout()
plt.show()

Figure 6: Posterior distributions of incremental sales lift from TV and digital channels at equal spend (30 units each), relative to a no-spend baseline.

Diminishing returns curve

The saturation transform means each additional unit of spend produces less incremental sales. We trace out this curve by running do() at multiple spend levels.

spend_levels = np.array([0, 5, 10, 20, 30, 40, 50, 60, 80])

tv_draws_list, dig_draws_list = [], []

for s in spend_levels:
    r_tv = model.do(
        set={"tv": float(s), "digital": 0.0}, simulate_over="time", kind="mean"
    )
    lift_tv = r_tv - r_no_spend
    tv_draws_list.append(lift_tv.draws("sales"))

    r_dig = model.do(
        set={"tv": 0.0, "digital": float(s)}, simulate_over="time", kind="mean"
    )
    lift_dig = r_dig - r_no_spend
    dig_draws_list.append(lift_dig.draws("sales"))

tv_draws_2d = np.column_stack(tv_draws_list)
dig_draws_2d = np.column_stack(dig_draws_list)

fig, ax = plt.subplots(figsize=(FIG_WIDTH, FIG_HEIGHT))
ax.plot(
    spend_levels, tv_draws_2d.mean(axis=0), "o-", color=COLOR_TV, label="TV", lw=2, ms=5
)
hdi_1 = az.hdi(tv_draws_2d, prob=0.94, axis=0)
ax.fill_between(spend_levels, hdi_1[:, 0], hdi_1[:, 1], alpha=0.15, color=COLOR_TV)
ax.plot(
    spend_levels,
    dig_draws_2d.mean(axis=0),
    "s-",
    color=COLOR_DIGITAL,
    label="Digital",
    lw=2,
    ms=5,
)
hdi_2 = az.hdi(dig_draws_2d, prob=0.94, axis=0)
ax.fill_between(spend_levels, hdi_2[:, 0], hdi_2[:, 1], alpha=0.15, color=COLOR_DIGITAL)
ax.set_xlabel("Channel spend (units)")
ax.set_ylabel("Incremental sales lift vs no spend")
ax.legend()
ax.axhline(0, color="black", ls=":", alpha=0.3)
plt.tight_layout()
plt.show()

/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:6: UserWarning: Intervention value 0.00 for 'tv' is outside the observed data range [10.71, 48.65]. Results are extrapolations and should be interpreted with caution.
  r_tv = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:6: UserWarning: Intervention value 0.00 for 'digital' is outside the observed data range [5.18, 29.81]. Results are extrapolations and should be interpreted with caution.
  r_tv = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:12: UserWarning: Intervention value 0.00 for 'tv' is outside the observed data range [10.71, 48.65]. Results are extrapolations and should be interpreted with caution.
  r_dig = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:12: UserWarning: Intervention value 0.00 for 'digital' is outside the observed data range [5.18, 29.81]. Results are extrapolations and should be interpreted with caution.
  r_dig = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:6: UserWarning: Intervention value 5.00 for 'tv' is outside the observed data range [10.71, 48.65]. Results are extrapolations and should be interpreted with caution.
  r_tv = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:6: UserWarning: Intervention value 0.00 for 'digital' is outside the observed data range [5.18, 29.81]. Results are extrapolations and should be interpreted with caution.
  r_tv = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:12: UserWarning: Intervention value 0.00 for 'tv' is outside the observed data range [10.71, 48.65]. Results are extrapolations and should be interpreted with caution.
  r_dig = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:12: UserWarning: Intervention value 5.00 for 'digital' is outside the observed data range [5.18, 29.81]. Results are extrapolations and should be interpreted with caution.
  r_dig = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:6: UserWarning: Intervention value 10.00 for 'tv' is outside the observed data range [10.71, 48.65]. Results are extrapolations and should be interpreted with caution.
  r_tv = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:6: UserWarning: Intervention value 0.00 for 'digital' is outside the observed data range [5.18, 29.81]. Results are extrapolations and should be interpreted with caution.
  r_tv = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:12: UserWarning: Intervention value 0.00 for 'tv' is outside the observed data range [10.71, 48.65]. Results are extrapolations and should be interpreted with caution.
  r_dig = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:6: UserWarning: Intervention value 0.00 for 'digital' is outside the observed data range [5.18, 29.81]. Results are extrapolations and should be interpreted with caution.
  r_tv = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:12: UserWarning: Intervention value 0.00 for 'tv' is outside the observed data range [10.71, 48.65]. Results are extrapolations and should be interpreted with caution.
  r_dig = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:6: UserWarning: Intervention value 0.00 for 'digital' is outside the observed data range [5.18, 29.81]. Results are extrapolations and should be interpreted with caution.
  r_tv = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:12: UserWarning: Intervention value 0.00 for 'tv' is outside the observed data range [10.71, 48.65]. Results are extrapolations and should be interpreted with caution.
  r_dig = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:12: UserWarning: Intervention value 30.00 for 'digital' is outside the observed data range [5.18, 29.81]. Results are extrapolations and should be interpreted with caution.
  r_dig = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:6: UserWarning: Intervention value 0.00 for 'digital' is outside the observed data range [5.18, 29.81]. Results are extrapolations and should be interpreted with caution.
  r_tv = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:12: UserWarning: Intervention value 0.00 for 'tv' is outside the observed data range [10.71, 48.65]. Results are extrapolations and should be interpreted with caution.
  r_dig = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:12: UserWarning: Intervention value 40.00 for 'digital' is outside the observed data range [5.18, 29.81]. Results are extrapolations and should be interpreted with caution.
  r_dig = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:6: UserWarning: Intervention value 50.00 for 'tv' is outside the observed data range [10.71, 48.65]. Results are extrapolations and should be interpreted with caution.
  r_tv = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:6: UserWarning: Intervention value 0.00 for 'digital' is outside the observed data range [5.18, 29.81]. Results are extrapolations and should be interpreted with caution.
  r_tv = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:12: UserWarning: Intervention value 0.00 for 'tv' is outside the observed data range [10.71, 48.65]. Results are extrapolations and should be interpreted with caution.
  r_dig = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:12: UserWarning: Intervention value 50.00 for 'digital' is outside the observed data range [5.18, 29.81]. Results are extrapolations and should be interpreted with caution.
  r_dig = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:6: UserWarning: Intervention value 60.00 for 'tv' is outside the observed data range [10.71, 48.65]. Results are extrapolations and should be interpreted with caution.
  r_tv = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:6: UserWarning: Intervention value 0.00 for 'digital' is outside the observed data range [5.18, 29.81]. Results are extrapolations and should be interpreted with caution.
  r_tv = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:12: UserWarning: Intervention value 0.00 for 'tv' is outside the observed data range [10.71, 48.65]. Results are extrapolations and should be interpreted with caution.
  r_dig = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:12: UserWarning: Intervention value 60.00 for 'digital' is outside the observed data range [5.18, 29.81]. Results are extrapolations and should be interpreted with caution.
  r_dig = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:6: UserWarning: Intervention value 80.00 for 'tv' is outside the observed data range [10.71, 48.65]. Results are extrapolations and should be interpreted with caution.
  r_tv = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:6: UserWarning: Intervention value 0.00 for 'digital' is outside the observed data range [5.18, 29.81]. Results are extrapolations and should be interpreted with caution.
  r_tv = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:12: UserWarning: Intervention value 0.00 for 'tv' is outside the observed data range [10.71, 48.65]. Results are extrapolations and should be interpreted with caution.
  r_dig = model.do(
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1829389347.py:12: UserWarning: Intervention value 80.00 for 'digital' is outside the observed data range [5.18, 29.81]. Results are extrapolations and should be interpreted with caution.
  r_dig = model.do(

Figure 7: Posterior mean incremental sales as a function of spend for each channel. The concave shape reflects logistic saturation.

---

Model 1 captures the nonlinear response of each channel — adstock and saturation transform raw spend into effective exposure. But it treats each channel as an independent input. In reality, upper-funnel spend (brand campaigns, TV) doesn’t just drive sales directly — it also lifts search traffic, website visits, and other lower-funnel indicators that drive purchases in turn. Ignoring this funnel structure means upper-funnel spend looks weak while the downstream indicators it generated get all the credit.

Model 2: Marketing funnel with mediation

Upper-funnel advertising (brand campaigns, TV, YouTube) doesn’t just drive sales directly. It also lifts awareness, which surfaces as increased organic search traffic, website visits, and lower-funnel engagement — and those indicators drive purchases in turn.

This causal chain is a mediation structure: brand spend affects sales both directly and indirectly through a lower-funnel mediator. Standard MMMs that ignore this funnel face an attribution problem:

• Brand spend looks weak — part of its effect is “stolen” by search/organic variables
• Search traffic looks strong — it’s capturing brand’s indirect contribution

By modelling the funnel explicitly as a path model, we can decompose brand’s total effect into its direct and indirect components, and make better budget allocation decisions.

This model uses cross-sectional data to focus on the mediation structure itself — the funnel logic applies equally to panel and cross-sectional settings.

The causal structure

Figure 8: Marketing funnel DAG. Brand spend drives sales directly (c) and indirectly through search traffic (a × b). Performance spend has a direct effect only. Seasonality is a common cause of both search traffic and sales.

Effect	Formula	Interpretation
Direct	c	Brand → Sales, not through search
Indirect	a × b	Brand → Search → Sales
Total	c + a × b	Full causal effect of brand spend

Seasonality is a confounder of search traffic and sales — it must appear in both equations to avoid bias, but it is not on the causal path from brand spend to sales.

Simulate data

rng = np.random.default_rng(42)
n = 500

true_a = 0.6  # brand → search
true_b = 0.5  # search → sales
true_c = 0.25  # brand → sales (direct)
true_d = 0.7  # perf → sales
true_season_search = 0.4
true_season_sales = 0.3

season = rng.normal(size=n)
brand_spend = rng.uniform(0, 10, size=n)
perf_spend = rng.uniform(0, 8, size=n)

search_traffic = (
    true_a * brand_spend + true_season_search * season + rng.normal(scale=0.8, size=n)
)

sales = (
    true_c * brand_spend
    + true_b * search_traffic
    + true_d * perf_spend
    + true_season_sales * season
    + rng.normal(scale=1.0, size=n)
)

df = pd.DataFrame({
    "brand_spend": brand_spend,
    "perf_spend": perf_spend,
    "search_traffic": search_traffic,
    "season": season,
    "sales": sales,
})

true_indirect = true_a * true_b
true_total_brand = true_c + true_indirect

print(f"True direct brand effect:   {true_c}")
print(f"True indirect brand effect: {true_indirect}")
print(f"True total brand effect:    {true_total_brand}")
print(f"Fraction through funnel:    {true_indirect / true_total_brand:.0%}")
df.head()

True direct brand effect:   0.25
True indirect brand effect: 0.3
True total brand effect:    0.55
Fraction through funnel:    55%

	brand_spend	perf_spend	search_traffic	season	sales
0	3.722616	5.996949	2.806107	0.304717	6.062537
1	1.536129	7.934462	0.597351	-1.039984	5.955883
2	6.008404	4.251310	4.169178	0.750451	6.535680
3	1.196726	5.279966	1.716349	0.940565	4.349525
4	3.649194	2.419842	0.761513	-1.951035	2.111758

Over half of brand’s total causal effect on sales flows through the search traffic mediator — invisible to a model that doesn’t represent the funnel.

Visualise the data

fig, axes = plt.subplots(1, 3, figsize=(FIG_WIDTH, FIG_HEIGHT))

axes[0].scatter(brand_spend, search_traffic, alpha=0.3, s=10, color=COLOR_BRAND)
axes[0].set_xlabel("Brand spend")
axes[0].set_ylabel("Search traffic")

axes[1].scatter(search_traffic, sales, alpha=0.3, s=10, color=COLOR_SEARCH)
axes[1].set_xlabel("Search traffic")
axes[1].set_ylabel("Sales")

axes[2].scatter(brand_spend, sales, alpha=0.3, s=10, color=COLOR_BRAND)
axes[2].set_xlabel("Brand spend")
axes[2].set_ylabel("Sales")

plt.tight_layout()
plt.show()

Figure 9: Pairwise scatter plots of brand spend, search traffic, and sales. The positive association between brand spend and search traffic reflects the upper-to-lower funnel link.

The attribution trap

WarningWhy naive regression misleads

Regressing sales ~ brand_spend + perf_spend + search_traffic + season estimates the direct effect of brand spend (≈ 0.25), not the total effect (≈ 0.55). By conditioning on the mediator search_traffic, the regression blocks the indirect path brand → search → sales.

This is correct if you want the direct effect. But for budget allocation — “how much total revenue does brand spend generate?” — you need the total effect, which requires modelling the funnel as a causal structure.

Specify and fit the path model

The spec encodes both equations in the funnel, with labeled coefficients for the paths we care about and a defined parameter for the indirect effect.

spec = """
search_traffic ~ a*brand_spend + season
sales ~ b*search_traffic + c*brand_spend + d*perf_spend + season
indirect := a*b
"""

model = pathmc.model(spec, data=df)

Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.

model.graph()

model.equations()

\begin{aligned} \beta_{search,traffic} &\sim \text{Normal}(mu=0,\, sigma=10) \\ \sigma_{search,traffic} &\sim \text{HalfNormal}(sigma=1) \\ \beta_{sales} &\sim \text{Normal}(mu=0,\, sigma=10) \\ \sigma_{sales} &\sim \text{HalfNormal}(sigma=1) \\[6pt] \mu_{search,traffic} &= \beta_{0,\,search,traffic} + a \cdot \mathrm{brand\_spend} + \mathrm{season} \\ \mathrm{search\_traffic} &\sim \text{Normal}(\mu_{search,traffic},\, \sigma_{search,traffic}) \\ \mu_{sales} &= \beta_{0,\,sales} \\ &\quad + b \cdot \mathrm{search\_traffic} \\ &\quad + c \cdot \mathrm{brand\_spend} \\ &\quad + d \cdot \mathrm{perf\_spend} \\ &\quad + \mathrm{season} \\ \mathrm{sales} &\sim \text{Normal}(\mu_{sales},\, \sigma_{sales}) \\ indirect &\equiv a \cdot b \end{aligned}

pm.model_to_graphviz(model.pymc_model)

PyMC plate diagram for the marketing funnel path model.

Sample

idata = model.fit(random_seed=42)

Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
NUTS[nutpie]: [beta_sales, sigma_search_traffic, beta_search_traffic, sigma_sales]
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.

Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.

Results

model.summary()

	mean	sd	eti89_lb	eti89_ub	ess_bulk	ess_tail	r_hat	mcse_mean	mcse_sd
beta_sales[Intercept]	0.082508	0.119687	-0.111247	0.272075	2353.949909	2589.308904	1.000677	0.002466	0.001742
beta_sales[search_traffic]	0.480797	0.056485	0.392478	0.572599	1214.665163	1425.753515	1.000525	0.001623	0.001184
beta_sales[brand_spend]	0.264196	0.037691	0.202557	0.323914	1176.472505	1546.068704	1.000928	0.001103	0.000806
beta_sales[perf_spend]	0.673370	0.019663	0.641903	0.705246	2652.147477	2713.501299	1.000843	0.000382	0.000275
beta_sales[season]	0.265064	0.050934	0.184085	0.348701	3282.209749	3022.040658	0.999869	0.000891	0.000624
...	...	...	...	...	...	...	...	...	...
mu_search_traffic[495]	4.298691	0.055507	4.209698	4.385626	4973.440710	3561.104854	1.000298	0.000787	0.000551
mu_search_traffic[496]	4.362853	0.087283	4.221610	4.506116	4597.811582	3109.240903	1.000339	0.001289	0.000887
mu_search_traffic[497]	-0.098858	0.067698	-0.205729	0.006947	2657.122110	2649.520260	1.000851	0.001313	0.000957
mu_search_traffic[498]	1.852374	0.082187	1.722139	1.985987	5396.942356	3271.651578	0.999540	0.001118	0.000816
mu_search_traffic[499]	4.769232	0.081254	4.639807	4.900159	4212.959167	3369.305338	1.000458	0.001251	0.000855

1010 rows × 9 columns

model.effects_summary()

	mean	sd	hdi_3%	hdi_97%
name
a	0.619181	0.012270	0.597171	0.643336
b	0.480797	0.056485	0.382269	0.595773
c	0.264196	0.037691	0.191312	0.334436
d	0.673370	0.019663	0.636183	0.709723
indirect	0.297681	0.035282	0.234848	0.366869

TipReading the effects table

• a: brand → search (true = 0.6)
• b: search → sales (true = 0.5)
• c: brand → sales, direct (true = 0.25)
• d: perf → sales (true = 0.7)
• indirect: a × b (true = 0.3)

The total brand effect is c + indirect ≈ 0.55.

Path-specific effects

The effect() method multiplies posterior draws along a specified path, giving the full posterior distribution of path-specific effects.

indirect_effect = model.effect("brand_spend -> search_traffic -> sales")
print(f"Indirect effect (brand → search → sales): {indirect_effect}")

direct_effect = model.effect("brand_spend -> sales")
print(f"Direct effect (brand → sales):            {direct_effect}")

Indirect effect (brand → search → sales): EffectResult('brand_spend -> search_traffic -> sales', mean=0.2977, 94% HDI=[0.2348, 0.3669])
Direct effect (brand → sales):            EffectResult('brand_spend -> sales', mean=0.2642, 94% HDI=[0.1913, 0.3344])

Decomposing brand’s total effect

fig, ax = plt.subplots(figsize=(FIG_WIDTH, FIG_HEIGHT))

a_draws = (
    idata
    .posterior["beta_search_traffic"]
    .sel(search_traffic_predictors="brand_spend")
    .values.flatten()
)
b_draws = (
    idata
    .posterior["beta_sales"]
    .sel(sales_predictors="search_traffic")
    .values.flatten()
)
c_draws = (
    idata.posterior["beta_sales"].sel(sales_predictors="brand_spend").values.flatten()
)
indirect_draws = a_draws * b_draws
total_draws = c_draws + indirect_draws

for draws, color, label, true_val in [
    (c_draws, COLOR_BRAND, "Direct (c)", true_c),
    (indirect_draws, COLOR_INDIRECT, "Indirect (a × b)", true_indirect),
    (total_draws, COLOR_SEARCH, "Total (c + a × b)", true_total_brand),
]:
    x_kde, y_kde, _ = az.kde(draws)
    ax.plot(x_kde, y_kde, color=color, lw=2, label=label)
    ax.fill_between(x_kde, y_kde, alpha=0.25, color=color)
    ax.axvline(true_val, color=color, ls="--", lw=1.5, alpha=0.7)

ax.axvline(0, color="black", ls=":", alpha=0.3)
ax.set_xlabel("Effect size (per unit brand spend)")
ax.set_ylabel("Density")
ax.legend()
plt.tight_layout()
plt.show()

Figure 10: Posterior distributions of brand spend’s direct effect (c), indirect effect (a × b), and total effect (c + a × b). True values shown as dashed lines.

The indirect effect is comparable in magnitude to the direct effect — a naive regression that conditions on search traffic would miss roughly half of brand’s total contribution.

Identification and standardized effects

print(
    f"Is brand_spend → sales identifiable? {model.is_identifiable('brand_spend', 'sales')}"
)
print(
    f"Valid adjustment sets:                {model.adjustment_sets('brand_spend', 'sales')}"
)

Is brand_spend → sales identifiable? True
Valid adjustment sets:                [set()]

The empty set is a valid adjustment set because brand spend has no confounders (it is exogenous in this DAG).

model.standardized()

	predictor	outcome	mean	sd	hdi_3%	hdi_97%
name
a	brand_spend	search_traffic	0.904982	0.017933	0.872812	0.940286
b	search_traffic	sales	0.378172	0.044428	0.300674	0.468606
c	brand_spend	sales	0.303721	0.043329	0.219934	0.384470
d	perf_spend	sales	0.611731	0.017863	0.577948	0.644756

Causal queries

Total causal effect of brand spend

The .ate() method computes the average treatment effect by contrasting two interventions — including both the direct path and the indirect path through search traffic:

brand_ate = model.ate("sales", "brand_spend", values=(2.0, 8.0))
brand_ate

ATE of brand_spend on sales

Mean	3.37
94% HDI	[3.19, 3.55]
P(> 0)	1.00
Draws	4000

Methods: .hdi() .prob() .summary() .draws() .by_time() .dataset

The ATE should be close to 0.55 × 6 = 3.3 — the total effect times the intervention size.

Probability query

p = model.prob("sales > 0", set={"brand_spend": 8.0})
print(f"P(sales > 0 | do(brand_spend=8)): {p:.2f}")

Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Sampling: [sales, search_traffic]

P(sales > 0 | do(brand_spend=8)): 1.00

Comparing channel effectiveness

Which channel produces more incremental sales per unit of spend?

r_base = model.do(set={"brand_spend": 0.0, "perf_spend": 0.0})

r_brand_only = model.do(set={"brand_spend": 5.0, "perf_spend": 0.0})
r_perf_only = model.do(set={"brand_spend": 0.0, "perf_spend": 5.0})

brand_lift = r_brand_only - r_base
perf_lift = r_perf_only - r_base

/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1173457870.py:1: UserWarning: Intervention value 0.00 for 'brand_spend' is outside the observed data range [0.01, 9.99]. Results are extrapolations and should be interpreted with caution.
  r_base = model.do(set={"brand_spend": 0.0, "perf_spend": 0.0})
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1173457870.py:1: UserWarning: Intervention value 0.00 for 'perf_spend' is outside the observed data range [0.00, 8.00]. Results are extrapolations and should be interpreted with caution.
  r_base = model.do(set={"brand_spend": 0.0, "perf_spend": 0.0})
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1173457870.py:3: UserWarning: Intervention value 0.00 for 'perf_spend' is outside the observed data range [0.00, 8.00]. Results are extrapolations and should be interpreted with caution.
  r_brand_only = model.do(set={"brand_spend": 5.0, "perf_spend": 0.0})
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1173457870.py:4: UserWarning: Intervention value 0.00 for 'brand_spend' is outside the observed data range [0.01, 9.99]. Results are extrapolations and should be interpreted with caution.
  r_perf_only = model.do(set={"brand_spend": 0.0, "perf_spend": 5.0})

Incremental sales from 5 units of brand spend:

brand_lift

DoResult — 4000 draws, 5 variables

variable	mean	94% HDI
brand_spend	5.00	[5.00, 5.00]
season	0.00	[0.00, 0.00]
perf_spend	0.00	[0.00, 0.00]
search_traffic	3.10	[2.99, 3.22]
sales	2.81	[2.65, 2.96]

Methods: .draws() .mean() .hdi() .by_time() .plot() .dataset

Incremental sales from 5 units of performance spend:

perf_lift

DoResult — 4000 draws, 5 variables

variable	mean	94% HDI
brand_spend	0.00	[0.00, 0.00]
season	0.00	[0.00, 0.00]
perf_spend	5.00	[5.00, 5.00]
search_traffic	0.00	[0.00, 0.00]
sales	3.37	[3.18, 3.55]

Methods: .draws() .mean() .hdi() .by_time() .plot() .dataset

fig, ax = plt.subplots(figsize=(FIG_WIDTH, FIG_HEIGHT))

for lift_obj, color, label in [
    (brand_lift, COLOR_BRAND, "Brand (total, incl. funnel)"),
    (perf_lift, COLOR_PERF, "Performance (direct only)"),
]:
    draws = lift_obj.draws("sales")
    x_kde, y_kde, _ = az.kde(draws)
    ax.plot(x_kde, y_kde, color=color, lw=2, label=label)
    ax.fill_between(x_kde, y_kde, alpha=0.25, color=color)
    ax.axvline(draws.mean(), color=color, ls="--", alpha=0.7, lw=1)

ax.axvline(0, color="black", ls=":", alpha=0.3)
ax.set_xlabel("Incremental sales vs zero spend")
ax.set_ylabel("Density")
ax.legend()
plt.tight_layout()
plt.show()

Figure 11: Posterior distributions of incremental sales from 5 units of brand spend vs 5 units of performance spend. Brand’s total effect includes the indirect path through search traffic.

What does brand spend do to search traffic?

Because the model includes the funnel equation, do() also gives us the downstream effect on the mediator.

r_low = model.do(set={"brand_spend": 2.0})
r_high = model.do(set={"brand_spend": 8.0})

r_low

DoResult — 4000 draws, 5 variables

variable	mean	94% HDI
brand_spend	2.00	[2.00, 2.00]
season	-0.01	[-0.01, -0.01]
perf_spend	4.04	[4.04, 4.04]
search_traffic	1.00	[0.91, 1.10]
sales	3.81	[3.67, 3.94]

Methods: .draws() .mean() .hdi() .by_time() .plot() .dataset

r_high

DoResult — 4000 draws, 5 variables

variable	mean	94% HDI
brand_spend	8.00	[8.00, 8.00]
season	-0.01	[-0.01, -0.01]
perf_spend	4.04	[4.04, 4.04]
search_traffic	4.72	[4.62, 4.81]
sales	7.18	[7.05, 7.31]

Methods: .draws() .mean() .hdi() .by_time() .plot() .dataset

The expected lift in search traffic from 6 units of brand spend is true_a × 6.

Budget allocation

spend_levels = np.array([0, 1, 2, 3, 5, 7, 10])

brand_draws_list, perf_draws_list = [], []

for s in spend_levels:
    r_b = model.do(set={"brand_spend": float(s), "perf_spend": 0.0})
    lift_b = r_b - r_base
    brand_draws_list.append(lift_b.draws("sales"))

    r_p = model.do(set={"brand_spend": 0.0, "perf_spend": float(s)})
    lift_p = r_p - r_base
    perf_draws_list.append(lift_p.draws("sales"))

brand_draws_2d = np.column_stack(brand_draws_list)
perf_draws_2d = np.column_stack(perf_draws_list)

fig, ax = plt.subplots(figsize=(FIG_WIDTH, FIG_HEIGHT))
ax.plot(
    spend_levels,
    brand_draws_2d.mean(axis=0),
    "o-",
    color=COLOR_BRAND,
    label="Brand (total)",
    lw=2,
    ms=5,
)
hdi_3 = az.hdi(brand_draws_2d, prob=0.94, axis=0)
ax.fill_between(spend_levels, hdi_3[:, 0], hdi_3[:, 1], alpha=0.15, color=COLOR_BRAND)
ax.plot(
    spend_levels,
    perf_draws_2d.mean(axis=0),
    "s-",
    color=COLOR_PERF,
    label="Performance",
    lw=2,
    ms=5,
)
hdi_4 = az.hdi(perf_draws_2d, prob=0.94, axis=0)
ax.fill_between(spend_levels, hdi_4[:, 0], hdi_4[:, 1], alpha=0.15, color=COLOR_PERF)

ax.plot(
    spend_levels,
    [true_total_brand * s for s in spend_levels],
    "--",
    color=COLOR_BRAND,
    alpha=0.5,
    label=f"True brand slope ({true_total_brand})",
)
ax.plot(
    spend_levels,
    [true_d * s for s in spend_levels],
    "--",
    color=COLOR_PERF,
    alpha=0.5,
    label=f"True perf slope ({true_d})",
)

ax.set_xlabel("Channel spend (units)")
ax.set_ylabel("Incremental sales vs zero spend")
ax.legend(fontsize=8)
ax.axhline(0, color="black", ls=":", alpha=0.3)
plt.tight_layout()
plt.show()

/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1959139574.py:6: UserWarning: Intervention value 0.00 for 'brand_spend' is outside the observed data range [0.01, 9.99]. Results are extrapolations and should be interpreted with caution.
  r_b = model.do(set={"brand_spend": float(s), "perf_spend": 0.0})
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1959139574.py:6: UserWarning: Intervention value 0.00 for 'perf_spend' is outside the observed data range [0.00, 8.00]. Results are extrapolations and should be interpreted with caution.
  r_b = model.do(set={"brand_spend": float(s), "perf_spend": 0.0})
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1959139574.py:10: UserWarning: Intervention value 0.00 for 'brand_spend' is outside the observed data range [0.01, 9.99]. Results are extrapolations and should be interpreted with caution.
  r_p = model.do(set={"brand_spend": 0.0, "perf_spend": float(s)})
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1959139574.py:10: UserWarning: Intervention value 0.00 for 'perf_spend' is outside the observed data range [0.00, 8.00]. Results are extrapolations and should be interpreted with caution.
  r_p = model.do(set={"brand_spend": 0.0, "perf_spend": float(s)})
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1959139574.py:6: UserWarning: Intervention value 0.00 for 'perf_spend' is outside the observed data range [0.00, 8.00]. Results are extrapolations and should be interpreted with caution.
  r_b = model.do(set={"brand_spend": float(s), "perf_spend": 0.0})
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1959139574.py:10: UserWarning: Intervention value 0.00 for 'brand_spend' is outside the observed data range [0.01, 9.99]. Results are extrapolations and should be interpreted with caution.
  r_p = model.do(set={"brand_spend": 0.0, "perf_spend": float(s)})
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1959139574.py:6: UserWarning: Intervention value 0.00 for 'perf_spend' is outside the observed data range [0.00, 8.00]. Results are extrapolations and should be interpreted with caution.
  r_b = model.do(set={"brand_spend": float(s), "perf_spend": 0.0})
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1959139574.py:10: UserWarning: Intervention value 0.00 for 'brand_spend' is outside the observed data range [0.01, 9.99]. Results are extrapolations and should be interpreted with caution.
  r_p = model.do(set={"brand_spend": 0.0, "perf_spend": float(s)})
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1959139574.py:6: UserWarning: Intervention value 0.00 for 'perf_spend' is outside the observed data range [0.00, 8.00]. Results are extrapolations and should be interpreted with caution.
  r_b = model.do(set={"brand_spend": float(s), "perf_spend": 0.0})
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1959139574.py:10: UserWarning: Intervention value 0.00 for 'brand_spend' is outside the observed data range [0.01, 9.99]. Results are extrapolations and should be interpreted with caution.
  r_p = model.do(set={"brand_spend": 0.0, "perf_spend": float(s)})
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1959139574.py:6: UserWarning: Intervention value 0.00 for 'perf_spend' is outside the observed data range [0.00, 8.00]. Results are extrapolations and should be interpreted with caution.
  r_b = model.do(set={"brand_spend": float(s), "perf_spend": 0.0})
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1959139574.py:10: UserWarning: Intervention value 0.00 for 'brand_spend' is outside the observed data range [0.01, 9.99]. Results are extrapolations and should be interpreted with caution.
  r_p = model.do(set={"brand_spend": 0.0, "perf_spend": float(s)})
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1959139574.py:6: UserWarning: Intervention value 0.00 for 'perf_spend' is outside the observed data range [0.00, 8.00]. Results are extrapolations and should be interpreted with caution.
  r_b = model.do(set={"brand_spend": float(s), "perf_spend": 0.0})
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1959139574.py:10: UserWarning: Intervention value 0.00 for 'brand_spend' is outside the observed data range [0.01, 9.99]. Results are extrapolations and should be interpreted with caution.
  r_p = model.do(set={"brand_spend": 0.0, "perf_spend": float(s)})
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1959139574.py:6: UserWarning: Intervention value 10.00 for 'brand_spend' is outside the observed data range [0.01, 9.99]. Results are extrapolations and should be interpreted with caution.
  r_b = model.do(set={"brand_spend": float(s), "perf_spend": 0.0})
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1959139574.py:6: UserWarning: Intervention value 0.00 for 'perf_spend' is outside the observed data range [0.00, 8.00]. Results are extrapolations and should be interpreted with caution.
  r_b = model.do(set={"brand_spend": float(s), "perf_spend": 0.0})
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1959139574.py:10: UserWarning: Intervention value 0.00 for 'brand_spend' is outside the observed data range [0.01, 9.99]. Results are extrapolations and should be interpreted with caution.
  r_p = model.do(set={"brand_spend": 0.0, "perf_spend": float(s)})
/var/folders/pd/p2qnky2x3xl4w3mgc4lct2200000gn/T/ipykernel_23316/1959139574.py:10: UserWarning: Intervention value 10.00 for 'perf_spend' is outside the observed data range [0.00, 8.00]. Results are extrapolations and should be interpreted with caution.
  r_p = model.do(set={"brand_spend": 0.0, "perf_spend": float(s)})

Figure 12: Sales response curves for brand and performance channels via do() at multiple spend levels. Brand’s curve includes the indirect effect through search. Dashed lines show true slopes.

In this linear model the response curves are straight lines. Performance spend has a steeper marginal effect (0.7 vs 0.55 per unit), but brand spend is not far behind — and a flat regression that conditions on search traffic would underestimate brand at 0.25.

---

Model 2 separates direct from indirect effects, solving the attribution problem. But it treats all observations identically — a single coefficient for TV effectiveness everywhere. In practice, TV might drive more sales in one region due to higher viewership, while digital outperforms in another where the audience skews younger. A hierarchical model lets each region have its own media coefficients while borrowing strength from the full panel.

Model 3: Hierarchical geo-varying effects

Marketing channels rarely have the same effectiveness everywhere. This model fits a hierarchical panel MMM with geo-varying media coefficients using random slopes. Each region gets its own TV and digital effect, partially pooled toward a shared mean — borrowing strength across geos while allowing genuine heterogeneity.

No transforms (adstock, saturation) are used here; the focus is on the hierarchical structure itself. In practice, you would combine random slopes with the transforms from Model 1.

The model

For region g at week t:

\text{sales}_{g,t} = \alpha_g + (\beta_{\text{tv}} + u_{g,\text{tv}}) \cdot \text{tv}_{g,t} + (\beta_{\text{dig}} + u_{g,\text{dig}}) \cdot \text{digital}_{g,t} + \beta_{\text{trend}} \cdot t + \varepsilon_{g,t}

where \alpha_g is a random intercept and u_{g,\text{tv}}, u_{g,\text{dig}} are random slopes — each drawn from a shared normal with estimated mean and scale.

Simulate panel data

rng = np.random.default_rng(42)

regions = ["North", "South", "East", "West", "Central", "Coast"]
n_weeks = 40

true_mu_tv = 0.6
true_sigma_tv = 0.15
true_mu_dig = 0.4
true_sigma_dig = 0.10

true_tv_effects = {r: rng.normal(true_mu_tv, true_sigma_tv) for r in regions}
true_dig_effects = {r: rng.normal(true_mu_dig, true_sigma_dig) for r in regions}
true_intercepts = {r: rng.uniform(40, 70) for r in regions}

rows = []
for region in regions:
    for week in range(1, n_weeks + 1):
        tv = rng.uniform(5, 40)
        digital = rng.uniform(3, 25)
        sales = (
            true_intercepts[region]
            + true_tv_effects[region] * tv
            + true_dig_effects[region] * digital
            + 0.1 * week
            + rng.normal(scale=2.0)
        )
        rows.append({
            "region": region,
            "week": week,
            "tv": tv,
            "digital": digital,
            "trend": week,
            "sales": sales,
        })

df = pd.DataFrame(rows)
print(f"Panel: {len(regions)} regions × {n_weeks} weeks = {len(df)} rows")
df.head()

Panel: 6 regions × 40 weeks = 240 rows

	region	week	tv	digital	trend	sales
0	North	1	33.967091	16.896617	1	87.953690
1	North	2	17.408409	24.355357	2	80.501138
2	North	3	32.243422	7.282052	3	84.506308
3	North	4	6.533132	6.394369	4	67.435582
4	North	5	31.066675	24.285214	5	88.876004

true_effects = pd.DataFrame({
    "region": regions,
    "tv_effect": [true_tv_effects[r] for r in regions],
    "digital_effect": [true_dig_effects[r] for r in regions],
})
true_effects

	region	tv_effect	digital_effect
0	North	0.645708	0.412784
1	South	0.444002	0.368376
2	East	0.712568	0.398320
3	West	0.741085	0.314696
4	Central	0.307345	0.487940
5	Coast	0.404673	0.477779

Visualise sales by region

fig, ax = plt.subplots(figsize=(FIG_WIDTH, FIG_HEIGHT))
palette = ["#2171b5", "#e6550d", "#31a354", "#756bb1", "#e7298a", "#636363"]
for region, color in zip(regions, palette):
    d = df[df["region"] == region]
    ax.plot(d["week"], d["sales"], label=region, color=color, alpha=0.7)
ax.set_xlabel("Week")
ax.set_ylabel("Sales")
ax.legend(ncol=3, fontsize=8)
plt.tight_layout()
plt.show()

Figure 13: Simulated weekly sales for six regions. Differences in level reflect random intercepts; differences in responsiveness to spend reflect random slopes.

Specify and fit the model

The pooling argument controls the hierarchical structure:

• "intercept": True gives each region its own baseline sales level. The compiler drops the fixed intercept either way; 0 + keeps equations() from printing a phantom 1 +
• "slopes": ["tv", "digital"] gives each region its own media coefficients, partially pooled toward a shared mean

spec = """
sales ~ 0 + tv + digital + trend
"""

model = pathmc.model(
    spec,
    data=df,
    panel={"unit": "region", "time": "week"},
    pooling={"intercept": True, "slopes": ["tv", "digital"]},
)

Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.

model.equations()

\begin{aligned} \beta_{sales} &\sim \text{Normal}(mu=0,\, sigma=10) \\ \sigma_{sales} &\sim \text{HalfNormal}(sigma=1) \\ \mu_{alpha,sales} &\sim \text{Normal}(mu=0,\, sigma=10) \\ \sigma_{alpha,sales} &\sim \text{HalfNormal}(sigma=1) \\ \alpha_{sales} &\sim \text{Normal}(mu\_alpha,\, sigma\_alpha) \\ \mu_{slope,sales,tv} &\sim \text{Normal}(mu=0,\, sigma=10) \\ \sigma_{slope,sales,tv} &\sim \text{HalfNormal}(sigma=1) \\ slope_{sales,tv} &\sim \text{Normal}(mu\_slope,\, sigma\_slope) \\ \mu_{slope,sales,digital} &\sim \text{Normal}(mu=0,\, sigma=10) \\ \sigma_{slope,sales,digital} &\sim \text{HalfNormal}(sigma=1) \\ slope_{sales,digital} &\sim \text{Normal}(mu\_slope,\, sigma\_slope) \\[6pt] \mu_{sales} &= \mathrm{tv} + \mathrm{digital} + \mathrm{trend} \\ \mathrm{sales} &\sim \text{Normal}(\mu_{sales},\, \sigma_{sales}) \end{aligned}

pm.model_to_graphviz(model.pymc_model)

PyMC plate diagram for the hierarchical panel MMM with random intercepts and random slopes for TV and digital.

Sample

This fit pins nuts_sampler="nutpie", which is the sampler that produced the published sample. Unit intercepts and slopes are sampled as a unit-normal offset, so a full model.summary() leads with the free *_offset rows and truncates the centered slopes. The nutpie run below has not converged; the centered digital-slope table in the results section is the one with the poor diagnostics.

idata = model.fit(
    draws=1000, tune=1000, chains=4, random_seed=42, nuts_sampler="nutpie"
)

Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
NUTS[nutpie]: [slope_sales_digital_offset, sigma_slope_sales_digital, mu_slope_sales_digital, slope_sales_tv_offset, sigma_slope_sales_tv, mu_slope_sales_tv, alpha_sales_offset, sigma_alpha_sales, mu_alpha_sales, beta_sales, sigma_sales]
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.

Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.

Results

Population-level coefficients

The centered mu_slope_sales_digital and slope_sales_digital rows in the table below have r-hat up to 1.10 and ESS down to 79. Those rows are not geo-level digital effects, and a mean of -0.13 is not digital effectiveness. The violins and the digital ATE below are the sums of the fixed digital coefficient and the region slope, not these rows and not the free offset rows.

dig_sum

	mean	sd	eti89_lb	eti89_ub	ess_bulk	ess_tail	r_hat	mcse_mean	mcse_sd
mu_slope_sales_digital	-0.1300	7.5	-12.0	12.0	79.0	130.0	1.1	0.85	0.6
slope_sales_digital[Central]	-0.0049	7.5	-12.0	12.0	79.0	130.0	1.1	0.85	0.6
slope_sales_digital[Coast]	-0.0650	7.5	-12.0	12.0	79.0	130.0	1.1	0.85	0.6
slope_sales_digital[East]	-0.2100	7.5	-12.0	12.0	79.0	130.0	1.1	0.85	0.6
slope_sales_digital[North]	-0.1100	7.5	-12.0	12.0	79.0	130.0	1.1	0.85	0.6
slope_sales_digital[South]	-0.1600	7.5	-12.0	12.0	79.0	130.0	1.1	0.85	0.6
slope_sales_digital[West]	-0.2600	7.5	-12.0	12.0	79.0	130.0	1.1	0.85	0.6

A full model.summary() leads with the free *_offset rows. The centered slopes above are Deterministics, so that truncation omits them.

model.summary()

	mean	sd	eti89_lb	eti89_ub	ess_bulk	ess_tail	r_hat	mcse_mean	mcse_sd
slope_sales_digital_offset[Central]	1.074407	0.625883	0.106653	2.091971	3800.410709	3184.459452	1.001309	0.010161	0.007379
slope_sales_digital_offset[Coast]	0.585805	0.590172	-0.321137	1.538282	4516.742085	3384.609401	1.000161	0.008802	0.006754
slope_sales_digital_offset[East]	-0.627445	0.592062	-1.606242	0.278678	4236.854533	3223.289025	1.002422	0.009161	0.006873
slope_sales_digital_offset[North]	0.184496	0.571846	-0.697980	1.121711	4751.143005	3281.436068	1.002060	0.008283	0.006448
slope_sales_digital_offset[South]	-0.232237	0.558328	-1.140999	0.641637	4793.347403	3518.234152	1.001893	0.008062	0.005867
...	...	...	...	...	...	...	...	...	...
alpha_sales[Coast]	41.712281	1.104655	39.931861	43.484389	4671.229048	3346.011794	1.000589	0.016186	0.011712
alpha_sales[East]	54.552780	0.869535	53.149499	55.963285	4684.781767	3535.545719	0.999776	0.012718	0.009202
alpha_sales[North]	59.412911	0.955428	57.896224	60.905045	4443.834403	3609.662074	0.999872	0.014362	0.010417
alpha_sales[South]	64.537690	1.016773	62.949577	66.166719	4586.194859	3757.129135	1.000583	0.015017	0.010751
alpha_sales[West]	48.953645	1.000628	47.331915	50.565678	4597.520207	3519.740580	1.001402	0.014783	0.010523

286 rows × 9 columns

Geo-level TV effects

The random slopes capture how each region deviates from the population mean. The posterior for each region’s TV effect is the sum of the fixed effect and the region-specific deviation.

fig, ax = plt.subplots(figsize=(FIG_WIDTH, FIG_HEIGHT * 1.2))

beta_tv = idata.posterior["beta_sales"].sel(sales_predictors="tv").values.flatten()
slope_tv = idata.posterior["slope_sales_tv"]

for i, region in enumerate(regions):
    region_slope = slope_tv.sel(unit=region).values.flatten()
    region_total = beta_tv + region_slope
    parts = ax.violinplot([region_total], positions=[i], showmedians=True, widths=0.7)
    for pc in parts["bodies"]:
        pc.set_facecolor(COLOR_TV)
        pc.set_alpha(0.4)
    for key in ["cmins", "cmaxes", "cbars", "cmedians"]:
        if key in parts:
            parts[key].set_color(COLOR_TV)

    ax.plot(i, true_tv_effects[region], "D", color="black", ms=6, zorder=5)

ax.set_xticks(range(len(regions)))
ax.set_xticklabels(regions)
ax.set_ylabel("TV effect on sales")
ax.axhline(
    true_mu_tv,
    color=COLOR_TV,
    ls="--",
    alpha=0.5,
    label=f"True pop. mean ({true_mu_tv})",
)
ax.legend(fontsize=8)
plt.tight_layout()
plt.show()

Figure 14: Posterior distributions of geo-level TV effects (violins) vs true values (black diamonds). Each region’s effect is the population mean plus a region-specific deviation, partially pooled toward the group mean.

Black diamonds show the true geo-level effects used in the DGP.

Geo-level digital effects

fig, ax = plt.subplots(figsize=(FIG_WIDTH, FIG_HEIGHT * 1.2))

beta_dig = (
    idata.posterior["beta_sales"].sel(sales_predictors="digital").values.flatten()
)
slope_dig = idata.posterior["slope_sales_digital"]

for i, region in enumerate(regions):
    region_slope = slope_dig.sel(unit=region).values.flatten()
    region_total = beta_dig + region_slope
    parts = ax.violinplot([region_total], positions=[i], showmedians=True, widths=0.7)
    for pc in parts["bodies"]:
        pc.set_facecolor(COLOR_DIGITAL)
        pc.set_alpha(0.4)
    for key in ["cmins", "cmaxes", "cbars", "cmedians"]:
        if key in parts:
            parts[key].set_color(COLOR_DIGITAL)

    ax.plot(i, true_dig_effects[region], "D", color="black", ms=6, zorder=5)

ax.set_xticks(range(len(regions)))
ax.set_xticklabels(regions)
ax.set_ylabel("Digital effect on sales")
ax.axhline(
    true_mu_dig,
    color=COLOR_DIGITAL,
    ls="--",
    alpha=0.5,
    label=f"True pop. mean ({true_mu_dig})",
)
ax.legend(fontsize=8)
plt.tight_layout()
plt.show()

Figure 15: Posterior of each region’s digital effect, the sum of the fixed digital coefficient and the region slope, against planted values (black diamonds). These sums are not the centered slope_sales_digital rows in the table above.

Shrinkage

The points below are posterior means of the fixed coefficient plus the region slope, against the planted geo effect. This panel is balanced, so the figure is not an example of shrinking regions that have less data. The horizontal lines are the fixed-coefficient means, not the group mean of those sums and not channel effectiveness.

fig, ax = plt.subplots(figsize=(FIG_WIDTH, FIG_HEIGHT))

for channel, beta_draws, slope_var, true_effects_dict, color, label in [
    ("tv", beta_tv, slope_tv, true_tv_effects, COLOR_TV, "TV"),
    ("digital", beta_dig, slope_dig, true_dig_effects, COLOR_DIGITAL, "Digital"),
]:
    true_vals = []
    post_means = []
    for region in regions:
        region_slope = slope_var.sel(unit=region).values.flatten()
        region_total = beta_draws + region_slope
        true_vals.append(true_effects_dict[region])
        post_means.append(region_total.mean())

    ax.scatter(true_vals, post_means, color=color, s=50, label=label, zorder=3)
    pop_mean = beta_draws.mean()
    ax.axhline(pop_mean, color=color, ls="--", alpha=0.3)

lims = ax.get_xlim()
ax.plot(lims, lims, "k:", alpha=0.3, label="Perfect recovery")
ax.set_xlabel("True geo-level effect")
ax.set_ylabel("Posterior mean")
ax.legend(fontsize=8)
plt.tight_layout()
plt.show()

Figure 16: True geo-level effect (x) versus posterior mean of the fixed coefficient plus the region slope (y), for TV (blue) and digital (orange). The diagonal is perfect recovery. The horizontal lines are fixed-coefficient means, not the group mean the points should shrink toward.

Hierarchical scale recovery

The estimated between-geo standard deviation tells us how much heterogeneity the model detected.

sigma_tv_post = idata.posterior["sigma_slope_sales_tv"].values.flatten()
sigma_dig_post = idata.posterior["sigma_slope_sales_digital"].values.flatten()

print(
    f"sigma_tv:  posterior mean = {sigma_tv_post.mean():.3f}  (true = {true_sigma_tv})"
)
print(
    f"sigma_dig: posterior mean = {sigma_dig_post.mean():.3f}  (true = {true_sigma_dig})"
)

sigma_tv:  posterior mean = 0.245  (true = 0.15)
sigma_dig: posterior mean = 0.139  (true = 0.1)

Causal queries with geo-varying effects

The do() operator propagates both random intercepts and random slopes, so interventional queries reflect the full heterogeneity across geos.

Average treatment effect

.ate() computes the population-average causal effect by contrasting two intervention levels:

ate_tv = model.ate("sales", "tv", values=(10.0, 30.0))
ate_tv
ate_tv_mean = float(ate_tv.mean())
ate_tv_hdi = np.asarray(ate_tv.hdi()).reshape(-1)
planted_tv_contrast = 0.6 * 20
ate_tv_cover_note = (
    "excludes"
    if not (ate_tv_hdi[0] <= planted_tv_contrast <= ate_tv_hdi[1])
    else "covers"
)

The planted contrast is 0.6 × 20 = 12.0. The printed mean is 10.41, and the 94% interval [9.91, 10.88] excludes 12.0.

Comparing channels

ate_dig = model.ate("sales", "digital", values=(5.0, 20.0))
ate_dig
ate_dig_mean = float(ate_dig.mean())
ate_dig_hdi = np.asarray(ate_dig.hdi()).reshape(-1)
planted_dig_contrast = 0.4 * 15
ate_dig_cover_note = (
    "covers" if ate_dig_hdi[0] <= planted_dig_contrast <= ate_dig_hdi[1] else "excludes"
)

The planted digital contrast is 0.4 × 15 = 6.0. The printed mean is 6.20, and the 94% interval [5.63, 6.78] covers 6.0.

Probability query

p = model.prob("sales > 10", set={"tv": 30.0}, kind="predictive")
print(f"P(sales > 10 | do(tv=30)): {p:.2f}")

Only clang++ is supported. With g++, we end up with strange g++/OSX bugs.
Sampling: [sales]

P(sales > 10 | do(tv=30)): 1.00

Summary

• Adstock and saturation transforms capture the nonlinear response of advertising channels — carry-over across time periods and diminishing returns at high spend levels.
• Posterior predictive checks verify that the model can reproduce observed data patterns before running causal queries.
• do() queries with transforms automatically recompute the full transform pipeline under the intervention, capturing how effects change nonlinearly with spend.
• The marketing funnel is a mediation structure. Upper-funnel spend drives sales both directly and through lower-funnel indicators like search traffic. Conditioning on the mediator blocks the indirect path and underestimates upper-funnel’s contribution.
• Path analysis with labeled coefficients gives the full posterior of direct, indirect, and total effects — the quantities that matter for budget allocation.
• Hierarchical panel models let each region have its own media coefficient while borrowing strength from the full panel through partial pooling.
• Random slopes are not a regional budget. TV intervals separate some regions. Digital intervals overlap, and the slope_sales_digital summary rows are not those effects.

NoteReflection

In your own marketing data:

• Which channels might have indirect effects through downstream indicators like search traffic, app installs, or social engagement? If your current model includes these as independent predictors alongside the spend variables, it may be underestimating the channels that drive them.
• Do channel effects vary by region, customer segment, or time period? Hierarchical models can capture this heterogeneity while borrowing strength where data is sparse.
• Could diminishing returns mean your highest-spend channels are past the point of efficient spending? Response curves from do() queries can reveal this.