pathmc

# pathmc `pathmc` lets you specify a system of structural equations as a string, compile it to a generative PyMC model, fit with MCMC, and reason about causal effects using the do-operator. ## Installation ```bash pip install pathmc # Optional faster samplers (nutpie, numpyro, jax): pip install "pathmc[samplers]" ``` ## Quick start ```python import pathmc spec = """ M ~ a*X Y ~ b*M + c*X indirect := a*b """ m = pathmc.model(spec, data=df) # returns a PathModel (NOT a fitted result) m.fit() # MCMC happens here m.effects_summary() # labeled coefficients + defined params m.ate("Y", "X", values=(0, 1)) # average treatment effect via do() m.adjustment_sets("X", "Y") # valid backdoor adjustment sets m.optimal_adjustment_set("X", "Y") # variance-optimal graphical O-set ``` The DSL is lavaan-inspired: - `Y ~ X` — regression - `Y ~~ X` — residual covariance - `indirect := a*b` — defined parameter - `a*X` — labeled coefficient - Transforms: `adstock(x, decay=...)`, `logistic_saturation(x, lam=...)` - Categorical predictors: string columns are treatment-coded automatically; `C(region, reference='north')` picks the reference and `C(store, prior='hierarchical')` pools one coefficient per level (no reference; the term takes over the intercept's role) ## Decision table | Need | Use | | ------------------------------------------------- | ------------------------------------------------------ | | Build a model from a spec + data | `m = pathmc.model(spec, data=df)` | | Explore the DAG without data | `m = pathmc.model(spec)` (data-free mode) | | Inspect causal DAG | `m.graph()` | | Inspect structural equations + priors | `m.equations()` | | Inspect priors only | `m.priors()` | | Refine priors | `m.set_priors({"beta_Y": Prior(...)})` | | Prior predictive check | `m.sample_prior_predictive()` | | Run MCMC | `m.fit()` | | Summarize posteriors | `m.summary()` or `m.effects_summary()` | | Standardized (stdyx) coefficients | `m.standardized()` | | Path-specific effect (e.g. `X -> M -> Y`) | `m.effect("X -> M -> Y")` | | Posterior predictions | `m.predict(...)` | | **Average treatment effect** | `m.ate(outcome, treatment, values=(0, 1))` | | **Conditional ATE** (effect modification) | `m.cate(outcome, treatment, condition={"Z": z0})` | | ATE on the treated / untreated | `m.att(...)` / `m.atu(...)` | | **Backdoor-adjusted outcome regression** | `adj = m.adjustment_model("X -> Y")` then `adj.fit()` | | Inspect adjustment set / formula before fit | `adj.adjustment_set`, `adj.formula` (before `adj.fit()`) | | Interventional / associational predictions | `m.predictions(outcome, set={...})` | | Interventional contrasts (structural model) | `m.comparisons(outcome, variable, contrast=(0, 1))` | | Marginal slopes under intervention | `m.slopes(outcome, wrt=...)` | | Same interpret API on adjustment model | `adj.comparisons(...)`, `adj.slopes(...)`, etc. | | Probability under intervention | `m.prob("Y > 0", set={"X": 1})` | | Manual intervention | `m.do(set={"X": 1})` | | Intervene on a categorical level | `m.do(set={"region": "south"})`, `m.ate("Y", "region", values=("north", "south"))` | | Declare an integer-coded column as labels | `Y ~ C(store_type)` | | Counterfactual / time-forward (panel) | `m.do(set={...}, simulate_over="time")` | | History-dependent panel treatment policy | `m.do(policy=TreatmentPolicy(...), simulate_over="time", horizon=n_times)` | | Adjustment sets for identification | `m.adjustment_sets(treatment, outcome)` | | Variance-optimal graphical adjustment set | `m.optimal_adjustment_set(treatment, outcome)` | | Yes/no identification check | `m.is_identifiable(treatment, outcome)` | | Front-door identification | `m.frontdoor_identifiable(treatment, mediator, outcome)` | | Flag colliders on any path (heuristic; validity is `is_valid_adjustment_set`) | `m.collider_warnings(adjust, treatment, outcome)` | | Enumerate implied conditional independences | `m.implied_independences()` | | Test DAG implications against data | `m.test_implications()` | | Falsify the whole DAG (permutation test) | `m.falsify()` | | Sensitivity analysis (unmeasured confounding) | `m.sensitivity(outcome, treatment)` | | Placebo refutation of an estimated effect | `m.refute_placebo(outcome, treatment)` | | Simulate from a fully-specified model | `pathmc.simulate(spec, data, params=...)` | ## Gotchas 1. **`pathmc.model(...)` returns a `PathModel`, not a fitted result.** You must call `.fit()` separately. `model()` only parses, builds the DAG, and compiles the PyMC graph — it does not sample. 2. **`m.do(...)` is a structural intervention, not conditioning.** It applies `pm.do()` graph surgery and forward-simulates from the intervened model, propagating posterior uncertainty through the causal chain (g-computation; Robins, 1986). It is **not** the same as conditioning on observed values. For typical user-facing queries, prefer the wrappers `m.ate()`, `m.cate()`, `m.att()`, `m.atu()`, `m.prob()`. 3. **`ate()`/`cate()`/`att()`/`atu()` return an `EstimandResult`, not a `DoResult`.** It knows the outcome, so `r.mean()`, `r.hdi()`, and `r.prob("> 0")` need no variable argument, `float(r)` gives the posterior mean, and printing it shows a tidy summary. `m.do(set=...)` returns a `DoResult` describing the whole system, where accessors still take a variable name (`r.mean("Y")`). `m.do(policy=...)` returns a `PolicyRun`, not a `DoResult`. 4. **The DSL is lavaan-*inspired*, not a 1:1 reimplementation.** `~`, `~~`, `:=`, and labeled coefficients all work. A `:=` name must not match a coefficient label (`a := a * 2` is rejected; use a distinct name such as `twice := a * 2`). Latent-variable measurement models (`=~`) are out of scope in v0.1 — see the user guide for the full operator list. 5. **`Prior` is re-exported from `pymc_extras` for convenience.** `from pathmc import Prior` is a shortcut for `from pymc_extras.prior import Prior`. The canonical reference and list of supported distributions live in `pymc_extras`. 6. **Panel lag terms are declared in the model spec.** Use `lag(sales)` directly in the DSL and pass `panel={"unit": "region", "time": "week"}` to `pathmc.model(...)`. pathmc builds the lagged design internally. 7. **Data-free models have a partial method surface.** When `data=None`, `graph()`, `equations()`, `priors()`, `adjustment_sets()`, `optimal_adjustment_set()`, `is_identifiable()`, `collider_warnings()`, `implied_independences()` all work. `fit()`, `do()`, `ate()`, `cate()`, `design()`, `sample_prior_predictive()`, `test_implications()`, `falsify()`, `sensitivity()`, `refute_placebo()` raise `RuntimeError` until the model is rebuilt with data (and `refute_placebo()` also needs a prior `.fit()`). 8. **`PathModel` is not in `pathmc.__all__`** — it's the class returned by `model()`. You don't import it directly; you receive it. Type annotations can use `pathmc.PathModel` (it is reachable as an attribute) but the public entrypoint is the `model()` function. 9. **`adjustment_model()` returns an `AdjustmentModel` facade.** Inspect `adj.adjustment_set` and `adj.formula` before calling `adj.fit()`. An empty set `{}` is valid when no covariates are needed to block backdoors; if no valid set exists, construction raises (effect not identifiable via backdoor). When several minimal sets exist, pass an explicit set or opt into `adjustment_set="optimal"`; the default does not pick among minimal sets. The O-set can include precision variables that are unnecessary for identification, and its efficiency guarantee is asymptotic for adjusted OLS in causal linear models, not a general Bayesian GLM posterior-variance guarantee. A required latent O-set member and a query with no directed causal path raise explicitly. Pass `data=` when the parent structural model is data-free. Panel models are not supported on the adjustment path. A shared prior key is inherited when its shape does not follow the predictor set or treatment coding: outcome dispersion (`sigma_Y`, `nu_Y`, `alpha_disp_Y`), transform parameters named in the reduced formula, and scalar hierarchical hyperpriors (`mu_beta_*`, `sigma_beta_*`). Otherwise a custom parent value must be passed in `priors=` or construction raises. That includes coefficient priors and non-scalar hyperpriors. Unknown shared keys fail closed the same way. An explicit `priors=` value still wins. 10. **`predictions()` / `comparisons()` / `slopes()` share one API on `PathModel` and `AdjustmentModel`.** On the structural model they use truncated-factorization g-computation; on `adjustment_model()` they delegate to the reduced outcome equation with `estimator="regression_adjustment"`. Read `result.causal`: only queries on the designated treatment support a causal reading; slopes or contrasts on adjustment covariates are interventional on the fitted surface, not causal effects of those covariates. See the user guide page *Predictions, Comparisons, and Slopes*. 11. **String columns become treatment-coded categoricals; integers stay continuous.** `Y ~ region` on a string column creates `beta_Y_region` indexed by level name (coordinate `Y_region_levels`), with the first sorted level as reference. Wrap integer codes in `C(...)` to treat them as labels. The level set is frozen at `model()` time: `do()`, `ate()`, `prob()`, and `predict(data=...)` accept labels, and an unseen label raises rather than re-inferring the coding. Categorical interactions, transforms, outcomes, and panel models are rejected with an explanatory error. 12. **`policy=` is a typed rule, not a Python callback.** `TreatmentPolicy` expressions (`current`, `previous`, `where`, `TimeSchedule`) are evaluated inside a Gaussian lag-panel scan. Pass `simulate_over="time"`, `horizon=` equal to the fitted time count, and an explicit pre-sample baseline. `horizon` is not inferred. `kind="predictive"` is rejected. Latent variables, non-Gaussian families, transforms, categoricals, and residual covariance are rejected rather than run on a weaker executor. 13. **Residual covariance works with Gaussian panel dynamics.** Complete observed Gaussian panels support `~~` together with `lag()` and geometric `adstock()`, including nested transforms. For example, `online ~ adstock(tv, decay=theta_tv) + lag(online); store ~ adstock(tv, decay=theta_tv) + lag(store); online ~~ store` uses `beta_online=[1, 0.8, 0.4]`, `beta_store=[2, 0.5, 0.2]`, `theta_tv=0.6`, and packed `chol_online_store=[1, 0.6, 0.8]` in `simulate()`. Members sort lexicographically, with row-major lower-triangle packing and positive diagonals. Covariance is global; innovations are independent across unit/time/block before temporal propagation. The template omits member `sigma`s; configure their marginal scale priors through `sigma_{member}`. Descendants consume realized states. One-step predictions condition on observed parents/history, while prior prediction, free-running prediction, and predictive hard `do(set=...)` generate correlated trajectories. Hard holds retain unheld members' marginal noise; NaN releases resume their equations. Mean `do()` uses zero process noise. Scaling converts covariance by `D @ Sigma @ D`. Latent nodes, any non-Gaussian equation, missing fitting outcomes, within-block directed dependencies (including cross-member lags), singleton/exogenous endpoints, and `policy=` are excluded. `~~` alone does not identify causal effects. ## Capabilities and boundaries **Agents using pathmc can:** - Write spec strings in the DSL (regressions, residual covariances, defined parameters, labeled coefficients, transforms, categorical predictors via inference or `C(...)`). - Configure custom priors via `Prior` objects from `pymc_extras`. - Run `fit()` with PyMC's NUTS sampler (or `nutpie` / `numpyro` via the `samplers` extra). - Query `ate`/`cate`/`att`/`atu`/`prob`/`effect` with full posterior uncertainty.