CausalEngine Tutorial

import numpy as np
import pandas as pd
import causalscience as ce
from causalscience.indexed_graph import Slot, ObsKernel, Product, Sum

assert ce.__version__ == "0.1.0a13"
rng = np.random.default_rng(27)
n = 2400
common = dict(seed=0, num_boost_round=30, params={"nthread": 1, "max_depth": 4})

print("CausalEngine:", ce.__version__)
print("NumPy:", np.__version__, "pandas:", pd.__version__)
print("Data seed: 27; rows per example:", n)
print("Estimator seed:", common["seed"], "boosting rounds:", common["num_boost_round"])

# Back-door mean: baseline Z and C suffice for adjustment in this SCM.
z = rng.binomial(1, .5, n)
c = rng.binomial(1, .2 + .5*z)
a = rng.binomial(1, .1 + .3*c + .2*z)
y = rng.binomial(1, .05 + .15*c + .1*z + a*(.2 + .2*z))
bd_data = pd.DataFrame(dict(Z=z, C=c, A=a, Y=y))
bd_schema = {v: "binary" for v in bd_data}
bd_options = dict(treatment="A", outcome="Y", adjustment=["Z", "C"], schema=bd_schema)
bd = ce.estimate_bd(bd_data, action=1, **bd_options, **common)

print("Status:", bd["status"], "/", bd["estimation_status"])
print("Interventional mean (A=1):", bd["estimate"])

# Sequential back-door: new blocks, in time order; L depends on A0.
c = rng.binomial(1, .5, n)
a0 = rng.binomial(1, .2 + .4*c)
l = rng.binomial(1, .1 + .2*c + .5*a0)
a1 = rng.binomial(1, .1 + .1*c + .6*l + .1*a0)
y = rng.binomial(1, .05 + .15*c + .3*l + .1*a0 + .25*a1)
sbd_data = pd.DataFrame(dict(C=c, A0=a0, L=l, A1=a1, Y=y))
sbd = ce.estimate_sbd(sbd_data, treatments=["A0", "A1"], outcome="Y",
    covariates=[["C"], ["L"]], actions={"A0": 1, "A1": 1},
    schema={v: "binary" for v in sbd_data}, **common)

print("Status:", sbd["status"], "/", sbd["estimation_status"])
print("Interventional mean (A0=1, A1=1):", sbd["estimate"])

# Conditional front-door: C enters every conditional and the outer average.
c = rng.binomial(1, .5, n)
u = rng.binomial(1, .5, n)
a = rng.binomial(1, .1 + .3*c + .4*u)
m = rng.binomial(1, .1 + .5*a + .25*c)
y = rng.binomial(1, .05 + .4*m + .25*c + .2*u)
fd_data = pd.DataFrame(dict(C=c, A=a, M=m, Y=y))
fd = ce.estimate_fd(fd_data, treatment="A", mediators=["M"], outcome="Y",
    covariates=["C"], action=1, schema={v: "binary" for v in fd_data}, **common)

print("Status:", fd["status"], "/", fd["estimation_status"])
print("Interventional mean (A=1):", fd["estimate"])

# UCA: explicitly supply a normalized expression, with no graph object.
Z, C, A, Y = map(Slot, ("Z", "C", "A", "Y"))
formula = Sum((Z, C), Product((
    ObsKernel((("Z", Z), ("C", C))),
    ObsKernel((("Y", Y),), (("A", A), ("Z", Z), ("C", C))),
)))
uca = ce.estimate_uca(bd_data, formula=formula, actions={"A": 1}, outcome="Y",
    schema=bd_schema, **common)
np.testing.assert_array_equal(uca["estimate"], bd["estimate"])

print("Status:", uca["status"], "/", uca["estimation_status"])
print("Interventional mean (A=1):", uca["estimate"])

# ATE: paired action 1 minus 0, retaining covariance of the two scores.
ate = ce.estimate_ate(bd_data, estimator="bd", actions={"A": 1},
    reference_actions={"A": 0}, **bd_options, **common)

print("Status:", ate["status"], "/", ate["estimation_status"])
print("Paired ATE (A=1 minus A=0):", ate["estimate"])

# CATE: regress paired pseudo-outcomes on the baseline modifier Z.
cate = ce.estimate_cate(bd_data, estimator="bd", modifiers=["Z"],
    actions={"A": 1}, reference_actions={"A": 0}, **bd_options, **common,
    cate_num_boost_round=30, cate_params={"nthread": 1, "max_depth": 4})

print("Status:", cate["status"], "/", cate["estimation_status"])
print("Marginal paired ATE summary (not a conditional CATE value):", cate["estimate"])

for name, result in dict(BD=bd, SBD=sbd, FD=fd, UCA=uca, ATE=ate, CATE=cate).items():
    assert result["estimation_status"] == "ACCEPT", result.get("reason")
    assert result["graph_verification"] == "not_run"
    assert result["causal_interval"] is None
print("Conditional CATE predictions at Z=0,1:", cate["model"].predict(pd.DataFrame({"Z": [0, 1]})))