Panel Data

pathmc supports panel (longitudinal) data — repeated observations of multiple units over time. Panel structure enables modeling temporal effects like carry-over and lagged dependencies, and borrowing strength across units via hierarchical pooling.

Temporal dependencies with `lag()`

Use the lag(var) term directly in your formula to include lag-1 values within each unit:

model = pathmc.model(
    "sales ~ lag(spend) + trend",
    data=df,
    panel={"unit": "region", "time": "week"},
)

The lag() term is a structural part of the formula, not a data preprocessing step. pathmc computes lag-1 values from the previous time step within the same unit — never from a different unit. The first time step of each unit has no valid lag and is dropped.

When you call model.graph(), lagged terms appear as dashed edges labeled with their coefficient and (t−1). A self-lag like lag(Y) in Y ~ lag(Y) renders as a dashed self-loop on Y. A cross-variable lag like lag(X) in Y ~ lag(X) renders as a dashed edge from X to Y, alongside any contemporaneous edge.

Figure 1: Lag-1 temporal structure respects panel boundaries. Dashed arrows connect consecutive time steps within each unit. No arrows cross unit boundaries — Region B’s week 1 does not leak into Region A’s lag.

Fitting with random intercepts

Partial pooling adds a hierarchical intercept per unit, borrowing strength across units:

model = pathmc.model(
    "sales ~ lag(spend) + trend",
    data=df,
    panel={"unit": "region", "time": "week"},
    pooling="partial",
)

Each region gets its own intercept, drawn from a shared distribution: alpha_region ~ Normal(mu_alpha, sigma_alpha).

Random slopes

For varying effects across units, specify which predictors get random slopes:

model = pathmc.model(
    "sales ~ spend + lag(spend)",
    data=df,
    panel={"unit": "region", "time": "week"},
    pooling={"intercept": True, "slopes": ["spend"]},
)

Multi-dimensional units

Marketing and longitudinal data often have more than one identifier per row — geo × brand, country × channel, store × product. Pass a list of columns as unit and pathmc pools over their Cartesian product:

model = pathmc.model(
    "sales ~ adstock(tv, decay=theta_tv)",
    data=df,
    panel={"unit": ["geo", "brand"], "time": "week"},
    pooling="partial",
)

The data must be rectangular: every geo × brand combination shares the same set of week values. Non-rectangular panels raise an error naming the offending combinations. Under the hood the composite key ("North|Acme" style labels) becomes the unit coordinate, so each cell gets its own random intercept under pooling="partial" — the same machinery as single-column panels. Because | separates the source labels, values in multi-dimensional unit columns cannot themselves contain |; recode those values before building the model. Per-dimension coords such as geo keep their native dtype (numeric columns stay numeric), while the composite unit labels are strings — cast when joining lookups across the two.

Structured pooling with `by_var`

Different parameters often need different pooling structure — a channel coefficient that varies by geo, an adstock decay that is fully unpooled per cell. The by_var entry of pooling controls this per parameter:

model = pathmc.model(
    "sales ~ 0 + tv + radio",
    data=df,
    panel={"unit": ["geo", "brand"], "time": "week"},
    pooling={
        "intercept": True,
        "by_var": {
            # coefficient of 'tv' varies per cell, pooled across geos:
            # beta_tv[c] ~ Normal(mu_tv_geo[geo(c)], sigma_tv_geo)
            "tv": {"coefficient": ("geo",)},
            # adstock decay is a separate parameter per panel cell
            "theta_tv": "none",
        },
    },
)

The tuple in "coefficient" names panel dimensions (the columns in panel['unit']) to pool over; a single string works for one dimension. Each pooled predictor gets a per-cell RV (beta_{var}) plus hyperpriors named after the requested dims (mu_{var}_{dim}, sigma_{var}_{dim}), and it replaces that predictor’s flat entry in beta_{lhs}. A bare string "none" forces unpooled per-cell parameters instead — for a predictor this means one free coefficient per cell with no hierarchy, and for a transform parameter (e.g. decay=theta_tv) it upgrades the shared scalar to one value per cell.

Use multiple names to pool over a dimension combination, for example "tv": {"coefficient": ("geo", "brand")} creates mu_tv_geo_brand[geo, brand]. If the observed panel omits some combinations, pathmc emits a compile-time warning: those unused hyperprior cells remain prior-only and can look like estimates in the posterior. Use a rectangular panel when every cell should be informed by data.

Posterior predictive checks with a lagged outcome

When the outcome depends on its own past — sales ~ lag(sales) — there are two different things predict() could mean, so it takes a flag.

By default (one_step_ahead=True) each time step conditions on the observed previous value of the outcome, matching the likelihood the model was fitted with. This is the posterior predictive distribution you want for a model check: residuals are one-step-ahead errors, and comparing them against the observed series is a fair test.

Passing one_step_ahead=False lets the recursion run free: step t conditions on the model’s own simulated value at t−1. Uncertainty compounds over the series, which is the honest picture of multi-step forecast error but a misleading basis for a fit diagnostic.

ppc = model.predict()  # one-step-ahead, for model checking
traj = model.predict(one_step_ahead=False)  # free-running trajectories

Note

The flag only matters for panel models with a lagged endogenous term. Everywhere else the two are identical and the argument is ignored. do() is unaffected — interventional simulation always runs the free-running recursion.

Time-forward simulation

When a model includes temporal state — adstock transforms or lagged variables — the do() operator must walk through time step by step to correctly re-compute the temporal dynamics under the intervention. Activate this with simulate_over="time":

scenario = model.do(
    set={"spend": 120},
    simulate_over="time",
    kind="mean",
)

Time-varying interventions

The set= parameter accepts arrays of shape (n_times,) for per-time-step interventions, enabling scenarios like a temporary spend increase:

import numpy as np

spend_scenario = np.full(n_weeks, 100.0)
spend_scenario[10:20] = 150.0  # boost during weeks 11-20

result = model.do(
    set={"spend": spend_scenario},
    simulate_over="time",
)

Per-time-step results

Panel do() results include per-time-step data accessible via .by_time():

contrast = scenario - baseline
curve = contrast.by_time("sales")  # shape: (n_times, n_samples)
mean_by_week = curve.mean(axis=1)  # posterior mean at each week
time_labels = contrast.time_index  # the time column values

When is simulate_over="time" needed?

Not always! If your panel model has only random intercepts, random slopes, or trend terms — but no adstock transforms or lagged variables — a regular do() call works fine. See Time-Forward Panel Simulation for the full explanation.

See the Panel Data Models, Difference-in-Differences, and Media Mix Models examples for panel mode in action.

Scaling heterogeneous units

When geos (or brands, or countries) differ in magnitude by orders of magnitude — a national channel vs. a regional one — a single pooled coefficient has to bridge wildly different scales, and default priors calibrated for one unit are badly suited to another. The scaling= argument divides affected columns by fitted scale factors before the model is compiled, so estimation happens on a common internal scale.

import xarray as xr

# population per geo, used as an external divisor for the channel column
population = xr.DataArray([10.0, 1.0], coords={"geo": ["north", "south"]}, dims=["geo"])

model = pathmc.model(
    spec,
    data=df,
    panel={"unit": "geo", "time": "week"},
    scaling=pathmc.Scaling(
        # divide each outcome by its per-geo maximum
        target={"method": "max", "dims": ("geo",)},
        # divide spend by population (like pymc-marketing's FixedScaling)
        channel={"method": "divide", "by": population, "dims": ("geo",)},
    ),
)

Each role takes a spec with a "method":

"max" / "mean": divide by the group-wise maximum/mean of that same column (groups defined by "dims", which must name panel["unit"] columns; omit "dims" for one global scale).
"fixed": divide by a supplied constant or grid ("value": ...).
"divide": divide by an external grid keyed by the unit dims ("by": ...) — an xarray.DataArray whose coordinates carry the dim names, or a dict mapping unit labels (or label tuples) to divisors.

The correspondence with pymc_marketing.mmm.scaling is direct: the target / channel slots mirror pymc-marketing’s identically-named MMM scaling slots, "max" with dims= matches Scaling(method="max", dims=...), and "divide" / "fixed" with a grid match FixedScaling(values=<DataArray>). pathmc just applies them to any panel outcome or predictor, not only MMM targets and media channels.

Fitted factors live on the model as model.fitted_scaling. When you simulate from known parameters that are expressed in these scaled units, pass the scaling object (or the fitted factors themselves) to simulate(): exogenous columns are divided before compilation and every generated outcome is multiplied back into business units.

sim = pathmc.simulate(
    spec,
    data=df,
    params=params_in_scaled_units,
    panel={"unit": "geo", "time": "week"},
    scaling=model.fitted_scaling,  # reuse the exact estimation-time scales
)

Omitting scaling= in simulate() changes nothing: parameters are then interpreted directly in raw data units.

User-facing inputs and outputs use business units: do(set=) values are divided by the fitted factor (including per-unit factors) before graph surgery, so do(set={"tv": 500}) means 500 of the original column, and predict(), do(), effects_summary(), and effect() return outcomes and coefficients in business units. predictions(newdata=) grids are also business units (the grid is divided before compilation, matching do(set=)). Coefficient rescaling depends on the term: a linear or adstock slope is f_out / f_pred; a saturating or HSGP coefficient is f_out only (the regressor is unitless); an interaction is f_out / prod(f_pred_i). Defined parameters (:=) are evaluated from those rescaled labeled draws, so indirect := a*b agrees with the product of the a and b rows. When factors vary by unit, a reported coefficient uses a data-weighted mean of the per-unit divisors — a ratio of means, not the mean of ratios, so it is not any one unit’s coefficient. Estimation still runs on scaled columns internally; model.fitted_scaling exposes the divisors for simulate() and manual transforms.

Residual covariance simulation

Cross-sectional residual-covariance blocks (Y1 ~~ Y2) are realized inside the generative graph, so descendants such as Z ~ Y1 receive the noisy block draws rather than mu_Y1. Directed edges within the same ~~ block (Y2 ~ Y1 together with Y1 ~~ Y2) still raise NotImplementedError, as do scan-compiled panel models that combine simulate() with residual covariances.

Latent dynamics with sparse measurements

Some quantities are never observed directly at every time step — brand awareness, true demand, sentiment. pathmc supports latent variables in panel models: declare them with latent=[...] and anchor them with a sparse measurement equation whose rows are NaN wherever no survey (or other measurement) exists:

model = pathmc.model(
    """
    survey ~ 0 + 1*awareness
    awareness ~ lag(awareness)
    """,
    data=df,  # 'survey' has NaN outside surveyed weeks
    panel={"unit": "market", "time": "week"},
    latent=["awareness"],
    families={"awareness": "latent_normal"},
)

The second line makes the latent state an AR(1) process: each week’s awareness is a coefficient times the previous week’s state plus process noise (families={"awareness": "latent_normal"}). Because latent variables have no data column, their initial condition cannot be read off the data — instead it is estimated: an init_awareness free parameter (one value per unit, default prior Normal(0, 1)) seeds the recursion. It is the state preceding the first observed period — the first latent value is computed as intercept + phi * init_awareness plus process noise. The default Normal(0, 1) prior is not scale-free — override it when the latent is not O(1), otherwise t=0 is pinned near zero and early-period trajectories are biased. Override it like any other prior when you have prior knowledge about the starting level:

from pathmc.priors import Prior

model = pathmc.model(
    spec,
    data=df,
    panel=panel,
    latent=["awareness"],
    families={"awareness": "latent_normal"},
    priors={"init_awareness": Prior("Normal", mu=0.5, sigma=0.2)},
)

After fitting, extract the inferred trajectory with latent_trajectory():

model.fit(draws=1000, tune=1000, chains=4)
traj = model.latent_trajectory("awareness")
traj.dims  # ('chain', 'draw', 'time', 'unit')
mean_traj = traj.mean(dim=("chain", "draw"))  # (time, unit) posterior mean

The result is an xarray.DataArray with dimensions (chain, draw, time, unit) and coordinates taken from the model’s panel structure. For stochastic latents it reads the realized state registered under the variable’s own name; for deterministic latents it reads the mu_{var} deterministic.

The same machinery works for forward simulation. pathmc.simulate() treats init_{var} as a required parameter alongside coefficients and scales — use pathmc.simulate_params_template() to discover every name and shape:

params = {
    "beta_awareness": [0.0, 0.9],  # intercept, AR coefficient
    "sigma_awareness": 0.15,
    "init_awareness": [0.5, 0.4],  # one per unit
    "sigma_survey": 0.05,
    "survey_unobserved": 0.0,
}
sim = pathmc.simulate(
    spec,
    data=df,
    params=params,
    panel=panel,
    latent=["awareness"],
    families={"awareness": "latent_normal"},
)

Temporal dependencies with lag()

Fitting with random intercepts

Random slopes

Multi-dimensional units

Structured pooling with by_var

Posterior predictive checks with a lagged outcome

Time-forward simulation

Time-varying interventions

Per-time-step results

Scaling heterogeneous units

Residual covariance simulation

Latent dynamics with sparse measurements

Temporal dependencies with `lag()`

Structured pooling with `by_var`