Predictions, Comparisons, and Slopes

ate() answers a focused question: what is the average difference in Y between two intervention levels of a treatment? That wrapper is enough for many binary-treatment studies, but real analyses often need response-scale predictions on a covariate grid, ratio or lift contrasts, local slopes for continuous treatments, or unit-level draws you can summarize yourself. This page introduces the shared interpret layer on both PathModel and AdjustmentModel: predictions(), comparisons(), and slopes(). Every result carries causal metadata so you can see whether a query is interventional, identifiable, and safe to read as a treatment effect. The terminology and workflow build on prior work in the R ecosystem, especially marginaleffects, which provides a general framework for predictions, comparisons, and slopes. Bambi’s interpret API offers related tools for Bayesian regression models and draws inspiration from marginaleffects; pathmc adapts these ideas to its structural and adjustment-model APIs.

Plots and worked figures live in the Examples gallery, not here:

Why ate() is not enough

Average treatment effects are one estimand in a larger vocabulary. The following situations show up often once you move beyond a single binary contrast.

Continuous treatments. You may want a local slope of Y with respect to X near observed covariate values, not only a two-point difference between fixed levels.
Non-additive contrasts. Uplift ((Y_hi - Y_lo) / Y_lo) and ratio (Y_hi / Y_lo) targets are common in marketing and epidemiology; they are not the same as a difference on the response scale.
Covariate grids. You may need predicted response means at many (X, Z) combinations, then contrast or average in a second step.
Unit-level uncertainty. ate() collapses to one draw per MCMC iteration; keeping the unit dimension lets you build subgroup summaries without re-running inference.

comparisons() generalizes the two-point contrast behind ate(). predictions() and slopes() cover grids and local derivatives. All three share the same draw contract and metadata fields.

The draw contract

Every interpret query returns posterior response-mean draws: the expected value of the outcome on the model’s link scale, after any inverse link, with residual noise turned off (the same kind="mean" semantics as do()).

Internal storage is an xarray dataset with dimensions chain, draw, and usually unit. Each row in the fitted data (or each row of newdata) is one unit. The workflow is always contrast first, then average: build unit-level contrasts or slopes, then call .mean("unit") when you want a population summary. HDIs and means are computed on those draws, not by averaging interval endpoints from separate fits.

Do not average endpoints

If you need an average contrast, pass average_by="all" or average the unit-level draws yourself. Do not fit separate models per grid point and average HDI bounds; that breaks the draw-wise logic this layer enforces.

With average_by="all", comparisons() and slopes() return an EstimandResult (one value per chain/draw). With average_by=None, they return an InterpretResult that keeps the unit dimension. predictions() always returns InterpretResult.

Shared example

We use a fork confounder DAG throughout: income affects both ad_spend and converted, and ad_spend affects converted. The same dataset keeps the snippets runnable in CI.

Figure 1: Fork confounder DAG. income confounds ad_spend and converted; ad_spend is the designated treatment.

import numpy as np
import pandas as pd
import pathmc

rng = np.random.default_rng(42)
n = 300
income = rng.normal(size=n)
ad_spend = 0.5 * income + rng.normal(scale=0.5, size=n)
converted = (0.4 * ad_spend + 0.6 * income + rng.normal(scale=0.5, size=n) > 0).astype(
    float
)
df = pd.DataFrame({"income": income, "ad_spend": ad_spend, "converted": converted})

spec = """
ad_spend ~ income
converted ~ ad_spend + income
"""
m = pathmc.model(spec, data=df, families={"converted": "bernoulli"})
m.fit(random_seed=42)

adjusted = m.adjustment_model("ad_spend -> converted")
adjusted.fit(random_seed=42)

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [beta_converted, sigma_ad_spend, beta_ad_spend]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 3 seconds.
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [beta_converted]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 2 seconds.

Identification for the total effect of ad_spend on converted requires adjusting for income; see adjustment_sets() on the causal inference page rather than re-deriving the backdoor argument here.

predictions()

predictions(outcome, set=..., newdata=...) evaluates the response mean for each unit in the fitted frame or on a covariate grid.

On PathModel, omitting set yields associational predictions: the fitted structural equations evaluated at observed covariates, with no graph surgery. Passing set={"ad_spend": 1.0} applies do(ad_spend=1) while holding other endogenous mechanisms fixed at each unit’s covariate profile.

On AdjustmentModel, the same signature applies to the reduced outcome equation. Without set, predictions are associational fits of converted ~ ad_spend + income. With set={"ad_spend": 1.0}, pathmc standardizes the outcome at the designated treatment value; causal=True only when the intervention targets that treatment alone.

Build grids with datagrid() (re-exported from pathmc and available on both model types). pathmc.datagrid(df, ...) uses the units in the frame you pass. On a scaled fitted model, m.datagrid(...) returns business units too: supplied values are left unchanged and columns held at their fitted mean or mode are converted from the internal scale.

grid = pathmc.datagrid(df, ad_spend=[0.0, 1.0], income=[-1.0, 0.0, 1.0])
pred_grid = m.predictions("converted", newdata=grid)

Associational baseline on the structural model:

pred_obs = m.predictions("converted")
pred_obs.interventional, pred_obs.causal

(False, False)

Interventional prediction at ad_spend=1 on the adjustment model:

pred_do = adjusted.predictions(set={"ad_spend": 1.0})
pred_do.estimator, pred_do.causal

('regression_adjustment', True)

See Conditional Predictions for response-scale curves, HDI bands, and datagrid() plots.

comparisons()

comparisons(outcome, variable, contrast=(lo, hi), comparison=..., conditional=..., average_by=...) contrasts two intervention values of variable under do(), holding any conditional covariates fixed.

The comparison argument selects the estimand applied draw-wise at each unit before averaging.

`comparison`	Unit-level contrast	Typical use
`"diff"`	`Y_hi - Y_lo`	Average treatment effect (same target as ate())
`"ratio"`	`Y_hi / Y_lo`	Multiplicative lift on the response scale
`"lift"`	`(Y_hi - Y_lo) / Y_lo`	Percent change relative to the baseline level

With comparison="diff", contrast=(0.0, 1.0), and average_by="all", comparisons() matches ate() on the same contrast. The convenience wrapper and the general method share one implementation path.

comp = m.comparisons(
    "converted", "ad_spend", contrast=(0.0, 1.0), comparison="diff", average_by="all"
)
ate = m.ate("converted", "ad_spend", values=(0.0, 1.0))
np.isclose(comp.mean(), ate.mean())

np.True_

The adjustment model exposes the same call shape; outcome and treatment default to the designated ad_spend -> converted query.

adj_comp = adjusted.comparisons(comparison="diff", average_by="all")
adj_ate = adjusted.ate(values=(0.0, 1.0))
np.isclose(adj_comp.mean(), adj_ate.mean())

np.True_

Set average_by=None when you need unit-level InterpretResult draws (for example, to inspect heterogeneity before averaging). See Interventional Contrasts for grid contrasts, ratio/lift tables, and unit-level distributions.

slopes()

slopes(outcome, wrt, slope="dydx", eps=..., conditional=..., average_by=...) approximates a local derivative by a symmetric finite difference: intervene at each unit’s observed wrt value and at wrt + eps, then divide the response contrast by eps. The default slope="dydx" is an ordinary derivative on the response scale. Elasticity-style summaries are available via slope="eyex" (percent change in Y per percent change in wrt); the other options (eydx, dyex) mix level and percent changes when you need them.

Local positivity matters: if wrt + eps leaves the support of the data, the finite difference can extrapolate far outside the region where the fit is informed.

Structural model, treatment slope (matches an interventional local effect when the DAG is linear):

slope_x = m.slopes("converted", "ad_spend", average_by="all")
slope_x.mean(), slope_x.causal

(0.09579199379829599, True)

Adjustment model defaults wrt to the designated treatment:

adj_slope = adjusted.slopes(average_by="all")
adj_slope.causal, adj_slope.estimator

(True, 'regression_adjustment')

See Local Slopes on the Response Scale for response-scale derivatives on GLMs and the adjustment-covariate warning.

Two estimators, one API

The method names and arguments are identical; only the estimator backend and causal stamping differ.

# Structural g-computation (truncated factorization on the full DAG)
m.comparisons("converted", "ad_spend", contrast=(0.0, 1.0), average_by="all")

# Outcome-regression standardization on the reduced equation
adjusted.comparisons(contrast=(0.0, 1.0), average_by="all")

On PathModel, draws propagate through every structural equation you specified. On AdjustmentModel, draws come from the single reduced outcome regression identified by the backdoor set, with estimator="regression_adjustment" stamped on the result. When both models are correctly specified for the total effect, their comparisons(..., comparison="diff") answers should agree; when they diverge, the mechanisms differ even though the API does not. See the worked SCM vs Regression Adjustment page for a side-by-side ATE comparison on one DAG.

Causal metadata on every result

Each InterpretResult and stamped EstimandResult exposes metadata you should read before reporting a number as causal.

Field	Meaning
interventional	`True` when `set` / do() graph surgery was used
causal	`True` when the query supports a causal reading of the designated estimand
identifiable	Backdoor identifiability for the treatment-outcome pair when defined
estimator	`"structural"` on PathModel, `"regression_adjustment"` on AdjustmentModel

On PathModel, any interventional query on the treatment-outcome pair sets causal=True when identification holds. On AdjustmentModel, causal=True only when the varied variable is the designated treatment from adjustment_model("ad_spend -> converted"). Intervening on or taking slopes with respect to a covariate in the adjustment set is still interventional on the reduced equation, but it is not the total effect of that covariate.

Adjustment covariates are not treatments

A slope or comparison with respect to income on an AdjustmentModel answers a what-if on the fitted outcome surface; it does not identify the causal effect of income on converted. Check result.causal before narrating the number as an effect of the variable you varied.

Confounder slope on the adjustment model (interventional fit, not a causal effect of income):

slope_z = adjusted.slopes(wrt="income", average_by="all")
slope_z.causal, slope_z.interventional

(False, True)

Next steps

Conditional Predictions, Interventional Contrasts, and Local Slopes for figure-first worked examples.
Causal Inference for the do-operator, ate() / cate() / att() / atu(), and identification helpers.
Estimation Approaches for when to fit the full structural model versus a reduced adjustment equation.
SCM vs Regression Adjustment for a worked comparison of both estimators on one dataset.
Treatment Effects with Non-Linear Models for why a GLM coefficient is not the response-scale effect.
API reference for PathModel, AdjustmentModel, InterpretResult, and EstimandResult.

Reflection

Pick a continuous treatment in your domain (dose, price, hours of training). Would a two-point ate() between fixed levels suffice, or do you need a grid of predictions(), a ratio contrast, or a local slopes() summary? If you fit an adjustment_model(), which covariates in the formula are adjustment variables rather than treatments you could ethically intervene on?