Code Snippets

This notebook contains all code snippets that appear in the forthcoming paper submission, available at [TODO].

from ananke import graphs, identification
from ananke.models import binary_nested
from ananke.estimation import CausalEffect
from ananke.datasets import load_wisconsin_health_study, load_frontdoor_data

Identification

vertices = ["A", "M", "Y"]
di_edges = [("A", "M"), ("M", "Y")]
bi_edges = [("A", "Y")]
front_door = graphs.ADMG(vertices, di_edges, bi_edges)
treatments = ["A"]
outcomes = ["Y"]
id_front_door = identification.OneLineID(graph=front_door, treatments=treatments, outcomes=outcomes)
print('Identifed =', id_front_door.id(), '; Functional =', id_front_door.functional())

{('Y',): ['M', 'A'], ('M',): ['Y', 'A']}
Identifed = True ; Functional = ΣM ΦAY(p(V);G) ΦAM(p(V);G)

Estimation through Mobius Parameterization

df = load_wisconsin_health_study().rename(columns={"E": "A"}).groupby(["A", "M", "Y"])["count"].sum().reset_index()
bnm = binary_nested.BinaryNestedModel(front_door)
bnm = bnm.fit(X=df, tol=1e-12)
theta_M0_A0 = bnm.fitted_params[(frozenset({'M'}), ('A',), (0,))]
theta_M0_A1 = bnm.fitted_params[(frozenset({'M'}), ('A',), (1,))]
theta_Y0_M0 = bnm.fitted_params[(frozenset({'Y'}), ('M',), (0,))]
theta_Y0_M1 = bnm.fitted_params[(frozenset({'Y'}), ('M',), (1,))]
pY1_A0 = bnm.estimate(treatment_dict={"A": 0}, outcome_dict={"Y": 1})
pY1_A1 = bnm.estimate(treatment_dict={"A": 1}, outcome_dict={"Y": 1})
print("ACE: ", pY1_A1 - pY1_A0)
# ACE: 0.004

/Users/jaron/Projects/ananke/ananke/models/binary_nested.py:507: RuntimeWarning: invalid value encountered in log
  logger.debug(np.log(A @ x - b))
/Users/jaron/Projects/ananke/ananke/models/binary_nested.py:508: RuntimeWarning: invalid value encountered in log
  lld = np.sum(counts * np.log(A @ x - b))

ACE:  0.004727830873286432

Estimation using Efficient Semiparametric Methods

data = load_frontdoor_data().rename(columns={"T": "A", "ViralLoad2": "M", "CD4": "Y"})[["A", "M", "Y"]]

ace_obj = CausalEffect(graph=front_door, treatment='A', outcome='Y')  # setting up the ACE object
ace = ace_obj.compute_effect(data, "eff-apipw")
print("ace = ", ace)

 Treatment is p-fixable and graph is mb-shielded.

 Available estimators are:

1. Primal IPW (p-ipw)
2. Dual IPW (d-ipw)
3. APIPW (apipw)
4. Efficient APIPW (eff-apipw)

Suggested estimator is Efficient APIPW

ace =  -0.9020293902163798