leonardodalinky/scider/causal-inference

# Simple DAG visualization with networkx
import networkx as nx
import matplotlib.pyplot as plt
G = nx.DiGraph()
G.add_edges_from([
    ("Education", "Income"),      # T → Y (causal path)
    ("Family_SES", "Education"),  # C → T (confounder)
    ("Family_SES", "Income"),     # C → Y (confounder)
])
pos = nx.spring_layout(G, seed=42)
nx.draw(G, pos, with_labels=True, node_color="lightblue",
        node_size=2000, arrows=True, arrowsize=20)

import pandas as pd
from scipy import stats
# Check covariate balance between treatment groups
for covariate in ["age", "income", "education"]:
    t_stat, p_val = stats.ttest_ind(
        df[df.treated == 1][covariate],
        df[df.treated == 0][covariate]
    )
    print(f"{covariate}: p={p_val:.3f} {'⚠️ IMBALANCED' if p_val < 0.05 else '✅'}")

from linearmodels.iv import IV2SLS
# Two-Stage Least Squares
res = IV2SLS.from_formula(
    "outcome ~ 1 + controls + [treatment ~ instrument]",
    data=df
).fit(cov_type="robust")
print(res.summary)
# Check first-stage relevance (F-stat > 10 rule of thumb)
from linearmodels.iv.model import _OLS
first_stage = _OLS.from_formula("treatment ~ 1 + controls + instrument", data=df).fit()
print(f"First-stage F-stat: {first_stage.f_statistic.stat:.2f}")

import statsmodels.formula.api as smf
# Two-period DiD
# treated: dummy for treatment group; post: dummy for post-treatment period
result = smf.ols(
    "outcome ~ treated + post + treated:post + controls",
    data=df
).fit(cov_type="HC3")
did_estimate = result.params["treated:post"]
print(f"DiD estimate: {did_estimate:.3f}")
print(result.summary())
# Parallel trends test: run DiD only on pre-treatment data
# The interaction should be near zero
pre_data = df[df.post == 0]
# ... (use period dummies instead of post)

from linearmodels.panel import PanelOLS
result = PanelOLS.from_formula(
    "outcome ~ treatment + EntityEffects + TimeEffects",
    data=df.set_index(["unit_id", "time_id"])
).fit(cov_type="clustered", cluster_entity=True)

import numpy as np
import statsmodels.formula.api as smf
cutoff = 50  # assignment threshold
df["above_cutoff"] = (df.running_var >= cutoff).astype(int)
df["centered"] = df.running_var - cutoff
# Sharp RDD: linear fit on each side of cutoff
result = smf.ols(
    "outcome ~ centered * above_cutoff",
    data=df[np.abs(df.centered) <= 10]  # bandwidth = 10
).fit(cov_type="HC3")
rdd_estimate = result.params["above_cutoff"]
print(f"RDD estimate at cutoff: {rdd_estimate:.3f}")
# Validity checks:
# 1. Continuity of density at cutoff (McCrary test — rddensity package)
# 2. Continuity of pre-determined covariates at cutoff (run RDD with covariate as outcome)

from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
import pandas as pd
import numpy as np
# Step 1: Estimate propensity score
X = df[["age", "income", "education", "prior_outcome"]]
y = df["treated"]
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
ps_model = LogisticRegression(max_iter=1000)
ps_model.fit(X_scaled, y)
df["propensity_score"] = ps_model.predict_proba(X_scaled)[:, 1]
# Step 2: Check overlap (common support)
import matplotlib.pyplot as plt
for t in [0, 1]:
    plt.hist(df[df.treated == t]["propensity_score"],
             alpha=0.5, label=f"Treated={t}", bins=30)
plt.xlabel("Propensity Score")
plt.legend()
plt.title("Common Support Check")
# Step 3: 1:1 nearest-neighbor matching
from scipy.spatial.distance import cdist
treated = df[df.treated == 1]
control = df[df.treated == 0]
distances = cdist(
    treated[["propensity_score"]].values,
    control[["propensity_score"]].values
)
matched_control_idx = distances.argmin(axis=1)
matched = pd.concat([
    treated.reset_index(drop=True),
    control.iloc[matched_control_idx].reset_index(drop=True).add_suffix("_ctrl")
], axis=1)
# Step 4: Check covariate balance after matching
att_estimate = (matched["outcome"] - matched["outcome_ctrl"]).mean()
print(f"ATT estimate: {att_estimate:.3f}")

# With dowhy (recommended for production analysis)
import dowhy
model = dowhy.CausalModel(
    data=df,
    treatment="treated",
    outcome="outcome",
    common_causes=["age", "income", "education"],
)
identified = model.identify_effect()
estimate = model.estimate_effect(identified, method_name="backdoor.propensity_score_matching")
refute = model.refute_estimate(identified, estimate, method_name="placebo_treatment_refuter")

# E-value formula for relative risks (VanderWeele & Ding, 2017)
def e_value(rr):
    """E-value for a relative risk estimate."""
    return rr + np.sqrt(rr * (rr - 1))
rr = 2.5  # observed relative risk
print(f"E-value: {e_value(rr):.2f}")
# Interpretation: an unmeasured confounder would need to have a risk ratio ≥ E-value
# with both treatment and outcome to explain away the result