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.

What `panel=` does to the frame

panel= is not a boolean. The values are column names:

panel = {"unit": "geo", "time": "week"}

geo and week are coordinates. X and sales stay measurements. The frame you hold stays long: one row is one unit at one time.

For carry-over, pathmc sorts by unit, then time, and reshapes each measurement column into a dense grid of shape (n_times, n_units). A cell is that column’s value for one unit at one time. The result is written back in the original row order. geo and week are not extra axes of that grid. They only place the rows.

This long frame is already in that sort order:

geo	week	X
A	0	0
A	1	1
B	0	0
B	1	0

It becomes this grid for X. Rows are time steps. Columns are units. week is not a column of the grid.

time step	A	B
0	0	0
1	1	0

unit may be several columns, as in panel={"unit": ["geo", "brand"], "time": "week"}. Those columns are joined into one unit id ("North|Acme"). The carry-over grid is still time by unit, not one axis per source column. See Multi-dimensional units.

The dense grid is required for carry-over, not only for lag():

lag() compiles to a panel scan. The lag of a unit is the previous time step of that same unit, which only exists once the rows are on the (n_times, n_units) grid. The scan state at each step is a vector of length n_units.
Geometric adstock() with panel= uses that same scan. The recursion y_t = x_t + decay * y_{t-1} runs inside each unit, including past lag 12.
delayed_adstock() and weibull_adstock() do not use the scan. They still sort, reshape to the same grid, and convolve down the time axis only. Their kernels stay finite.

A panel with no lag() and no convolution adstock (adstock(), delayed_adstock(), or weibull_adstock()) does not take this reshape. Random intercepts, random slopes, and a trend stay row-wise. model() accepts that frame without a dense rectangle. simulate() still requires every unit to share the same times.

On a scan model, the dense grid is a compile-time shape. n_units and n_times are baked in. pm.set_data() itself does not reject a different shape. predict() rejects a scan model whose 2-D data no longer has that (n_times, n_units) shape, rather than silently mis-reshaping. Build a new model on the new frame instead.

do() keeps the compiled horizon. It is not a shape check. pm.set_data() on pymc_model reaches do(kind="mean") when no outcome was observed: a latent model, or a model whose every outcome column has a missing value. A longer grid is then truncated. If any outcome column is complete, that call does not reach do(). do(kind="predictive") on an observed Gaussian model does not read the resized array.

Without panel=, convolution adstock treats row order as one series (the non-panel warning in #655).

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 is retained. It is scored with a pre-sample lag, not 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.

Pre-sample initial state

The first fitted time uses a pre-sample state as its lag() input. That state is not an extra outcome period. Omitting initial keeps the legacy fill: an endogenous lag reads that unit’s first observed value; an exogenous lag reads a <variable>_lag1 column when one is present, otherwise the first observed row, otherwise zeros. Those fills are asymmetric, so an adaptive policy should pass the same explicit baseline to model() and to do().

initial names the base variable used by lag(), not a generated column. A value is a finite number broadcast to every fitted unit, or a mapping of exactly the fitted unit labels. pathmc aligns the mapping by label, not by row order, and scales business units once with the fitted factor. Latent lag bases stay on their init_{variable} prior and are excluded. model.initial_baseline records the normalized model-unit values and whether the fill was explicit. do(initial=) overrides that baseline for one call. It does not refit the model and does not change the stored record. do() without initial reuses the compiled baseline.

initial = {"spend": {"North": 120.0, "South": 80.0}}
model = pathmc.model(
    "sales ~ lag(spend) + trend",
    data=df,
    panel={"unit": "region", "time": "week"},
    initial=initial,
)
model.do(set={"spend": 0.0}, simulate_over="time", initial=initial)

Fitting with random intercepts

Partial pooling adds a hierarchical intercept per unit, borrowing strength across units. The compiler already drops the fixed intercept. Write 0 + so equations() does not print a phantom 1 + intercept that is not in the design matrix:

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

Each region gets its own intercept, drawn from a shared distribution: alpha_{outcome} ~ Normal(mu_alpha, sigma_alpha). Sampling uses a unit-normal offset, alpha_{outcome}_offset, so the sampler does not walk the centered funnel. Posterior tables still report the centered alpha_{outcome}. simulate() takes that centered vector, not the offset.

Random slopes

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

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

Random slopes use the same non-centered form: slope_{outcome}_{var} is the centered effect, and slope_{outcome}_{var}_offset is the unit-normal free parameter.

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 ~ 0 + adstock(tv, decay=theta_tv)",
    data=df,
    panel={"unit": ["geo", "brand"], "time": "week"},
    pooling="partial",
)

That adstock formula takes the carry-over reshape, so the frame must be rectangular: every geo × brand combination shares the same set of week values. A non-rectangular panel raises 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",
)

A NaN in an endogenous array releases that step. The scan runs the variable’s own equation there, and the reported trajectory keeps the finite held values on the other steps, so a temporary hold is “set it, then let it go”. NaN does not release an exogenous input: an input has no equation to fall back on, and those values must be finite.

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.

History-dependent treatment policies

A static set= schedule is fixed before the simulation starts. A treatment policy is a rule evaluated inside the panel scan, so the action at time t can depend on that unit’s simulated history. That is the g-computation path for dynamic regimes: the same posterior draws and the same process noise are reused across rules, and the result keeps the posterior dimension instead of collapsing to a noise-free mean.

This example uses the fitted greenhouse model from Water on a schedule—or when plants need it?, with equations stress ~ lag(stress) + lag(water) + heat and growth ~ stress + water + stress:water. Its pre-sample states are stress and water; growth has no lag and needs no initial value. Here data is the pandas frame used to fit model.

from pathmc import TimeSchedule, TreatmentPolicy, previous, where

n_times = data["time"].nunique()
every_third = TreatmentPolicy(
    "every_third",
    "water",
    TimeSchedule(tuple(1.0 if i % 3 == 0 else 0.0 for i in range(n_times))),
)
responsive = TreatmentPolicy(
    "responsive",
    "water",
    where(previous("stress") > 0.5, 1.0, 0.0),
)
baseline = {"water": 0.0, "stress": 0.0}

fixed = model.do(
    policy=every_third,
    simulate_over="time",
    horizon=n_times,
    initial=baseline,
    random_seed=7,
)
adaptive = model.do(
    policy=responsive,
    simulate_over="time",
    horizon=n_times,
    initial=baseline,
    random_seed=7,
)
contrast = adaptive - fixed
growth = contrast.aggregate_time("growth", reduction="sum")
water = contrast.aggregate_time("water", reduction="sum")

Rules are expressions, not Python callbacks. current(name) reads the simulated value at the decision time and is legal only when that measurement is available before the action. previous(name) reads the preceding period, which at the first fitted time is the explicit pre-sample state. where(condition, yes, no) is the threshold. Comparisons follow the operator exactly: > is false at the threshold. Python if, and, or, and not are rejected.

policy= requires a Gaussian lag panel, simulate_over="time", and an explicit pre-sample baseline covering every observed lag base plus any previous reference the rule reads. do(policy=...) returns a PolicyRun, not a DoResult. kind="predictive" is rejected because the policy path always integrates process noise; omitting random_seed records the seed that was used. horizon is required and must equal the fitted time count; it is not inferred. n_simulations defaults to 256 and must be at least 2. Unit paths are omitted unless keep_trajectories=True. Retained paths use dims (chain, draw, simulation, unit, time).

A fixed TimeSchedule with no residual noise matches do(set=..., kind="mean"). With noise, a threshold on a stochastic state does not match that noise-free path: the reported mean is the average over posterior draws and simulations. Contrast two results only when they share the posterior, the baseline, the seed, the horizon, and the simulation count. The contrast is the difference; it does not pick a winner. simulation_mcse(var) is the per-posterior-draw Monte Carlo standard error. A time-indexed result keeps time; call aggregate_time first for a cumulative or terminal estimand. The chain and draw dimensions are not averaged away.

Latent variables, non-Gaussian families, adstock and other transforms, categoricals, and residual covariance are rejected rather than evaluated by a weaker executor. simulation_info() records the target, the fitted horizon, the unit population, and the seed.

The reference page states the contract. Water on a schedule—or when plants need it? fits the greenhouse panel and contrasts watering every third period with thresholds 0 and 0.5. A causal reading needs sequential exchangeability, consistency, treatment support over relevant histories, model adequacy, and no cross-unit interference. Static backdoor checks do not prove longitudinal identification. A marginal min/max warning does not establish sequential positivity.

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.

The same convention applies to causal estimands and refutation results. cate() interprets its intervention values and condition in business units. att() and atu() interpret their subgroup and intervention values in business units, converting the subgroup value to the fitted scale before selecting observations; these methods currently support cross-sectional models only. refute_placebo() reports both the observed and placebo effects in business units. When do() fills an unspecified exogenous variable from the fitted data, it converts each row back to business units before calculating the fill value for the relevant subgroup.

Convert values explicitly

Use the fitted factors to cross the unit boundary when working directly with scaled data. to_internal() divides business-unit values by their fitted factors; to_business() reverses the conversion. The kind argument identifies the value’s role ("outcome", "regressor", or "coefficient"); "elasticity" and "derived" values are dimensionless and remain unchanged. Repeating a conversion in the same direction on its direct result leaves the values unchanged. Scalar results retain that state through arithmetic with Python int and float values, but round() and NumPy operations can return an untagged value. Arithmetic that changes units also retains the tag, so explicitly cross the boundary before using a converted value in a different unit role.

factors = model.fitted_scaling
internal_channels = factors.to_internal(df, kind="regressor")
business_channels = factors.to_business(internal_channels, kind="regressor")

For a frame, fitted role metadata selects the columns to convert. If you create ScalingFactors manually without roles, select the columns explicitly with dims={"columns": ("tv",)}; otherwise frame conversion raises an error rather than guessing. Unnamed individual values need a column in dims, for example dims={"term": "tv"}. When factors vary by unit, provide row-aligned values and the corresponding data; a single scalar cannot represent different per-unit conversions. Row arrays must end in the fitted row count, and scan-panel arrays must end in (n_times, n_units) order. Ambiguous shapes raise an alignment error instead of silently reordering factors.

Residual covariance with temporal dynamics

Two outcomes can share unmeasured shocks while retaining adstock and autoregressive persistence. For regional weekly sales, use the same TV adstock in both equations and declare a residual block:

spec = """
online ~ adstock(tv, decay=theta_tv) + lag(online)
store ~ adstock(tv, decay=theta_tv) + lag(store)
online ~~ store
"""
params = {
    "beta_online": [1.0, 0.8, 0.4],
    "beta_store": [2.0, 0.5, 0.2],
    "theta_tv": 0.6,
    "chol_online_store": [1.0, 0.6, 0.8],
}
# grid contains region, week, and tv; one row per region/week.
simulated = pathmc.simulate(
    spec,
    grid,
    params=params,
    panel={"unit": "region", "time": "week"},
    random_seed=42,
)
m = pathmc.model(spec, simulated, panel={"unit": "region", "time": "week"})
m.fit()

simulate_params_template(spec, grid, panel=...) reports precisely these parameter names and shapes. The coefficient order is intercept, adstock, self-lag. Block members have no separate sigma_online or sigma_store simulation parameter: chol_online_store supplies both scales and correlation. Members are sorted by name; for two outcomes the packed factor is [L00, L10, L11]. Here L = [[1, 0], [0.6, 0.8]], so Sigma = L @ L.T = [[1, 0.6], [0.6, 1]]. Supply a finite packed vector with positive diagonal entries. A connected component such as a ~~ b; b ~~ c is one full three-outcome covariance block.

The covariance is global, shared across units and periods. Each unit/time receives a fresh independent white innovation vector, colored by L; there is no shared regional shock or additional serial error process. Temporal dependence comes from the model’s dynamics. Generated states already contain their correlated innovations when descendants and later lags read them. A descendant such as total ~ 0 + 1*online + 1*store therefore consumes those same realized sales. Coefficient pooling changes structural means while leaving the covariance global.

Fitting evaluates one joint Gaussian density per block using observed contemporaneous parents and observed lag history. predict(one_step_ahead=True) uses those conditional means and jointly samples the block outcomes under their original names. Prior prediction and predict(one_step_ahead=False) generate coherent free-running trajectories. do(kind="predictive", simulate_over="time") does the same under hard interventions; holding one member leaves the other members’ marginal innovation laws intact. A NaN step releases a held endogenous outcome back to its equation using the previously carried state. do(kind="mean") instead follows noise-free structural recursion, independent of covariance parameters for fixed mean parameters.

Factors and parameters are in model units. With outcome scaling matrix D, the business-unit innovation covariance is D @ Sigma @ D; generated outcomes are unscaled once. Configure a member’s marginal scale prior through sigma_{member} in priors= or .set_priors(). These priors feed the LKJ covariance, rather than adding another scale RV. chol_{block} is a simulation parameter and an introspection entry, not a Prior override key.

Scan + ~~ requires complete finite observed Gaussian outcomes in every equation. Latent nodes, non-Gaussian equations anywhere, and missing fitting outcomes remain unsupported in this combination. Residual endpoints must each have an equation, and each block must contain at least two distinct outcomes. Directed dependencies within a block, including cross-member lags and transforms, remain rejected; self-lags are supported. Prediction uses the fitted panel layout, and policy= remains unsupported with ~~. Declaring residual covariance represents shared unexplained variation; it does not itself identify causal effects.

Cross-sectional and non-scan panel residual blocks also propagate realized noisy members to descendants, with the same within-block directed-edge restriction.

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()
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"},
)

What panel= does to the frame

Temporal dependencies with lag()

Pre-sample initial state

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

History-dependent treatment policies

Scaling heterogeneous units

Convert values explicitly

Residual covariance with temporal dynamics

Latent dynamics with sparse measurements

What `panel=` does to the frame

Temporal dependencies with `lag()`

Structured pooling with `by_var`