Causal Inference

Overview

Causal inference determines cause-and-effect relationships and estimates treatment effects, going beyond correlation to understand what causes what.

When to Use

• Evaluating the impact of policy interventions or business decisions

• Estimating treatment effects when randomized experiments aren't feasible

• Controlling for confounding variables in observational data

• Determining if a marketing campaign or product change caused an outcome

• Analyzing heterogeneous treatment effects across different user segments

• Making causal claims from non-experimental data using propensity scores or instrumental variables

Key Concepts

• Treatment: Intervention or exposure

• Outcome: Result or consequence

• Confounding: Variables affecting both treatment and outcome

• Causal Graph: Visual representation of relationships

• Treatment Effect: Impact of intervention

• Selection Bias: Non-random treatment assignment

Causal Methods

• Randomized Controlled Trials (RCT): Gold standard

• Propensity Score Matching: Balance treatment/control

• Difference-in-Differences: Before/after comparison

• Instrumental Variables: Handle endogeneity

• Causal Forests: Heterogeneous treatment effects

Implementation with Python

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

import seaborn as sns

from sklearn.linear_model import LinearRegression, LogisticRegression

from sklearn.preprocessing import StandardScaler

from scipy import stats

# Generate observational data with confounding

np.random.seed(42)

n = 1000

# Confounder: Age (affects both treatment and outcome)

age = np.random.uniform(25, 75, n)

# Treatment: Training program (more likely for younger people)

treatment_prob = 0.3 + 0.3 * (75 - age) / 50  # Inverse relationship with age

treatment = (np.random.uniform(0, 1, n) < treatment_prob).astype(int)

# Outcome: Salary (affected by both treatment and age)

# True causal effect of treatment: +$5000

salary = 40000 + 500 * age + 5000 * treatment + np.random.normal(0, 10000, n)

df = pd.DataFrame({

    'age': age,

    'treatment': treatment,

    'salary': salary,

})

print("Observational Data Summary:")

print(df.describe())

print(f"\nTreatment Rate: {df['treatment'].mean():.1%}")

print(f"Average Salary (Control): ${df[df['treatment']==0]['salary'].mean():.0f}")

print(f"Average Salary (Treatment): ${df[df['treatment']==1]['salary'].mean():.0f}")

# 1. Naive Comparison (BIASED - ignores confounding)

naive_effect = df[df['treatment']==1]['salary'].mean() - df[df['treatment']==0]['salary'].mean()

print(f"\n1. Naive Comparison: ${naive_effect:.0f} (BIASED)")

# 2. Regression Adjustment (Covariate Adjustment)

X = df[['treatment', 'age']]

y = df['salary']

model = LinearRegression()

model.fit(X, y)

regression_effect = model.coef_[0]

print(f"\n2. Regression Adjustment: ${regression_effect:.0f}")

# 3. Propensity Score Matching

# Estimate probability of treatment given covariates

ps_model = LogisticRegression()

ps_model.fit(df[['age']], df['treatment'])

df['propensity_score'] = ps_model.predict_proba(df[['age']])[:, 1]

print(f"\n3. Propensity Score Matching:")

print(f"PS range: [{df['propensity_score'].min():.3f}, {df['propensity_score'].max():.3f}]")

# Matching: find control for each treated unit

matched_pairs = []

treated_units = df[df['treatment'] == 1].index

for treated_idx in treated_units:

    treated_ps = df.loc[treated_idx, 'propensity_score']

    treated_age = df.loc[treated_idx, 'age']

    # Find closest control unit

    control_units = df[(df['treatment'] == 0) &#x26;

                      (df['propensity_score'] >= treated_ps - 0.1) &#x26;

                      (df['propensity_score'] <= treated_ps + 0.1)].index

    if len(control_units) > 0:

        closest_control = min(control_units,

                             key=lambda x: abs(df.loc[x, 'propensity_score'] - treated_ps))

        matched_pairs.append({

            'treated_idx': treated_idx,

            'control_idx': closest_control,

            'treated_salary': df.loc[treated_idx, 'salary'],

            'control_salary': df.loc[closest_control, 'salary'],

        })

matched_df = pd.DataFrame(matched_pairs)

psm_effect = (matched_df['treated_salary'] - matched_df['control_salary']).mean()

print(f"PSM Effect: ${psm_effect:.0f}")

print(f"Matched pairs: {len(matched_df)}")

# 4. Stratification by Propensity Score

df['ps_stratum'] = pd.qcut(df['propensity_score'], q=5, labels=False, duplicates='drop')

stratified_effects = []

for stratum in df['ps_stratum'].unique():

    stratum_data = df[df['ps_stratum'] == stratum]

    if (stratum_data['treatment'] == 0).sum() > 0 and (stratum_data['treatment'] == 1).sum() > 0:

        treated_mean = stratum_data[stratum_data['treatment'] == 1]['salary'].mean()

        control_mean = stratum_data[stratum_data['treatment'] == 0]['salary'].mean()

        effect = treated_mean - control_mean

        stratified_effects.append(effect)

stratified_effect = np.mean(stratified_effects)

print(f"\n4. Stratification by PS: ${stratified_effect:.0f}")

# 5. Visualization

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Treatment distribution by age

ax = axes[0, 0]

treated = df[df['treatment'] == 1]

control = df[df['treatment'] == 0]

ax.hist(control['age'], bins=20, alpha=0.6, label='Control', color='blue')

ax.hist(treated['age'], bins=20, alpha=0.6, label='Treated', color='red')

ax.set_xlabel('Age')

ax.set_ylabel('Frequency')

ax.set_title('Age Distribution by Treatment')

ax.legend()

ax.grid(True, alpha=0.3, axis='y')

# Salary vs Age (colored by treatment)

ax = axes[0, 1]

ax.scatter(control['age'], control['salary'], alpha=0.5, label='Control', s=30)

ax.scatter(treated['age'], treated['salary'], alpha=0.5, label='Treated', s=30, color='red')

ax.set_xlabel('Age')

ax.set_ylabel('Salary')

ax.set_title('Salary vs Age by Treatment')

ax.legend()

ax.grid(True, alpha=0.3)

# Propensity Score Distribution

ax = axes[1, 0]

ax.hist(df[df['treatment'] == 0]['propensity_score'], bins=20, alpha=0.6, label='Control', color='blue')

ax.hist(df[df['treatment'] == 1]['propensity_score'], bins=20, alpha=0.6, label='Treated', color='red')

ax.set_xlabel('Propensity Score')

ax.set_ylabel('Frequency')

ax.set_title('Propensity Score Distribution')

ax.legend()

ax.grid(True, alpha=0.3, axis='y')

# Treatment Effect Comparison

ax = axes[1, 1]

methods = ['Naive', 'Regression', 'PSM', 'Stratified']

effects = [naive_effect, regression_effect, psm_effect, stratified_effect]

true_effect = 5000

ax.bar(methods, effects, color=['red', 'orange', 'yellow', 'lightgreen'], alpha=0.7, edgecolor='black')

ax.axhline(y=true_effect, color='green', linestyle='--', linewidth=2, label=f'True Effect (${true_effect:.0f})')

ax.set_ylabel('Treatment Effect ($)')

ax.set_title('Treatment Effect Estimates by Method')

ax.legend()

ax.grid(True, alpha=0.3, axis='y')

for i, effect in enumerate(effects):

    ax.text(i, effect + 200, f'${effect:.0f}', ha='center', va='bottom')

plt.tight_layout()

plt.show()

# 6. Doubly Robust Estimation

from sklearn.ensemble import RandomForestRegressor

# Propensity score model

ps_model_dr = LogisticRegression().fit(df[['age']], df['treatment'])

ps_scores = ps_model_dr.predict_proba(df[['age']])[:, 1]

# Outcome model

outcome_model = RandomForestRegressor(n_estimators=50, random_state=42)

outcome_model.fit(df[['treatment', 'age']], df['salary'])

# Doubly robust estimator

treated_mask = df['treatment'] == 1

control_mask = df['treatment'] == 0

# Adjust for propensity score

treated_adjusted = (treated_mask.astype(int) * df['salary']) / (ps_scores + 0.01)

control_adjusted = (control_mask.astype(int) * df['salary']) / (1 - ps_scores + 0.01)

# Outcome predictions

pred_treated = outcome_model.predict(df[['treatment', 'age']].replace({'treatment': 0, 1: 1}))

pred_control = outcome_model.predict(df[['treatment', 'age']].replace({'treatment': 1, 0: 0}))

dr_effect = treated_adjusted.sum() / treated_mask.sum() - control_adjusted.sum() / control_mask.sum()

print(f"\n6. Doubly Robust Estimation: ${dr_effect:.0f}")

# 7. Heterogeneous Treatment Effects

print(f"\n7. Heterogeneous Treatment Effects (by Age Quartile):")

for age_q in pd.qcut(df['age'], q=4, duplicates='drop').unique():

    mask = (df['age'] >= age_q.left) &#x26; (df['age'] < age_q.right)

    stratum_data = df[mask]

    if (stratum_data['treatment'] == 0).sum() > 0 and (stratum_data['treatment'] == 1).sum() > 0:

        treated_mean = stratum_data[stratum_data['treatment'] == 1]['salary'].mean()

        control_mean = stratum_data[stratum_data['treatment'] == 0]['salary'].mean()

        effect = treated_mean - control_mean

        print(f"  Age {age_q.left:.0f}-{age_q.right:.0f}: ${effect:.0f}")

# 8. Sensitivity Analysis

print(f"\n8. Sensitivity Analysis (Hidden Confounder Impact):")

# Vary hidden confounder correlation with outcome

for hidden_effect in [1000, 2000, 5000, 10000]:

    adjusted_effect = regression_effect - hidden_effect * 0.1

    print(f"  If hidden confounder worth ${hidden_effect}: Effect = ${adjusted_effect:.0f}")

# 9. Summary Table

print(f"\n" + "="*60)

print("CAUSAL INFERENCE SUMMARY")

print("="*60)

print(f"True Treatment Effect: ${true_effect:,.0f}")

print(f"\nEstimates:")

print(f"  Naive (BIASED): ${naive_effect:,.0f}")

print(f"  Regression Adjustment: ${regression_effect:,.0f}")

print(f"  Propensity Score Matching: ${psm_effect:,.0f}")

print(f"  Stratification: ${stratified_effect:,.0f}")

print(f"  Doubly Robust: ${dr_effect:,.0f}")

print("="*60)

# 10. Causal Graph (Text representation)

print(f"\n10. Causal Graph (DAG):")

print(f"""

Age → Treatment ← (Selection Bias)

  ↓        ↓

  └─→ Salary

Interpretation:

- Age is a confounder

- Treatment causally affects Salary

- Age directly affects Salary

- Age affects probability of Treatment

""")

Causal Assumptions

• Unconfoundedness: No unmeasured confounders

• Overlap: Common support on propensity scores

• SUTVA: No interference between units

• Consistency: Single version of treatment

Treatment Effect Types

• ATE: Average Treatment Effect (overall)

• ATT: Average Treatment on Treated

• CATE: Conditional Average Treatment Effect

• HTE: Heterogeneous Treatment Effects

Method Strengths

• RCT: Gold standard, controls all confounders

• Matching: Balances groups, preserves overlap

• Regression: Adjusts for covariates

• Instrumental Variables: Handles endogeneity

• Causal Forests: Learns heterogeneous effects

Deliverables

• Causal graph visualization

• Treatment effect estimates

• Sensitivity analysis

• Heterogeneous treatment effects

• Covariate balance assessment

• Propensity score diagnostics

• Final causal inference report