Causal Inference

Causal Inference

Dependencies

This skill requires the bayesian-workflow skill for all PyMC modeling steps (priors, sampling, diagnostics, calibration, reporting).

Detect it:

ls ~/.claude/skills/bayesian-workflow/SKILL.md 2>/dev/null || ls .claude/skills/bayesian-workflow/SKILL.md 2>/dev/null

If not found, install it:

git clone https://github.com/Learning-Bayesian-Statistics/baygent-skills.git /tmp/baygent-skills
cp -r /tmp/baygent-skills/bayesian-workflow ~/.claude/skills/

For all PyMC modeling steps (priors, sampling, diagnostics, calibration, reporting), follow the bayesian-workflow skill.

Workflow overview

Every causal analysis follows this sequence. Steps 1-4 are the thinking phase (no code). Steps 5-8 are the doing phase. Think before you do.

1. Formulate the causal question — Propose precise estimand (ATE, ATT, LATE, etc.). ⚠️ ASK USER TO CONFIRM. 2. Draw the DAG — Propose causal graph with nodes, edges, and explicit non-edges. ⚠️ ASK USER TO CONFIRM. See references/dags-and-identification.md 3. Identify — Determine identification strategy (backdoor, front-door, IV, RDD, DiD). ⚠️ ASK USER TO CONFIRM untestable assumptions. See references/dags-and-identification.md 4. Choose design — Match problem to method using table below. ⚠️ ASK USER TO CONFIRM. See references/quasi-experiments.md or references/structural-models.md 5. Estimate — Build and fit the model. Delegate all PyMC mechanics to bayesian-workflow skill. 6. Refute — MANDATORY. Run design-specific robustness checks. See references/refutation.md 7. Interpret — Effect size + decision-relevant HDIs + probability of direction. 8. Report — Generate causal analysis report. See references/reporting.md

Design selection guide

Critical rules

• No estimation without a confirmed DAG. A causal graph is not optional decoration — it makes

assumptions explicit and determines the adjustment set. If the user resists, explain why the DAG is non-negotiable before proceeding.

• No causal claims without refutation. Every design has failure modes. Run at minimum one

design-specific robustness check (placebo test, sensitivity analysis, falsification test) before reporting results. See references/refutation.md.

• State assumptions before results. Lead with what must be true for the estimate to be causal.

Bury the estimate after the assumptions, not before. This is not optional politeness — it prevents misuse of results.

• Adapt HDIs to the decision context. The bayesian-workflow skill's 94% HDI is a sensible

default; adapt it with explicit explanation when the decision stakes warrant it (e.g., 89% for exploratory, 97% for high-stakes policy). Report multiple intervals when the decision threshold matters.

• Downgrade causal language when warranted. If identification assumptions are unverifiable or

refutation raises flags, soften claims: "consistent with a causal effect" not "causes", "estimated effect" not "true effect". Flag uncertainty loudly in the report.

• Ask the user when domain knowledge is needed. You cannot know whether an instrument is valid,

whether parallel trends holds, or whether a confounder exists without domain expertise. Ask before assuming.

• Delegate PyMC mechanics to bayesian-workflow. This skill handles causal structure and design.

The bayesian-workflow skill handles priors, sampling, diagnostics, calibration, and reporting format. Don't duplicate those rules here.

Common gotchas

These are battle-tested lessons that save hours of debugging:

• CausalPy formula syntax uses `C()` for categoricals. Passing a string column directly without

C() will silently produce wrong dummy coding. Always wrap categorical treatment and group variables: "y ~ C(treatment) + C(group)".

• DoWhy requires explicit `U` nodes for unobserved confounders. Omitting them from the graph

will make DoWhy treat your model as fully identified when it isn't. Add latent nodes explicitly and mark them as unobserved.

• CausalPy's PyMC models don't auto-store log-likelihood. Same issue as bayesian-workflow:

nutpie silently drops it. Call pm.compute_log_likelihood(idata, model=model) after sampling if you need it for model comparison.

• Parallel trends is untestable in the post-treatment period. Pre-treatment trend tests are

necessary but not sufficient — passing them doesn't prove the assumption holds after treatment. State this explicitly in every DiD report.

• Synthetic control requires the treated unit to lie within the convex hull of donors. If the

treated unit is an outlier (highest GDP, largest city), no weighted combination of donors can approximate its counterfactual. Check this before running — if violated, the design is invalid.

• DiD group variable must be dummy-coded (0/1). CausalPy rejects string labels like "treatment"/"control". Use integers: 1 = treatment, 0 = control. Data also requires a unit column.
• SyntheticControl expects wide-format data. Index = time, columns = unit names, values = outcome. If your data is long format, pivot first: df.pivot(index="date", columns="unit", values="outcome").

When things go wrong

causal-inference

An opinionated Agent Skill for production-grade Bayesian causal inference using PyMC, CausalPy, and DoWhy.

Compatible with Claude Code, Kimi Code, Cursor, Gemini CLI, and any agent that supports the Agent Skills spec.

Full breakdown here.

What it does

Guides your coding agent through the full causal inference workflow, enforcing DAG-first thinking and mandatory assumption checkpoints:

1. Formulate the causal question — Precise estimand (ATE, ATT, LATE, etc.) 2. Draw the DAG — Explicit causal graph with nodes, edges, and non-edges 3. Identify — Backdoor, front-door, IV, RDD, DiD 4. Choose design — Match problem to method (DiD, synthetic control, ITS, RDD, IV, IPSW, structural) 5. Estimate — Build and fit the model (delegates PyMC mechanics to bayesian-workflow) 6. Refute — Mandatory design-specific robustness checks (placebo tests, sensitivity analysis, falsification) 7. Interpret — Effect size + decision-relevant HDIs + probability of direction 8. Report — Defensible causal language with assumption-first structure

The skill enforces guardrails that agents won't apply on their own: no estimation without a confirmed DAG, no causal claims without refutation, assumptions stated before results, and automatic downgrading of causal language when warranted.

Install

This skill requires the bayesian-workflow skill for all PyMC modeling steps (priors, sampling, diagnostics, calibration, reporting). Install both together.

Claude Code

git clone https://github.com/Learning-Bayesian-Statistics/baygent-skills.git /tmp/baygent-skills
mkdir -p ~/.claude/skills
cp -r /tmp/baygent-skills/bayesian-workflow ~/.claude/skills/
cp -r /tmp/baygent-skills/causal-inference ~/.claude/skills/

For project-level installation (available only in that project), copy into .claude/skills/ at the project root instead.

Other compatible agents (Kimi Code, Cursor, etc.)

git clone https://github.com/Learning-Bayesian-Statistics/baygent-skills.git /tmp/baygent-skills
cp -r /tmp/baygent-skills/bayesian-workflow/ ~/.config/agents/skills/bayesian-workflow/
cp -r /tmp/baygent-skills/causal-inference/ ~/.config/agents/skills/causal-inference/

Python dependencies

mamba install -c conda-forge pymc nutpie arviz arviz-stats causalpy dowhy

Example prompts

Once installed, just ask your agent naturally:

• "We ran a marketing campaign in 3 cities starting in March. I have monthly revenue data for those cities plus 10 control cities. Did the campaign work?"
• "I have observational data on a drug treatment. Help me estimate the causal effect controlling for confounders."
• "Does X cause Y? I have panel data with a treatment that rolled out at different times across regions."
• "I need to estimate the effect of a policy change using regression discontinuity — students above a test score threshold got a scholarship."
• "Build a synthetic control for California's tobacco tax using other states as donors."

Design selection guide

What's included

causal-inference/
├── SKILL.md                          # Main workflow instructions
└── references/
    ├── dags-and-identification.md    # DAG construction, backdoor/front-door criteria
    ├── quasi-experiments.md          # DiD, synthetic control, ITS, RDD, IV, IPSW
    ├── structural-models.md          # pm.do(), pm.observe(), counterfactuals
    ├── refutation.md                 # Design-specific robustness checks
    └── reporting.md                  # Causal language guardrails, report templates

License

MIT - see LICENSE.

DAGs and Causal Identification

Contents

1. Drawing causal graphs 2. DAG specification in code 3. Common causal structures 4. Identification strategies 5. Finding adjustment sets 6. Collider bias warning 7. Asking the user: prompt templates

---

Drawing causal graphs

A Directed Acyclic Graph (DAG) represents causal structure with nodes (variables) and directed edges (direct causal effects). Arrows point from cause to effect.

The most important rule: missing edges are the STRONGEST assumptions. Omitting an edge asserts that one variable has no direct causal effect on another — a strong claim that must be justified by domain knowledge, not by convenience or by the data.

Every edge and every non-edge requires a justification. Before touching data, you should be able to answer for each pair of variables: "Why do I believe X does (or does not) directly cause Y?"

The DAG encodes domain knowledge, not statistical patterns. It cannot be learned from data alone. Two datasets with identical joint distributions can correspond to different causal DAGs. Fitting models to data and then reading off a DAG is not causal inference — it is pattern matching.

---

DAG specification in code

import networkx as nx
from dowhy import CausalModel

dag = nx.DiGraph()
dag.add_edges_from([
    ("treatment", "outcome"),
    ("confounder", "treatment"),
    ("confounder", "outcome"),
])

model = CausalModel(
    data=df,
    treatment="treatment",
    outcome="outcome",
    graph=dag,
)

Warning: DoWhy requires explicit U nodes for unobserved confounders. If you omit them, DoWhy assumes there is no unobserved confounding — a very strong assumption that is almost never defensible in observational data. Always add unobserved nodes explicitly:

dag.add_edges_from([
    ("U_unobserved", "treatment"),
    ("U_unobserved", "outcome"),
])

---

Common causal structures

1. Confounder

C → T
C → Y
T → Y

C affects both treatment and outcome. Without adjusting for C, the T→Y estimate is biased. Adjust for C using regression, matching, IPW, or any valid backdoor adjustment method.

Example: age confounds the relationship between exercise (T) and cardiovascular health (Y) — older people exercise less and have worse health independently.

2. Mediator

T → M → Y
T → Y  (possibly)

M lies on a causal path from T to Y. Do NOT adjust for M if you want the total effect of T on Y. Adjusting for a mediator blocks the indirect path and gives you only the direct effect — which may not be what you want.

Example: a job training program (T) improves income (Y) partly through improved employment (M). Conditioning on employment removes the mechanism you care about.

3. Collider

T → C ← Y

C is caused by both T and Y. NEVER adjust for C. Conditioning on a collider opens a non-causal path between T and Y, inducing spurious association even when T and Y are independent. See Collider bias warning.

4. Instrument

Z → T → Y
Z ⊥ Y | T  (Z affects Y only through T)

Z is a valid instrument when it (1) affects T (relevance), (2) is independent of all confounders of T→Y (exogeneity), and (3) affects Y only via T (exclusion restriction). Use IV estimation — see references/quasi-experiments.md.

---

Identification strategies

Backdoor criterion

A set of variables S satisfies the backdoor criterion if: 1. No variable in S is a descendant of T. 2. S blocks every path between T and Y that has an arrow into T (backdoor paths).

When S satisfies the backdoor criterion, adjusting for S identifies the causal effect of T on Y. This is the most common identification strategy in observational studies.

Frontdoor criterion

When all paths from T to Y go through a mediator M, and M has no unobserved confounders with Y, the frontdoor criterion applies. Rare in practice, but it allows identification even when T and Y have unobserved common causes. The estimator chains two regressions: T→M and M→Y, adjusting for T in the second.

Instrumental variables

When a valid instrument Z exists (relevance + exogeneity + exclusion), use two-stage least squares (2SLS) or a Bayesian IV model. IV gives a Local Average Treatment Effect (LATE) — the effect for compliers only, not the full population. Be explicit about this restriction when reporting results.

See references/quasi-experiments.md for implementation details.

Design-based identification

Randomized assignment, regression discontinuity (RDD), difference-in-differences (DiD), and interrupted time series (ITS) provide identification through study design rather than adjustment. These require their own assumptions (parallel trends for DiD, continuity for RDD) which must be stated and tested where possible.

See references/quasi-experiments.md.

---

Finding adjustment sets

identified_estimand = model.identify_effect(
    proceed_when_unidentifiable=False,  # NEVER set to True silently
)
print(identified_estimand)

DoWhy will enumerate valid adjustment sets given your DAG. Review the output carefully — it tells you which variables to include and which identification strategy it found.

Warning: proceed_when_unidentifiable=True instructs DoWhy to continue even when causal identification fails. Never use this option without explicitly warning the user that the resulting estimate is not causally identified and should be interpreted as an association, not an effect. Using it silently is a form of scientific misconduct.

When DoWhy reports the effect is not identifiable, the correct response is to revise the DAG with the user (add instruments, rethink confounders, consider a design-based approach) — not to override the check.

---

Collider bias warning

Adjusting for a collider opens a non-causal path between treatment and outcome, creating spurious associations that do not exist in the population. This is one of the most common and most damaging mistakes in causal inference.

Classic example: you study the effect of a disease (T) on recovery (Y). You condition on hospitalization (H) because your data comes from a hospital. But both disease severity (a confounder) and recovery both affect who gets hospitalized — making H a collider. Conditioning on H induces a spurious negative association between disease and recovery even if none exists.

Severity → H ← Recovery
T (disease) → H
T → Y (recovery)

This mistake is especially easy to make when:

• You filter your sample on a post-treatment variable.
• You include a variable "to control for sample selection."
• You condition on any variable that is a common effect of T and Y (or their causes).

Always trace paths in your DAG before adding any variable to an adjustment set.

---

Asking the user: prompt templates

Before proceeding with any causal analysis, confirm the DAG and its assumptions with the user. Use these templates verbatim or adapt them.

DAG confirmation prompt

"Here is the causal graph I'm proposing:

>

Nodes: [list all variables]

Edges (direct causal effects): [list all arrows, e.g., 'age → treatment', 'treatment → outcome']

Non-edges (explicit no-direct-effect assumptions): [list key omitted edges, e.g., 'income does NOT directly affect outcome, only through treatment']

>

I'm assuming [X] does NOT directly cause [Y] — is that correct? Please confirm or correct each assumption before I proceed. Changing even one edge can change which variables to adjust for and whether the effect is identified at all."

Assumption confirmation prompt

"This analysis rests on the following untestable assumptions:

>

1. [Assumption 1, e.g., 'No unobserved confounders of the treatment–outcome relationship other than those listed.']

2. [Assumption 2, e.g., 'The instrument Z affects Y only through T (exclusion restriction).']

3. [Assumption 3, e.g., 'The DAG is acyclic — no feedback loops exist between any variables.']

>

Are you comfortable defending these in a peer review or stakeholder context? If any feel fragile, I will flag them prominently in the report and, where possible, run a sensitivity analysis to quantify how much hidden confounding would be needed to overturn the conclusion."

Quasi-Experimental Designs

Contents

• Design selection guide
• Difference-in-Differences (DiD)
• Staggered DiD
• Synthetic Control
• Interrupted Time Series (ITS)
• Piecewise ITS
• Regression Discontinuity (RDD)
• Regression Kink Design
• Instrumental Variables (IV)
• Inverse Propensity Score Weighting (IPSW)
• Accessing results (all designs)

---

Design selection guide

If none of these fit, see references/structural-models.md for the pm.do() / pm.observe() path.

Always prefer PyMC-backed models over sklearn. PyMC models give full posterior uncertainty; sklearn gives point estimates only. Every template below uses cp.pymc_models.*.

---

Difference-in-Differences (DiD)

Use when: Treatment switches on at a known calendar time, and an untreated control group is available.

Key assumption: Parallel trends — treatment and control would have followed the same trajectory absent treatment. Partially testable (check pre-treatment trends) but never fully verifiable post-treatment.

CausalPy class: cp.DifferenceInDifferences

import causalpy as cp
import numpy as np

rng = np.random.default_rng(sum(map(ord, "did-analysis")))

result = cp.DifferenceInDifferences(
    data=df,
    formula="outcome ~ 1 + post_treatment + group + post_treatment:group",
    time_variable_name="time",
    group_variable_name="group",
    # group must be dummy-coded: 1 = treatment, 0 = control (NOT string labels)
    # post_treatment must also be 0/1 (defaults to column named "post_treatment")
    # data must also have a "unit" column labeling unique units (used for plotting)
    model=cp.pymc_models.LinearRegression(
        sample_kwargs={"nuts_sampler": "nutpie", "random_seed": rng}
    ),
)
fig, ax = result.plot()
es = result.effect_summary(direction="two-sided")
print(es.text)
print(es.table)

The group column must be dummy-coded (0/1 integers), not string labels — CausalPy will reject string-valued group variables. The data must also have a unit column identifying individual units (used for plotting). The interaction post_treatment:group is the DiD estimator.

What can go wrong: Parallel trends violated (different pre-treatment trajectories); compositional changes in group membership over time; anticipation effects (units respond before the official treatment date); SUTVA violations (treated units affect controls).

Refutation: See references/refutation.md — placebo treatment time, parallel trends test (visual + formal).

---

Staggered DiD

Use when: Treatment rolls out at different times for different units. Standard DiD with a single date would be misspecified.

Key assumption: Parallel trends for each cohort — each treated cohort's counterfactual is well-approximated by not-yet-treated and never-treated units.

CausalPy class: cp.StaggeredDifferenceInDifferences

import causalpy as cp
import numpy as np

rng = np.random.default_rng(sum(map(ord, "staggered-did")))

result = cp.StaggeredDifferenceInDifferences(
    data=df,
    formula="outcome ~ 1 + C(unit)",
    unit_variable_name="unit",
    time_variable_name="time",
    treated_variable_name="treated",
    treatment_time_variable_name="treatment_date",
    # never_treated_value=np.inf    # default; change if data uses -1 or 0 as sentinel
    # event_window=(-5, 10)         # optional: event-study style (periods relative to treatment)
    model=cp.pymc_models.LinearRegression(
        sample_kwargs={"nuts_sampler": "nutpie", "random_seed": rng}
    ),
)
fig, ax = result.plot()
es = result.effect_summary(direction="two-sided")
print(es.text)

What can go wrong: TWFE DiD uses already-treated units as controls for later-treated units — biased when treatment effects evolve. CausalPy's staggered implementation mitigates this; verify the comparisons being made. Few cohorts (< 5) produce noisy cohort-specific estimates.

Refutation: See references/refutation.md — Bacon decomposition, cohort-specific parallel trends checks.

---

Synthetic Control

Use when: You have a single treated unit and a pool of untreated donors whose pre-treatment history can be weighted to approximate the counterfactual.

Key assumption: The treated unit's pre-treatment outcomes can be reproduced as a weighted combination of donors. The treated unit must lie within (or near) the convex hull of donors. CausalPy 0.8+ warns about convex hull violations.

CausalPy class: cp.SyntheticControl

import causalpy as cp
import pandas as pd
import numpy as np

rng = np.random.default_rng(sum(map(ord, "synthetic-control")))

result = cp.SyntheticControl(
    data=df_wide,  # MUST be wide format: index=time, columns=unit names, values=outcome
    treatment_time=pd.Timestamp("2020-01-01"),
    control_units=["unit_A", "unit_B", "unit_C"],
    treated_units=["unit_treated"],
    # SC does NOT use a formula — it finds optimal weights over donors directly
    # If your data is long format, pivot first:
    # df_wide = df.pivot(index="date", columns="unit", values="outcome")
    model=cp.pymc_models.WeightedSumFitter(
        sample_kwargs={"nuts_sampler": "nutpie", "random_seed": rng}
    ),
)
fig, ax = result.plot()
es = result.effect_summary(
    direction="two-sided",
    cumulative=True,   # total cumulative effect over post-treatment window
    relative=True,     # effect as % of counterfactual
)
print(es.text)
print(es.table)

What can go wrong: Convex hull violation (treated unit is an outlier — inspect pre-treatment fit RMSE); donor contamination (donors affected by spillovers — remove them); fewer than ~5 donors (poorly constrained weights); poor pre-treatment fit means the counterfactual is unreliable.

Refutation: See references/refutation.md — leave-one-out donors, pre-treatment RMSE, placebo treatments on each control unit.

---

Interrupted Time Series (ITS)

Use when: Single time series, no control group, intervention at a known time. The pre-intervention period defines the counterfactual trend. Weaker than DiD — relies on the pre-trend extrapolating forward with no contaminating events.

Key assumption: No confounding event at the treatment time. Requires an explicit user checkpoint.

CausalPy class: cp.InterruptedTimeSeries

import causalpy as cp
import pandas as pd
import numpy as np

rng = np.random.default_rng(sum(map(ord, "its-analysis")))

result = cp.InterruptedTimeSeries(
    data=df,
    treatment_time=pd.Timestamp("2020-01-01"),
    formula="y ~ 1 + t + C(month)",   # t = numeric time index; C(month) absorbs seasonality
    model=cp.pymc_models.LinearRegression(
        sample_kwargs={"nuts_sampler": "nutpie", "random_seed": rng}
    ),
)
fig, ax = result.plot()
es = result.effect_summary(
    direction="two-sided",
    cumulative=True,
    window="post",   # "post" = all post-treatment points (default); or pass a tuple/slice
)
print(es.text)
print(es.table)

⚠️ MANDATORY USER CHECKPOINT: Before fitting, ask: "Did anything else change at [treatment_time] that could explain a break in the series?" If yes, the design may not identify the causal effect.

What can go wrong: Confounding event (the most common and most serious ITS failure — no statistical test can catch it); autocorrelation in residuals; non-linear pre-trend (add polynomial time terms if needed); too few pre-treatment observations (aim for 12–24+ periods).

Refutation: See references/refutation.md — placebo treatment time, autocorrelation diagnostics, confounding event search (user checkpoint).

---

Piecewise ITS

Use when: Multiple interventions occurred at different times, or the series has known structural breaks creating distinct segments.

Key assumption: Same as ITS — no confounding events at each breakpoint. Assumptions multiply with the number of breaks.

CausalPy class: cp.PiecewiseITS (requires CausalPy 0.8+)

import causalpy as cp
import numpy as np

rng = np.random.default_rng(sum(map(ord, "piecewise-its")))

# step() = level shift (0 before date, 1 after)
# ramp() = slope change (0 before date, increases linearly after)
result = cp.PiecewiseITS(
    data=df,
    formula=(
        "y ~ 1 + t"
        " + step(t, change_date_1) + ramp(t, change_date_1)"
        " + step(t, change_date_2)"   # add ramp() if slope also changes at date 2
    ),
    model=cp.pymc_models.LinearRegression(
        sample_kwargs={"nuts_sampler": "nutpie", "random_seed": rng}
    ),
)
fig, ax = result.plot()

step() and ramp() are CausalPy formula transforms exported as cp.step and cp.ramp.

Fallback if PiecewiseITS is unavailable: Use pymc_extras.statespace (https://github.com/pymc-devs/pymc-extras/tree/main/pymc_extras/statespace) — PyMC's state-space submodule. Fall back to a manual PyMC model only if that doesn't fit. Document whichever path is taken.

What can go wrong: All ITS failures, multiplied by the number of breaks. Overfitting with too many segments; break point uncertainty (if timing is approximate, fixing exact dates introduces misspecification).

Refutation: Placebo break points, autocorrelation checks, confounding event search at each break.

---

Regression Discontinuity (RDD)

Use when: Treatment is assigned by whether a running variable (score, age, index) crosses a threshold. Units just below and just above the threshold are similar on all other characteristics — the threshold creates local random assignment.

Key assumption: No manipulation at the threshold — units cannot sort themselves precisely to one side. Requires a smooth density of the running variable through the threshold (McCrary density test).

CausalPy class: cp.RegressionDiscontinuity

import causalpy as cp
import numpy as np

rng = np.random.default_rng(sum(map(ord, "rdd-analysis")))

result = cp.RegressionDiscontinuity(
    data=df,
    formula="outcome ~ 1 + running_var + treated + running_var:treated",
    running_variable_name="running_var",
    treatment_threshold=cutoff_value,
    model=cp.pymc_models.LinearRegression(
        sample_kwargs={"nuts_sampler": "nutpie", "random_seed": rng}
    ),
    bandwidth=bandwidth_value,   # restrict to obs within ±bandwidth of threshold
)
fig, ax = result.plot()
es = result.effect_summary(direction="two-sided")
print(es.text)
print(es.table)

running_var:treated lets the slope differ on each side. For non-linear relationships, add polynomial terms: "outcome ~ 1 + running_var + I(running_var**2) + treated + running_var:treated". Start bandwidth selection with a data-driven approach (Imbens-Kalyanaraman), then check sensitivity in refutation.

What can go wrong: Bunching at threshold (manipulation — caught by McCrary test); misspecified functional form (linear fit attributes non-linearity to a spurious discontinuity); bandwidth too wide (weakens local randomization argument); bandwidth too narrow (noisy estimates); covariates jumping at the threshold (selection on observables).

Refutation: See references/refutation.md — McCrary density test, bandwidth sensitivity, covariate balance at threshold, placebo thresholds.

---

Regression Kink Design

Use when: Treatment intensity (not binary assignment) changes slope at a known threshold. You look for a change in the slope of the outcome, not a level jump.

Key assumption: No manipulation at the kink point, smooth density of running variable. The relationship between the running variable and outcome must be smooth everywhere except at the kink.

CausalPy class: cp.RegressionKink

import causalpy as cp
import numpy as np

rng = np.random.default_rng(sum(map(ord, "rk-analysis")))

result = cp.RegressionKink(
    data=df,
    formula="outcome ~ 1 + running_var + treated + running_var:treated",
    running_variable_name="running_var",
    kink_point=kink_value,
    model=cp.pymc_models.LinearRegression(
        sample_kwargs={"nuts_sampler": "nutpie", "random_seed": rng}
    ),
    bandwidth=bandwidth_value,
)
fig, ax = result.plot()
es = result.effect_summary(direction="two-sided")
print(es.text)

RDD vs. Regression Kink: RDD detects a level jump (treatment switches 0→1). Regression Kink detects a slope change (treatment intensity changes continuously). Example: a benefit program that phases out linearly above an income threshold — the kink is where the phase-out begins.

What can go wrong: Same as RDD, plus a pre-existing kink in the outcome coinciding with the kink point; small kink angle makes the design underpowered.

Refutation: Same as RDD — bandwidth sensitivity, covariate balance, placebo kink points.

---

Instrumental Variables (IV)

Use when: Treatment is endogenous (confounded), but you have an instrument that (a) affects treatment and (b) affects the outcome only through treatment (exclusion restriction). Classic instruments: lottery assignment, distance to facility, policy eligibility cutoffs.

Key assumptions: 1. Relevance: Instrument strongly predicts treatment. Weak instruments → severely biased estimates. 2. Exclusion restriction: Instrument affects outcome only through treatment. Untestable — must be defended on substantive grounds. 3. Independence: Instrument is as-good-as-randomly assigned.

CausalPy class: cp.InstrumentalVariable

import causalpy as cp
import numpy as np

rng = np.random.default_rng(sum(map(ord, "iv-analysis")))

result = cp.InstrumentalVariable(
    instruments_data=df,
    data=df,
    instruments_formula="treatment ~ 1 + instrument",   # first stage
    formula="outcome ~ 1 + treatment",                  # second stage
    model=cp.pymc_models.InstrumentalVariableRegression(
        sample_kwargs={"nuts_sampler": "nutpie", "random_seed": rng}
    ),
)
fig, ax = result.plot()

IV estimates the Local Average Treatment Effect (LATE) — the effect for compliers (units whose treatment status is changed by the instrument). Report this distinction; LATE ≠ ATE unless all units comply.

What can go wrong: Weak instruments (first-stage F < 10 is a warning; F < 5 is very weak — IV estimates are then more biased than OLS); exclusion restriction violated (instrument has a direct effect on outcome — untestable); LATE ≠ ATE when effect heterogeneity is large.

Note: Design-specific refutation for IV (weak instrument tests, overidentification) is deferred to v1.1. Use general DoWhy refuters for v1.0. See references/refutation.md.

---

Inverse Propensity Score Weighting (IPSW)

Use when: Observational data, no natural experiment. Reweight the sample so treatment and control groups have similar covariate distributions, by modeling the probability of treatment given observed covariates (propensity score).

Key assumptions: 1. No unmeasured confounders (strong ignorability): All variables that affect both treatment and outcome are observed. Untestable. 2. Positivity: Every unit has non-zero probability of either treatment status. Extreme propensities → extreme weights → unreliable estimates.

CausalPy class: cp.InversePropensityWeighting

import causalpy as cp
import numpy as np

rng = np.random.default_rng(sum(map(ord, "ipsw-analysis")))

result = cp.InversePropensityWeighting(
    data=df,
    formula="treatment ~ 1 + confounder1 + confounder2 + confounder3",
    outcome_variable="outcome",
    weighting_scheme="robust",   # "raw": 1/p weights — extreme weights possible
                                 # "robust": Hajek normalized weights — recommended default
                                 # "doubly robust": IPW + outcome model — consistent if either is correct
                                 # "overlap": trims extreme propensities — reduces variance
    model=cp.pymc_models.PropensityScore(
        sample_kwargs={"nuts_sampler": "nutpie", "random_seed": rng}
    ),
)
fig, ax = result.plot()
es = result.effect_summary(direction="two-sided")
print(es.text)
print(es.table)

Start with "robust". Use "doubly robust" when you also have a plausible outcome model. Use "overlap" when propensity scores cluster near 0 or 1.

What can go wrong: Unmeasured confounders (the most critical threat — propensity methods only adjust for observed covariates); positivity violations (extreme weights dominate the estimate — inspect propensity distribution); propensity model misspecification (check covariate balance after weighting); limited overlap (IPSW extrapolates if groups don't share covariate support).

Note: Design-specific refutation for IPSW (covariate balance after weighting, extreme weight diagnostics) is deferred to v1.1. Use general DoWhy refuters for v1.0. See references/refutation.md.

---

Accessing results (all designs)

All CausalPy result objects share a common interface.

# Visualization — always available
fig, ax = result.plot()

# Effect summary — parameters vary by design:
#   DiD, RDD, IV, IPSW: direction, alpha, min_effect
#   SC and ITS additionally: cumulative, relative, window
es = result.effect_summary(direction="two-sided")
print(es.text)    # prose narrative
print(es.table)   # structured DataFrame summary

# SC / ITS with cumulative effects:
es = result.effect_summary(direction="two-sided", cumulative=True, relative=True)

# Underlying InferenceData for deeper analysis:
idata = result.idata
# Delegate all diagnostics to bayesian-workflow:
# import arviz_stats as azs; azs.diagnose(idata)

# Log-likelihood is NOT stored automatically — compute explicitly if needed for LOO:
# with result.model:
#     pm.compute_log_likelihood(idata)

Always check the class docstring (help(cp.DifferenceInDifferences)) for the full list of effect_summary parameters — they differ across designs.

Refutation and Sensitivity Analysis

Refutation is MANDATORY for every causal analysis. No exceptions. If you skip this step, the analysis cannot support causal claims.

Contents

1. Refutation principles 2. DiD refutation recipes 3. Synthetic Control refutation recipes 4. RDD refutation recipes 5. ITS refutation recipes 6. General refutation (all designs — DoWhy) 7. Sensitivity to unobserved confounding 8. Pass/fail interpretation 9. What to do when refutation fails

---

Refutation principles

Every causal claim rests on assumptions that cannot be directly tested with the available data. Refutation does not test whether the assumptions are true — it tests whether your conclusions are robust to plausible violations of those assumptions.

The logic of refutation:

• Pass does NOT prove causality. It means you have tried to falsify your result and have not succeeded. The result remains plausible, not proven.
• Fail means the causal claim is suspect. The data are consistent with a world where your assumptions are violated and the apparent effect is spurious.
• The value of refutation is asymmetric: a pass gives limited reassurance, a fail gives strong warning.

Mandatory refutation checklist — run every test for your design, report all results:

Report results as a pass/fail table. Do not bury failures in footnotes.

---

DiD refutation recipes

Placebo treatment time

Shift the treatment to a point in the pre-treatment period. If your parallel trends assumption holds and the real effect is real, the "effect" estimated at the placebo time should be indistinguishable from zero.

import causalpy as cp
import numpy as np

# Use only pre-treatment data
df_pre_only = df[df["time"] < treatment_time].copy()
placebo_time = df_pre_only["time"].median()  # midpoint of pre-period

result_placebo = cp.DifferenceInDifferences(
    data=df_pre_only,
    formula=same_formula,
    time_variable_name="time",
    group_variable_name="group",
    model=cp.pymc_models.LinearRegression(
        sample_kwargs={"random_seed": rng, "nuts_sampler": "nutpie"}
    ),
)

# PASS if: placebo effect HDI includes zero
print(result_placebo.summary())

Interpretation: If the placebo effect is large and its HDI excludes zero, the parallel trends assumption is likely violated — units were already diverging before treatment. The real effect estimate is unreliable.

Parallel trends test

Before any modeling, plot pre-treatment outcomes for treatment and control groups. They should track each other closely. Divergence in the pre-period is the strongest visual signal that DiD is invalid.

import matplotlib.pyplot as plt

pre = df[df["time"] < treatment_time].copy()

fig, ax = plt.subplots(figsize=(9, 4))
for g in pre["group"].unique():
    gd = pre[pre["group"] == g]
    ax.plot(gd["time"], gd["outcome"], label=str(g), marker="o")

ax.axvline(treatment_time, color="red", linestyle="--", alpha=0.6, label="Treatment")
ax.set_title("Pre-treatment trends (should be parallel)")
ax.set_xlabel("Time")
ax.set_ylabel("Outcome")
ax.legend()
plt.tight_layout()

What to look for: Trend lines should be approximately parallel — same slope, even if different levels. If they are converging, diverging, or crossing, parallel trends is violated.

If the trends are not parallel: DiD is not valid for this data. Consider a matching pre-trend model, or switch designs (synthetic control handles heterogeneous pre-trends better).

Bacon decomposition (staggered DiD)

When treatment is staggered (units adopt treatment at different times), the standard DiD estimator is a weighted average of many 2×2 DiD sub-estimates — some of which use already-treated units as controls. This can produce negatively-weighted comparisons that flip the sign of the true effect.

# Use the bacon_decomposition package or inspect manually
# Each 2x2 estimate: early-treated vs control, late-treated vs control,
# early-treated vs late-treated (problematic if effects are heterogeneous)

# Flag: if any 2x2 comparison uses a treated unit as control, report it
# Flag: if weights are negative, the overall estimate may be uninformative

Interpretation: Large variation across sub-estimates, or negative weights on sub-estimates, indicates treatment effect heterogeneity that the standard TWFE estimator cannot handle. Use a staggered DiD estimator (Callaway-Sant'Anna, Sun-Abraham) instead.

---

Synthetic Control refutation recipes

Leave-one-out donors

Re-fit the synthetic control model removing each donor unit one at a time. If the post-treatment trajectory changes dramatically when any single donor is removed, the result is fragile and dependent on that one unit.

donors = df[df["group"] == "donor"]["unit"].unique()

loo_estimates = {}
for drop_unit in donors:
    df_loo = df[df["unit"] != drop_unit].copy()
    control_units_loo = [u for u in control_units if u != drop_unit]
    sc_loo = cp.SyntheticControl(
        data=df_loo,
        treatment_time=treatment_time,
        control_units=control_units_loo,
        treated_units=treated_units,
        model=cp.pymc_models.WeightedSumFitter(
            sample_kwargs={"random_seed": rng, "nuts_sampler": "nutpie"}
        ),
    )
    loo_estimates[drop_unit] = sc_loo.effect_summary()

# PASS if: post-treatment effect estimates are stable across all LOO runs
# FAIL if: removing any single donor shifts the estimate substantially

Pre-treatment fit quality

The synthetic control must closely track the treated unit in the pre-treatment period. Poor pre-treatment fit means the synthetic control is not a valid counterfactual, and the post-treatment gap is not interpretable as a causal effect.

import numpy as np

pre_actual = df[(df["group"] == "treated") & (df["time"] < treatment_time)]["outcome"].values
# Pre-period synthetic values are stored in sc_result.datapre (an xarray Dataset)
pre_synthetic = sc_result.datapre["synthetic"].values

rmse = np.sqrt(np.mean((pre_actual - pre_synthetic) ** 2))
print(f"Pre-treatment RMSE: {rmse:.4f}")

# PASS if: RMSE is small relative to the scale of the outcome and the post-treatment gap
# Rule of thumb: RMSE < 10% of the post-treatment effect estimate warrants confidence

If pre-treatment fit is poor: Do not interpret the post-treatment gap causally. Add more predictors to the formula, expand the donor pool, or switch to a Bayesian structural time series model.

Placebo treatments on controls

Apply the synthetic control procedure to each control unit as if it were the treated unit, using the remaining controls as donors. Compute the post-treatment "effect" for each placebo. If many placebo units show effects as large as the real treated unit, the result is not statistically meaningful.

placebo_effects = {}
for placebo_treated in donors:
    placebo_donors = [u for u in donors if u != placebo_treated]
    sc_placebo = cp.SyntheticControl(
        data=df[df["unit"] != treated_unit].copy(),
        treatment_time=treatment_time,
        control_units=placebo_donors,
        treated_units=[placebo_treated],
        model=cp.pymc_models.WeightedSumFitter(
            sample_kwargs={"random_seed": rng, "nuts_sampler": "nutpie"}
        ),
    )
    summary = sc_placebo.effect_summary()
    placebo_effects[placebo_treated] = summary["mean"].values[0]

# Compute p-value analogue: fraction of placebos with effect >= real effect
real_summary = sc_result.effect_summary()
real_effect = real_summary["mean"].values[0]
p_value_analogue = np.mean([abs(v) >= abs(real_effect) for v in placebo_effects.values()])
print(f"Fraction of placebos as extreme as treated: {p_value_analogue:.2f}")

# PASS if: p_value_analogue < 0.10 (fewer than 1 in 10 placebos match the real effect)

---

RDD refutation recipes

McCrary density test

Test for bunching at the threshold in the running variable. If units can manipulate whether they fall just above or just below the threshold, the "as-good-as-random" assumption of RDD is violated. A sharp spike in density just below (for a threshold where being above means treatment) is the classic manipulation signature.

# Visual inspection
import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax.hist(df["running_var"], bins=50, edgecolor="white")
ax.axvline(threshold, color="red", linestyle="--", label="Threshold")
ax.set_title("Running variable density (check for bunching at threshold)")
ax.legend()

# Formal test: use rddensity (R) or a Python port
# If density shows a discontinuity at the threshold, RDD is invalid

If bunching is detected: Report that the RDD design assumption is violated. Units near the threshold are not comparable. Do not proceed with causal claims. Consider restricting to units far from the threshold (fuzzy analysis) or switching designs.

Bandwidth sensitivity

The RDD estimate should be robust to reasonable choices of bandwidth. If the estimate changes dramatically as bandwidth varies, the result is sensitive to an arbitrary analyst choice.

import causalpy as cp

base_bw = rdd_result.bandwidth  # whatever the default or IK-optimal bandwidth was
bandwidths = [base_bw * f for f in [0.5, 0.75, 1.0, 1.25, 1.5, 2.0]]

bw_estimates = []
for bw_val in bandwidths:
    r = cp.RegressionDiscontinuity(
        data=df,
        formula=formula,
        running_variable_name="running_var",
        treatment_threshold=threshold,
        model=cp.pymc_models.LinearRegression(
            sample_kwargs={"random_seed": rng, "nuts_sampler": "nutpie"}
        ),
        bandwidth=bw_val,
    )
    summary = r.summary()
    bw_estimates.append({"bandwidth": bw_val, "estimate": summary["mean"], "hdi_lo": summary["hdi_3%"], "hdi_hi": summary["hdi_97%"]})

# Plot bandwidth sensitivity
import pandas as pd
bw_df = pd.DataFrame(bw_estimates)
fig, ax = plt.subplots()
ax.plot(bw_df["bandwidth"], bw_df["estimate"], marker="o")
ax.fill_between(bw_df["bandwidth"], bw_df["hdi_lo"], bw_df["hdi_hi"], alpha=0.2)
ax.axhline(0, color="gray", linestyle="--")
ax.set_xlabel("Bandwidth")
ax.set_ylabel("Effect estimate")
ax.set_title("Bandwidth sensitivity (should be stable)")

# PASS if: estimates are stable and HDIs overlap across bandwidths

Covariate balance at threshold

Covariates (pre-determined characteristics) should not jump discontinuously at the threshold. A covariate jump indicates that units on either side of the threshold differ in ways unrelated to treatment — violating the continuity assumption.

# Run a separate RDD with each covariate as the outcome
covariates = ["age", "income", "prior_outcome"]  # replace with your pre-treatment covariates

for cov in covariates:
    r_cov = cp.RegressionDiscontinuity(
        data=df,
        formula=f"{cov} ~ 1 + running_var",
        running_variable_name="running_var",
        treatment_threshold=threshold,
        model=cp.pymc_models.LinearRegression(
            sample_kwargs={"random_seed": rng, "nuts_sampler": "nutpie"}
        ),
        bandwidth=base_bw,
    )
    print(f"Covariate balance — {cov}:")
    print(r_cov.summary())
    # PASS if: effect HDI includes zero for all covariates

Placebo thresholds

Move the threshold to other values in the pre-treatment distribution. No effect should appear at values that do not correspond to the true treatment threshold.

running_var_quantiles = df["running_var"].quantile([0.25, 0.40, 0.60, 0.75]).tolist()

for placebo_thresh in running_var_quantiles:
    if abs(placebo_thresh - threshold) < base_bw:
        continue  # skip placebo thresholds too close to the real one
    r_placebo = cp.RegressionDiscontinuity(
        data=df[df["running_var"] < threshold],  # use only pre-threshold data
        formula=formula,
        running_variable_name="running_var",
        treatment_threshold=placebo_thresh,
        model=cp.pymc_models.LinearRegression(
            sample_kwargs={"random_seed": rng, "nuts_sampler": "nutpie"}
        ),
        bandwidth=base_bw,
    )
    print(f"Placebo threshold {placebo_thresh:.2f}: {r_placebo.summary()}")
    # PASS if: effect HDI includes zero at all placebo thresholds

---

ITS refutation recipes

Placebo treatment time

Same logic as the DiD placebo. Shift treatment to the middle of the pre-period and estimate the "effect." If real, it should be near zero.

df_pre_only = df[df["time"] < treatment_time].copy()
placebo_time = df_pre_only["time"].median()

result_placebo = cp.InterruptedTimeSeries(
    data=df_pre_only,
    treatment_time=placebo_time,
    formula=formula,
    model=cp.pymc_models.LinearRegression(
        sample_kwargs={"random_seed": rng, "nuts_sampler": "nutpie"}
    ),
)

# PASS if: placebo effect HDI includes zero
print(result_placebo.effect_summary())

Autocorrelation diagnostics

ITS fits a regression to time-series data. Residuals from time-series regressions are often autocorrelated. If they are, the effective sample size is smaller than the nominal sample size, and standard errors are underestimated — potentially producing a spurious significant result.

import statsmodels.stats.stattools as sms

# Extract posterior mean residuals
residuals = its_result.get_residuals()  # or compute manually

# Durbin-Watson statistic: 2 = no autocorrelation, < 2 = positive, > 2 = negative
dw = sms.durbin_watson(residuals)
print(f"Durbin-Watson: {dw:.3f}")
# Values between 1.5 and 2.5 are generally acceptable

# Plot ACF
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
plot_acf(residuals, lags=20)
plt.title("Residual autocorrelation (spikes outside band = problem)")
plt.tight_layout()

# If autocorrelation is present: add AR terms to the model formula,
# or use a model with explicit autocorrelation structure

Confounding event search

ITS is especially vulnerable to events that happened at the same time as the treatment but are not the treatment. There is no statistical test for this — it requires domain knowledge.

Ask the user: "Did anything else change at or near the treatment time that could explain the observed trend change? Examples: a policy change elsewhere, a market shock, a shift in data collection practices, a simultaneous intervention on a related variable."

If the user identifies a potential confound, either (a) control for it by adding it to the model, (b) collect data on a control unit that experienced the confound but not the treatment (converting ITS to DiD), or (c) downgrade the causal claim.

---

General refutation (all designs — DoWhy)

These three refuters apply to any analysis where you have a DAG and a DoWhy model. Run all three. They probe complementary failure modes.

import dowhy

dowhy_model = dowhy.CausalModel(
    data=df,
    treatment="treatment",
    outcome="outcome",
    graph=dag,  # your nx.DiGraph from dags-and-identification.md
)

identified = dowhy_model.identify_effect(proceed_when_unidentifiable=False)
estimate = dowhy_model.estimate_effect(
    identified,
    method_name="backdoor.linear_regression",
)

# 1. Random common cause
# Adds a random covariate as an unobserved common cause.
# If the estimate changes substantially, the model is fragile.
ref_rcc = dowhy_model.refute_estimate(
    identified, estimate, method_name="random_common_cause"
)
print(ref_rcc)
# PASS if: new estimate ≈ original estimate (within sampling noise)

# 2. Placebo treatment
# Replaces the real treatment with a random permutation.
# The effect should vanish — a real causal effect cannot survive random treatment assignment.
ref_placebo = dowhy_model.refute_estimate(
    identified, estimate, method_name="placebo_treatment_refuter"
)
print(ref_placebo)
# PASS if: refuted estimate ≈ 0

# 3. Data subset
# Re-estimates on random 80% subsets of the data.
# A robust effect should be stable across subsets.
ref_subset = dowhy_model.refute_estimate(
    identified, estimate, method_name="data_subset_refuter", subset_fraction=0.8
)
print(ref_subset)
# PASS if: estimate stable (HDIs overlap across subsets)

---

Sensitivity to unobserved confounding

This is the most important sensitivity analysis for observational data. Even if all other refutations pass, an unobserved confounder could explain your result. The question is: how strong would that confounder need to be?

A weak effect that can be explained by a modest confounder is fragile. A large effect that requires an implausibly strong confounder to be explained away is credible.

# Start with a small effect strength and increase until the estimate is "explained away"
for strength in [0.05, 0.1, 0.2, 0.3, 0.5]:
    ref_confound = dowhy_model.refute_estimate(
        identified,
        estimate,
        method_name="add_unobserved_common_cause",
        confounders_effect_on_treatment="binary_flip",
        confounders_effect_on_outcome="linear",
        effect_strength_on_treatment=strength,
        effect_strength_on_outcome=strength,
    )
    print(f"Confounder strength {strength}: refuted estimate = {ref_confound.new_effect:.4f}")

# The "tipping point" is the strength at which the effect is driven to zero
# Report: "The effect is robust to unobserved confounders up to strength X.
#          A confounder of strength > X would be required to explain away the result.
#          For context, [observed covariate Z] has an association strength of ~Y."

Interpreting the tipping point: Compare the tipping-point strength to the observed associations of your measured confounders. If explaining away the result requires a confounder much stronger than anything you have measured, the result is credible. If it requires only a modest confounder, report with appropriate humility.

Reporting language template:

"The estimated effect of [treatment] on [outcome] is [X] (95% HDI: [lo, hi]). Sensitivity analysis shows this estimate is robust to unobserved confounders with effect strengths up to approximately [tipping point]. Given that the strongest observed confounder ([variable]) has an association of [observed strength], an unobserved confounder capable of explaining away this result would need to be approximately [ratio]× stronger than any variable we have measured."

---

Pass/fail interpretation

Critical tests (failure of any one is disqualifying):

• Placebo treatment time shows large non-zero effect
• McCrary density test shows bunching at RDD threshold
• Pre-treatment RMSE (SC) is comparable in size to the post-treatment gap
• DoWhy placebo treatment does not drive effect to zero
• Unobserved confounder tipping point is implausibly low

Marginal tests (failure warrants caveats but not disqualification):

• Bandwidth sensitivity shows modest variation within overlapping HDIs
• Durbin-Watson indicates mild autocorrelation in ITS residuals
• Leave-one-out SC shows one influential donor
• Data subset refuter shows estimate drifts but remains on same side of zero

---

What to do when refutation fails

Rule: do not hide failures. A buried failure in an appendix is still a failure. Report it prominently, at the top of the findings section.

Step-by-step failure response

1. Identify which assumption failed. Is it a testable design assumption (parallel trends, density continuity) or an untestable structural assumption (no unobserved confounding)? Testable failures are more serious — they are direct evidence of a problem, not just a risk.

2. Determine if it is fixable. Some failures suggest a different model or design. Others are fundamental to the data-generating process.

3. Downgrade the language. Use "associated with" instead of "causes." Use "suggestive" instead of "demonstrates." Never use causal language when the design has failed its own refutation tests.

4. Consider alternative designs. If observational adjustment fails, is there a quasi-experimental design available? Can you find an instrument, a discontinuity, or a comparison group that restores identification?

5. Mark the assumption as fragile in the report. Use explicit flagging:

"Warning: the parallel trends assumption is not supported by pre-treatment data. The DiD estimate below should be interpreted as associational, not causal. We retain the analysis for transparency but recommend collecting additional comparison units before drawing policy conclusions."

6. Do not iterate until you get a pass. Re-running refutation with slightly different parameters until it passes is p-hacking by another name. Run the tests once, with pre-specified parameters, and report honestly.

Reporting Causal Analyses

Contents

1. Causal analysis report template 2. Causal language guardrails 3. Decision-relevant HDIs 4. Audience adaptation 5. Common reporting mistakes

---

Causal analysis report template

Every causal analysis produces a report with this mandatory structure. Adapt sections as needed, but do not drop sections 1, 7, or 8 — they are non-negotiable.

1. Causal question

One sentence: "What is the effect of [treatment] on [outcome] in [population]?"

Write this before touching data. If you cannot write this sentence, you do not yet have a causal question — you have a dataset.

2. DAG and assumptions

Include the causal graph (generated with model.plot() or drawn explicitly) and an assumption transparency table:

Every assumption in the table must be discussed — not just listed. For each fragile assumption, state what evidence, if any, supports it.

3. Identification strategy

"We use [method] to identify the causal effect. This is valid because [justification]."

Be explicit: which identification result applies (backdoor, frontdoor, IV, RDD continuity, DiD parallel trends)? Which variables are adjusted for and why? Reference dags-and-identification.md for identification criteria and quasi-experiments.md for design-based methods.

If the effect is not point-identified, say so. Partial identification — bounding the effect rather than pinning it down — is a legitimate and honest result.

4. Estimation

State the model specification: likelihood, priors, and any structural constraints. Summarize diagnostics (R-hat, ESS, divergences). Defer to bayesian-workflow/references/diagnostics.md for diagnostic standards and thresholds. Do not re-explain those standards here; link to them.

For structural causal models, report the DoWhy estimand and estimation method:

estimate = model.estimate_effect(
    identified_estimand,
    method_name="backdoor.linear_regression",
)

5. Results with uncertainty

Report effect size with full posterior distribution and multiple HDIs. Never report only a point estimate.

Example:

"We estimate the policy increased test scores by 4.2 points (50% HDI: [3.1, 5.3]; 95% HDI: [-0.3, 8.9]). The most likely effect is 3–5 points, but we cannot rule out a null effect at 95% credibility. There is an 87% posterior probability the effect is positive."

Include:

• A posterior density plot or forest plot of the causal effect
• The probability of direction: P(effect > 0) or P(effect < threshold)
• Effect size in domain-relevant units, not just standardized coefficients

6. Refutation results

Run all applicable refutation tests and report every result — including failures. Do not cherry-pick passing tests.

For DoWhy refutations:

refute = model.refute_estimate(
    identified_estimand,
    estimate,
    method_name="placebo_treatment_refuter",
)
print(refute)

If a refutation test fails, downgrade the causal language accordingly — see Causal language guardrails.

7. Limitations and threats to validity

This section is mandatory and must be prominent — not buried in an appendix. Decision-makers must see it.

Rank threats by severity. For each:

• State the assumption that might be violated
• Explain the direction of bias if it is violated
• Quantify if possible (E-value, sensitivity analysis, bounding)
• State what additional data or design would resolve the threat

Example:

"The main threat to this conclusion is unobserved confounding. An unobserved variable would need to be associated with both treatment assignment and outcomes with a relative risk greater than 2.3 to explain away the estimated effect (E-value = 2.3). Given the rich covariate set and the DiD design, we consider this unlikely but cannot rule it out without experimental data."

8. Plain-language conclusion

Close every report with one paragraph in plain language, regardless of audience:

"We estimate [treatment] causes [outcome] to change by [effect] ([HDI]), assuming [key assumptions hold]. There is a [P]% probability the effect is positive. The main threat to this conclusion is [biggest weakness]. If that assumption is violated, the true effect could be [direction and magnitude of bias]."

---

Causal language guardrails

The strength of your language must match the strength of your identification. Using causal language without identification is not just imprecise — it is misleading.

Default to the more conservative language when in doubt. Overclaiming in a causal analysis is a more serious error than underclaiming — it can drive bad decisions.

Never use the word "effect" when the analysis state is Descriptive. "Association," "correlation," and "relationship" are correct.

---

Decision-relevant HDIs

Do not default to a fixed HDI width. Choose widths that map to intuitive probabilities for the decision context.

Report multiple widths when useful, especially when the conclusion changes across them:

"The effect is 4.2 points (50% HDI: [3.1, 5.3]; 95% HDI: [-0.3, 8.9]). The most likely effect is 3–5 points, but we cannot rule out a null at 95% credibility."

Always state why you chose the reported HDI width. "We report the 75% HDI because this decision requires us to act if the effect is positive with 3-in-4 confidence" is better than silently presenting a number.

For causal analyses specifically, also report the probability of direction:

# Probability effect is positive
p_positive = (idata.posterior["causal_effect"] > 0).mean().item()
print(f"P(effect > 0) = {p_positive:.2f}")

This is more interpretable than any fixed HDI when the question is directional ("does the policy help or hurt?").

---

Audience adaptation

Technical audience (researchers, analysts)

• Full DAG with node and edge justifications
• Formal identification result (backdoor/frontdoor/IV/RDD/DiD) with citation if applicable
• Posterior plots and full diagnostic summary
• Complete refutation table with test statistics
• Sensitivity analysis (E-values, partial R² bounds, or Rosenbaum bounds)
• Code or link to code repository in appendix

Decision-makers (executives, policymakers)

• Causal question in plain language — one sentence
• Effect size translated to natural frequencies: "For every 100 people exposed, we estimate 8 more would [outcome]"
• Key threats in 1–2 sentences, stated simply: "The main reason this estimate could be wrong is [X]. If so, the true effect is likely [smaller/larger]."
• Actionable recommendation with explicit uncertainty: "Given the uncertainty, we recommend [action] if the cost of a false positive is less than [threshold]."
• Technical details (DAG, model specification, diagnostics, refutation table) in a clearly labeled appendix

Both audiences get the limitations section. Never hide limitations from decision-makers on the grounds that they are "too technical." Translate, do not omit.

Presenting to mixed audiences

Open with the plain-language conclusion and effect size. Then walk through the DAG visually — most people understand arrows even without training. Reserve equations and diagnostic plots for Q&A or written appendices.

---

Common reporting mistakes

1. Using causal language without identification. If there is no identification strategy, the word "effect" is wrong. Use "association."

2. Reporting only the point estimate. The posterior is the result. Always show the full distribution or at minimum multiple HDIs.

3. Hiding refutation failures. A failed refutation is information. Report it, downgrade your language, and explain what the failure means for the conclusion.

4. Burying limitations. Threats to validity belong in the body of the report, ranked by severity — not in an appendix labeled "caveats."

5. Conflating LATE with ATE. IV estimates give the Local Average Treatment Effect for compliers only. DiD estimates are often ATT (Average Treatment Effect on the Treated). Be explicit about whose effect you are estimating and whether it answers the question.

6. Ignoring spillovers. If SUTVA is violated — units affect each other — the estimated effect conflates direct and spillover effects. State whether spillovers are plausible and, if so, what direction they push the estimate.

7. Omitting the E-value or sensitivity analysis. For observational studies, always quantify how much unobserved confounding would be needed to overturn the conclusion. This anchors the limitations discussion in something concrete rather than vague hedging.

Structural Causal Models

Contents

• When to use structural models
• pm.do() for interventions
• pm.observe() for conditioning
• Counterfactual queries
• Mediation analysis
• Structural vs. quasi-experimental: decision guide

---

When to use structural models

Choose a structural causal model (SCM) when:

• You have a full causal theory and want to model mechanisms explicitly
• The problem does not fit a quasi-experimental template (no natural experiment, no discontinuity, no instrument)
• You need counterfactual queries: "what would have happened to this specific unit if X had been different?"
• You need to decompose effects into direct and indirect paths (mediation)

Advantage: flexible, answers counterfactuals, decomposes effects, transparent about assumptions via the DAG.

Disadvantage: requires specifying a full structural model — more assumptions, more ways to be wrong. If a quasi-experimental design is available, prefer it; fewer assumptions is better.

---

pm.do() for interventions

pm.do() implements the do-calculus: it cuts incoming edges to the intervened variable, simulating an ideal randomized experiment on the existing model graph.

import pymc as pm
import numpy as np

RANDOM_SEED = sum(map(ord, "causal-structural-v1"))
rng = np.random.default_rng(RANDOM_SEED)

# 1. Define the generative (observational) model
with pm.Model() as scm:
    Z = pm.Normal("Z", mu=0, sigma=1)                    # confounder
    X = pm.Normal("X", mu=0.5 * Z, sigma=0.5)            # treatment (caused by Z)
    Y = pm.Normal("Y", mu=0.3 * X + 0.7 * Z, sigma=0.5) # outcome

# 2. Intervene: set X = 1 (do-calculus)
scm_do_1 = pm.do(scm, {"X": 1})
with scm_do_1:
    idata_do_1 = pm.sample_prior_predictive(samples=2000, random_seed=rng)

# 3. Compare do(X=1) vs do(X=0) for Average Treatment Effect
scm_do_0 = pm.do(scm, {"X": 0})
with scm_do_0:
    idata_do_0 = pm.sample_prior_predictive(samples=2000, random_seed=rng)

ate = idata_do_1.prior["Y"].mean() - idata_do_0.prior["Y"].mean()
print(f"ATE: {ate.values:.3f}")  # Should recover ~0.3 (the direct X -> Y coefficient)

Note: Using sample_prior_predictive gives the ATE under the prior (forward-simulating from structural equations). For posterior-based counterfactuals — using observed data to sharpen estimates — use pm.observe() first, then pm.do().

---

pm.observe() for conditioning

pm.observe() conditions the model on observed data, fixing variables to their measured values. This is the standard way to incorporate data into an SCM before doing counterfactual reasoning.

import pymc as pm

# Condition on observed outcome to infer latent variables
scm_obs = pm.observe(scm, {"Y": y_observed})
with scm_obs:
    idata_obs = pm.sample(nuts_sampler="nutpie", random_seed=rng)

Use pm.observe() when you want to:

• Fit the structural model to data (standard Bayesian inference)
• Perform the abduction step of a counterfactual query (see below)
• Infer latent confounders from observed variables

---

Counterfactual queries

A counterfactual query asks: "What would Y have been for this specific unit if X had been different?" This is a unit-level question, not a population average — quasi-experimental methods cannot answer it.

The answer requires the three-step abduction-action-prediction procedure:

Step 1 — Abduction: condition on the unit's factual observations to infer their latent factors.

Step 2 — Action: intervene on the treatment using pm.do().

Step 3 — Prediction: forward-simulate the outcome under the counterfactual treatment, using the inferred latent factors from Step 1.

# Step 1: Abduction — infer this unit's latent confounder Z
scm_unit = pm.observe(scm, {"X": x_factual, "Y": y_factual})
with scm_unit:
    idata_abduction = pm.sample(nuts_sampler="nutpie", random_seed=rng)

# Steps 2 & 3: Action + Prediction
# Use the posterior of Z (the inferred latent noise) to forward-simulate Y
# under the counterfactual treatment value
z_posterior = idata_abduction.posterior["Z"]
y_counterfactual = 0.3 * x_counterfactual + 0.7 * z_posterior

print(f"Counterfactual Y mean: {float(y_counterfactual.mean()):.3f}")

Critical assumption: the structural equations and noise distributions are the same in the factual and counterfactual worlds. If your model is misspecified, counterfactuals inherit that misspecification.

---

Mediation analysis

Mediation decomposes the total treatment effect into:

• Total Effect (TE): the full effect of T on Y through all paths
• Natural Direct Effect (NDE): the effect of T on Y that does not pass through the mediator M
• Natural Indirect Effect (NIE): the effect of T on Y that operates through M

Additive decomposition: TE ≈ NDE + NIE

# Schematic DAG: T -> Y (direct), T -> M -> Y (indirect)
# NDE: intervene do(T=1) but hold M fixed at its value under do(T=0)
# NIE: intervene do(T=0) but let M take its value under do(T=1)
#
# Preferred: use pm.do() to compute interventional distributions for each path.
# Acceptable alternative: Bayesian SEM with chained structural equations, computing
# NDE/NIE from posterior draws of structural coefficients (a * b for indirect, c for direct).
# pm.do() is preferred because it makes the interventional logic explicit and avoids
# manual algebra on coefficients, but both approaches give equivalent results when
# the structural model is correctly specified.

with pm.Model() as mediation_scm:
    T = pm.Bernoulli("T", p=0.5)
    M = pm.Normal("M", mu=0.6 * T, sigma=0.5)       # mediator
    Y = pm.Normal("Y", mu=0.4 * T + 0.5 * M, sigma=0.5)  # outcome

# NDE: fix M to its T=0 distribution, intervene T=1
m_under_t0 = pm.do(mediation_scm, {"T": 0})
# ... sample M under T=0, then use that M in Y equation under T=1
# Full implementation depends on model structure; use pm.do() to block/unblock paths

Strong assumptions for mediation — make these explicit:

1. No unmeasured T-M confounding 2. No unmeasured M-Y confounding 3. No unmeasured T-Y confounding 4. No effect of T on M-Y confounders

Violation of any of these invalidates the NDE/NIE decomposition. Always communicate these assumptions to the user before reporting mediated effects.

---

Structural vs. quasi-experimental: decision guide

When in doubt, ask: "Do I have a credible natural experiment?" If yes, start there. If no, build the SCM and be explicit about every assumption in the DAG.