GuidesFit a structural causal model

Fit a structural causal model

Fit a Bayesian structural causal model (SCM) over a graph and ask it for an average effect. You'll also see what it returns when the graph says the data can't answer.

Run the program#

The data is simulated, so you know the answer. A hidden variable U drives both X and Y, and X affects Y only through M. The true effect of X on Y is 0.8 × 1.5 = 1.2.

import numpy as np
import pandas as pd
from ergodic import SCM, admg, dag, tabular

# 1. Simulate 500 units. A hidden U drives X and Y; X acts on Y only through M.
rng = np.random.default_rng(0)
n = 500
u = rng.normal(size=n)
x = u + rng.normal(size=n)
m = 0.8 * x + rng.normal(size=n)
y = 1.5 * m + 2.0 * u + rng.normal(size=n)
data = tabular(pd.DataFrame({"X": x, "M": m, "Y": y}))
print(f"true effect of X on Y: {0.8 * 1.5:.2f}")
print(f"slope of Y on X:       {np.polyfit(x, y, 1)[0]:.2f}")

# 2. Draw the graph. X <-> Y says that a hidden common cause links X and Y.
graph = admg(["X -> M", "M -> Y", "X <-> Y"])

# 3. Fit one mechanism per node.
scm = SCM(graph, data, draws=500, tune=500, chains=1, random_state=0).fit()
for node, mechanism in scm.mechanisms.items():
    print(f"{node}: {type(mechanism).__name__}, parents {mechanism.parents}")

# 4. Ask for the average effect of a one-unit increase in X.
effect = scm.ate("X", "Y")
print(effect)

# 5. Drop the mediator from the graph. Now no formula identifies the effect.
bow = dag(["U -> X", "U -> Y", "X -> Y"], latent=["U"])
bounds = SCM(bow, data, draws=500, tune=500, chains=1, random_state=0).fit().ate("X", "Y")
low, high = bounds.range()
print(f"{type(bounds).__name__}: effect between {low:.2f} and {high:.2f}")
print(f"effect if the hidden path carries covariance 2.0: {bounds.at(2.0):.2f}")

true effect of X on Y: 1.20
slope of Y on X:       2.32
X: GaussianAdditive, parents ()
M: GaussianAdditive, parents ('X',)
Y: GaussianAdditive, parents ('M',)
SCMResult(ate, identified: Y=1.265, sd=0.0772, 95% CrI [1.12, 1.42])
Sensitivity: effect between 0.00 and 2.32
effect if the hidden path carries covariance 2.0: 1.28

• Step 1. The raw slope is almost twice the true effect. The hidden U inflates it.
• Step 3. Each continuous node gets a GaussianAdditive mechanism: a Bayesian linear equation plus its own noise term.
• Step 4. The graph allows the front-door formula, so ate returns an SCMResult: a posterior mean of about 1.27, with a 95% credible interval that contains the true 1.2.
• Step 5. Without M, the effect isn't identified. ate returns a Sensitivity instead of a number: the effect for each assumed strength of the hidden path.

Read the two return types#

ate asks identify_effect first, and the answer decides the return type:

• Identified by adjustment, by the front-door formula, or because no directed path exists: an SCMResult. Read .value, .std_error, .interval(level) and .draws().
• Not identified: a Sensitivity. It has no .value, so you can't read a single number by mistake.

For a continuous treatment and outcome, a Sensitivity holds the line b(c) = (cov(X, Y) − c) / var(X), where c is the covariance the hidden path adds. Read it with .range(), .at(c) and the .confounding and .effect arrays:

low, high = bounds.range()
print(f"{type(bounds).__name__}: effect between {low:.2f} and {high:.2f}")
print(f"effect if the hidden path carries covariance 2.0: {bounds.at(2.0):.2f}")

The range runs from no confounding (c = 0, the raw slope) to all of the covariance being confounding (an effect of zero). Here U adds a covariance of 2.0, so .at(2.0) lands near the truth. With real data you don't know c.

The range covers the truth only if the hidden path pushes in the same direction as the observed association, and adds no more than all of it. For any other c, .at(c) extrapolates exactly along the line.

For a binary or categorical treatment with a continuous outcome, the band comes from refitting the model under a grid of confounding priors instead. Then .confounding holds prior scales, not covariances.

Choose the mechanisms#

SCM assigns a mechanism to each node from its column type:

• continuous: GaussianAdditive;
• binary: BinaryGumbelMax, a logistic choice with explicit noise;
• categorical: CategoricalGumbelMax.

To override one, pass mechanisms={"Y": "linear"}. The names are "linear", "binary" and "categorical". You can also pass a mechanism instance from ergodic.scm, for example GaussianAdditive(name="Y", parents=("M",), coef_prior_sd=1.0).

draws, tune and chains set the sampler budget. The defaults are 1,000 draws, 1,000 tuning steps and 2 chains. This page uses a smaller budget so it runs in seconds. cores (default 1) caps the processes PyMC may use for the chains, so the chains run one after another unless you raise it; the draws are the same either way.

Next steps#