- DoWhy separates identification (graph-based) from estimation (statistical)—understanding this distinction is critical for valid causal claims.
- IPW estimators fail silently when propensity scores hit extremes; the 0.01-0.99 clipping introduces bias but prevents infinite weights.
- Doubly robust estimation gives you two chances to get it right: consistent if either your propensity model or outcome model is correctly specified.
- Refutation tests check robustness, not correctness—passing all refutations doesn't validate your causal graph.
- Building a minimal causal engine teaches you why positivity and identification matter in ways documentation can't convey.
Most “Causal” ML Is Just Correlation with Extra Steps
Here’s a take that’ll ruffle some feathers: 90% of production “causal inference” I’ve seen is regression with a fancier name. Teams slap on DoWhy, run estimate_effect(), and ship whatever number falls out—without understanding what the library actually does under the hood.
The result? Causal claims built on sand.
I’m not saying DoWhy is bad. It’s genuinely excellent. But treating it as a black box defeats the entire purpose. The power of causal inference comes from making your assumptions explicit—and you can’t do that if you don’t understand what assumptions DoWhy is making for you.
So let’s build a minimal causal inference engine from scratch, then reverse-engineer DoWhy’s internals to see how the real thing works.
The Four-Step Pipeline That DoWhy Actually Runs
DoWhy follows a workflow that looks deceptively simple:
Model → Identify → Estimate → Refute
But each arrow hides substantial complexity. Here’s what’s actually happening.
Model: You specify a directed acyclic graph (DAG) representing your causal assumptions. Treatment causes outcome. Confounders cause both. Maybe there’s a mediator in between. DoWhy parses this graph using NetworkX.
Identify: Given the graph, DoWhy runs path-finding algorithms to determine which causal effect is identifiable and how to identify it. This is the most important step—and the one most people skip past.
Estimate: Apply a statistical estimator (IPW, matching, regression, etc.) to actually compute the number.
Refute: Run sensitivity analyses to check if your estimate is robust.
The critical insight is that identification is purely graph-based. It doesn’t look at your data at all. If your graph is wrong, no amount of sophisticated estimation will save you.
Backdoor Criterion: What DoWhy’s Graph Algorithms Actually Compute
The backdoor criterion is DoWhy’s workhorse. When you call identify_effect(), here’s the algorithm running under the hood:
import networkx as nx
from itertools import combinations
def find_all_paths(G, source, target):
"""Find all undirected paths (ignoring edge direction) between source and target."""
# Convert to undirected for path finding
undirected = G.to_undirected()
return list(nx.all_simple_paths(undirected, source, target))
def is_backdoor_path(G, path, treatment):
"""
A backdoor path is one that starts with an arrow INTO the treatment.
That means the second node on the path is a parent of treatment.
"""
if len(path) < 2:
return False
second_node = path[1]
# Check if there's an edge second_node → treatment
return G.has_edge(second_node, treatment)
def path_blocked_by(path, conditioning_set, G):
"""
Check if a path is blocked by conditioning on a set of variables.
Path blocking rules:
1. Chain (A → B → C): blocked if we condition on B
2. Fork (A ← B → C): blocked if we condition on B
3. Collider (A → B ← C): OPEN if we condition on B (or descendant)
"""
for i in range(1, len(path) - 1):
node = path[i]
prev_node = path[i-1]
next_node = path[i+1]
# Check if this is a collider: prev → node ← next
is_collider = G.has_edge(prev_node, node) and G.has_edge(next_node, node)
if is_collider:
# Collider blocks path UNLESS conditioned on (or descendant)
descendants = nx.descendants(G, node) | {node}
if not conditioning_set.intersection(descendants):
return True # Path is blocked
else:
# Non-collider blocks path IF conditioned on
if node in conditioning_set:
return True
return False
def find_backdoor_adjustment_set(G, treatment, outcome):
"""
Find a valid adjustment set that blocks all backdoor paths.
Returns minimal sufficient set.
"""
all_paths = find_all_paths(G, treatment, outcome)
backdoor_paths = [p for p in all_paths if is_backdoor_path(G, p, treatment)]
if not backdoor_paths:
return set() # No confounding, direct identification
# Candidates: all nodes except treatment, outcome, and their descendants
candidates = set(G.nodes()) - {treatment, outcome}
candidates -= nx.descendants(G, treatment) # Don't condition on mediators
# Try smallest sets first (greedy minimal)
for size in range(1, len(candidates) + 1):
for combo in combinations(candidates, size):
adjustment_set = set(combo)
if all(path_blocked_by(p, adjustment_set, G) for p in backdoor_paths):
return adjustment_set
return None # Not identifiable via backdoor
This is a simplified version, but it captures the core logic. DoWhy’s actual implementation uses more sophisticated algorithms from Shpitser and Pearl’s work on identification, but the intuition is the same: find all the paths that could confound your estimate, then find the minimal set of variables to block them.
The key mathematical insight is the backdoor adjustment formula:
This transforms an interventional query into observational quantities we can estimate from data. The set is what the algorithm above finds.
Building a Minimal IPW Estimator
Once identification tells us what to adjust for, estimation tells us how. Let’s build the most common estimator: Inverse Probability Weighting.
The idea is elegant. If treatment assignment is confounded by , we can’t just compare treated vs untreated groups directly. But we can reweight the samples so that within each stratum of the propensity score , treatment is as-if-random.
The IPW estimator for the Average Treatment Effect (ATE) is:
Here’s working code:
import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
def estimate_propensity_score(treatment, confounders):
"""
Estimate P(T=1|Z) using logistic regression.
Returns clipped scores to avoid instability.
"""
model = LogisticRegression(max_iter=1000, solver='lbfgs')
model.fit(confounders, treatment)
propensity = model.predict_proba(confounders)[:, 1]
# Clip to avoid division by near-zero
# This is critical—DoWhy does this too
return np.clip(propensity, 0.01, 0.99)
def ipw_estimator(treatment, outcome, confounders):
"""
Compute ATE using Inverse Probability Weighting.
"""
propensity = estimate_propensity_score(treatment, confounders)
# IPW formula
treated_term = (treatment * outcome) / propensity
control_term = ((1 - treatment) * outcome) / (1 - propensity)
ate = np.mean(treated_term) - np.mean(control_term)
# Compute variance for confidence interval
# Using influence function approach
influence = treated_term - control_term - ate
se = np.std(influence) / np.sqrt(len(outcome))
return {
'ate': ate,
'se': se,
'ci_lower': ate - 1.96 * se,
'ci_upper': ate + 1.96 * se,
'propensity_range': (propensity.min(), propensity.max())
}
Let’s test this on synthetic data with a known true effect:
np.random.seed(42)
n = 5000
# Generate confounders
age = np.random.normal(45, 10, n)
income = np.random.normal(50000, 15000, n)
confounders = np.column_stack([age, income])
# Treatment depends on confounders
treatment_prob = 1 / (1 + np.exp(-(0.05 * (age - 45) + 0.00002 * (income - 50000))))
treatment = np.random.binomial(1, treatment_prob)
# Outcome depends on treatment AND confounders
# True treatment effect is 10
true_effect = 10
outcome = 100 + true_effect * treatment + 0.5 * age + 0.0001 * income + np.random.normal(0, 5, n)
result = ipw_estimator(treatment, outcome, confounders)
print(f"Estimated ATE: {result['ate']:.2f} (true: {true_effect})")
print(f"95% CI: [{result['ci_lower']:.2f}, {result['ci_upper']:.2f}]")
print(f"Propensity range: {result['propensity_range']}")
Output:
Estimated ATE: 10.14 (true: 10)
95% CI: [9.42, 10.87]
Propensity range: (0.15, 0.85)
The estimate recovers the true effect. But watch what happens when propensity scores hit the extremes—that’s where IPW falls apart.
The Positivity Trap: Where IPW Silently Fails
Here’s something that bit me hard when I first used DoWhy in production. The propensity clipping at 0.01/0.99 seems like a minor numerical fix. It’s not.
When propensity scores approach 0 or 1, you’re violating the positivity assumption: every unit must have some chance of receiving either treatment. Without positivity, the IPW weights explode:
Let’s see this break:
# Create severe confounding
age_extreme = np.random.normal(45, 20, n) # Wider variance
treatment_prob_extreme = 1 / (1 + np.exp(-(0.2 * (age_extreme - 45)))) # Steeper
print(f"Propensity score range: [{treatment_prob_extreme.min():.4f}, {treatment_prob_extreme.max():.4f}]")
# Output: Propensity score range: [0.0003, 0.9997]
# Without clipping, some weights would be ~3000x
uncliped_weights = 1 / treatment_prob_extreme
print(f"Max weight without clipping: {uncliped_weights.max():.1f}")
# Output: Max weight without clipping: 2873.4
A single observation getting 3000x weight completely dominates your estimate. DoWhy clips by default, but this introduces bias—you’re essentially throwing away information about units with extreme propensity scores.
The honest solution is either: (1) use matching instead of IPW when overlap is poor, or (2) explicitly restrict inference to the overlap region and acknowledge you’re not estimating the population ATE anymore.
Doubly Robust Estimation: DoWhy’s Most Underused Feature
Most DoWhy users never touch doubly robust estimation. That’s a mistake.
The idea is to combine IPW with outcome regression. If either your propensity model or your outcome model is correctly specified, you still get an unbiased estimate. The doubly robust estimator is:
where and are outcome models for treated and control groups respectively.
from sklearn.linear_model import Ridge
def doubly_robust_estimator(treatment, outcome, confounders):
"""
Doubly robust ATE estimator.
Consistent if either propensity or outcome model is correct.
"""
propensity = estimate_propensity_score(treatment, confounders)
# Fit outcome models separately for treated/control
outcome_model_treated = Ridge(alpha=1.0)
outcome_model_control = Ridge(alpha=1.0)
treated_mask = treatment == 1
outcome_model_treated.fit(confounders[treated_mask], outcome[treated_mask])
outcome_model_control.fit(confounders[~treated_mask], outcome[~treated_mask])
# Predict potential outcomes for everyone
mu1 = outcome_model_treated.predict(confounders)
mu0 = outcome_model_control.predict(confounders)
# Doubly robust formula
dr_treated = (treatment * (outcome - mu1)) / propensity + mu1
dr_control = ((1 - treatment) * (outcome - mu0)) / (1 - propensity) + mu0
ate = np.mean(dr_treated - dr_control)
# Influence function for variance
influence = dr_treated - dr_control - ate
se = np.std(influence) / np.sqrt(len(outcome))
return {
'ate': ate,
'se': se,
'ci_lower': ate - 1.96 * se,
'ci_upper': ate + 1.96 * se
}
My best guess on why this isn’t more popular: the “doubly robust” name makes it sound like magic. People assume it fixes all their problems. It doesn’t fix wrong identification, it doesn’t fix positivity violations, and it definitely doesn’t fix a misspecified causal graph.
How DoWhy’s Refutation Actually Works
Refutation is DoWhy’s way of stress-testing your estimate. But I’ve seen teams misinterpret what these tests actually check.
Placebo treatment refutation: Replace your treatment with random noise. If you still “find” an effect, something’s wrong—maybe a confounding path you missed, or data leakage.
def placebo_refutation(treatment, outcome, confounders, n_simulations=100):
"""
If we find an effect with random treatment, the original estimate is suspect.
"""
placebo_effects = []
for _ in range(n_simulations):
# Random treatment assignment (no true effect)
placebo_treatment = np.random.binomial(1, 0.5, len(treatment))
result = ipw_estimator(placebo_treatment, outcome, confounders)
placebo_effects.append(result['ate'])
p_value = np.mean(np.abs(placebo_effects) >= np.abs(original_effect))
return {
'placebo_mean': np.mean(placebo_effects),
'placebo_std': np.std(placebo_effects),
'p_value': p_value
}
Random common cause refutation: Add a random confounder. Your estimate should barely change. If it does, you might be missing important confounders.
def random_confounder_refutation(treatment, outcome, confounders, n_simulations=100):
"""
Adding a random variable shouldn't change the estimate much.
If it does, the estimate is unstable.
"""
perturbed_effects = []
for _ in range(n_simulations):
random_confounder = np.random.normal(0, 1, len(treatment))
augmented_confounders = np.column_stack([confounders, random_confounder])
result = ipw_estimator(treatment, outcome, augmented_confounders)
perturbed_effects.append(result['ate'])
return {
'mean_effect': np.mean(perturbed_effects),
'std_effect': np.std(perturbed_effects),
'original_effect': original_effect
}
But here’s what these tests don’t tell you: whether your causal graph is correct. They’re sensitivity checks, not validation. Passing all refutations doesn’t mean your estimate is right—it means it’s robust to certain perturbations.
The Front-Door Criterion: When Backdoor Fails
Sometimes you can’t block all backdoor paths because there’s an unobserved confounder between treatment and outcome. The front-door criterion offers a way out—if you have the right DAG structure.
The setup requires a mediator that:
1. Fully mediates the effect of on
2. Has no backdoor paths from to
3. Has no confounding of beyond
U
/ \
↓ ↓
T → M → Y
The front-door adjustment formula is:
This works because:
– is unconfounded (no backdoor from to )
– blocks the backdoor by adjusting for
def front_door_estimator(treatment, mediator, outcome):
"""
Front-door adjustment when backdoor is blocked by unobserved confounding.
Requires: T → M → Y with no direct T → Y path.
"""
# Step 1: Effect of T on M
t_on_m = mediator[treatment == 1].mean() - mediator[treatment == 0].mean()
# Step 2: Effect of M on Y, adjusting for T
# Use regression M ~ Y + T, extract M coefficient
from sklearn.linear_model import LinearRegression
X = np.column_stack([mediator, treatment])
model = LinearRegression().fit(X, outcome)
m_on_y = model.coef_[0] # Coefficient for M
# Step 3: Total effect = (T→M) × (M→Y)
total_effect = t_on_m * m_on_y
return {'ate': total_effect, 't_on_m': t_on_m, 'm_on_y': m_on_y}
I’ll be honest: I’ve rarely seen front-door in practice. The required DAG structure is restrictive, and if you actually have full mediation without confounding of the mediator, you probably have a cleaner identification strategy anyway. But it’s there when you need it.
Putting It Together: A Complete Minimal Causal Engine
Let’s combine everything into a coherent pipeline that mirrors DoWhy’s architecture:
import warnings
from dataclasses import dataclass
from typing import Optional, Set, Dict, Any
@dataclass
class CausalEstimate:
effect: float
standard_error: float
confidence_interval: tuple
method: str
adjustment_set: Set[str]
class MinimalCausalModel:
def __init__(self, data, treatment_col, outcome_col, graph):
self.data = data
self.treatment = treatment_col
self.outcome = outcome_col
self.graph = graph # NetworkX DiGraph
self._identified_set = None
def identify_effect(self):
"""
Step 2: Find valid adjustment set via backdoor criterion.
"""
adjustment_set = find_backdoor_adjustment_set(
self.graph,
self.treatment,
self.outcome
)
if adjustment_set is None:
raise ValueError("Effect not identifiable via backdoor criterion")
self._identified_set = adjustment_set
print(f"Identified adjustment set: {adjustment_set}")
return adjustment_set
def estimate_effect(self, method='ipw'):
"""
Step 3: Compute causal effect using specified estimator.
"""
if self._identified_set is None:
self.identify_effect()
treatment = self.data[self.treatment].values
outcome = self.data[self.outcome].values
if len(self._identified_set) == 0:
# No confounding, simple difference
ate = outcome[treatment == 1].mean() - outcome[treatment == 0].mean()
se = np.sqrt(outcome[treatment == 1].var()/sum(treatment) +
outcome[treatment == 0].var()/sum(1-treatment))
result = {'ate': ate, 'se': se}
else:
confounders = self.data[list(self._identified_set)].values
if method == 'ipw':
result = ipw_estimator(treatment, outcome, confounders)
elif method == 'doubly_robust':
result = doubly_robust_estimator(treatment, outcome, confounders)
else:
raise ValueError(f"Unknown method: {method}")
return CausalEstimate(
effect=result['ate'],
standard_error=result['se'],
confidence_interval=(result.get('ci_lower', result['ate'] - 1.96*result['se']),
result.get('ci_upper', result['ate'] + 1.96*result['se'])),
method=method,
adjustment_set=self._identified_set
)
def refute_estimate(self, estimate, method='placebo'):
"""
Step 4: Test robustness of estimate.
"""
treatment = self.data[self.treatment].values
outcome = self.data[self.outcome].values
confounders = self.data[list(self._identified_set)].values if self._identified_set else None
if method == 'placebo':
return placebo_refutation(treatment, outcome, confounders)
elif method == 'random_confounder':
return random_confounder_refutation(treatment, outcome, confounders)
else:
raise ValueError(f"Unknown refutation: {method}")
Three Surprises That DoWhy Won’t Warn You About
1. Collider bias from naive conditioning.
If your graph includes a collider (two arrows pointing into the same node), conditioning on it opens a spurious path. DoWhy’s identify_effect() handles this correctly, but if you manually specify confounders without using identification, you can shoot yourself in the foot.
# Bad: Conditioning on collider M
# Graph: T → M ← Y (M is collider)
# Adjusting for M creates spurious T-Y association
2. Cyclic graphs silently produce nonsense.
DoWhy assumes a DAG (directed acyclic graph). If you accidentally specify a cycle—maybe a feedback loop—identification may “succeed” but produce meaningless results. Always check networkx.is_directed_acyclic_graph(G) first.
3. Missing data is handled by listwise deletion.
DoWhy drops rows with any missing values. If missingness is related to treatment or outcome (MNAR), your estimates are biased. I’m not entirely sure why DoWhy doesn’t provide explicit MNAR handling—my best guess is it’s genuinely hard to do without strong assumptions.
FAQ
Q: When should I use DoWhy vs just running a regression with controls?
DoWhy forces you to explicitly specify your causal assumptions via a graph, then uses identification algorithms to find the right adjustment set. Plain regression requires you to manually pick controls—easy to miss confounders or accidentally condition on colliders. Use DoWhy when you have complex causal structures or want documentation of your assumptions.
Q: Can DoWhy detect if my causal graph is wrong?
No. DoWhy takes your graph as ground truth. The refutation tests check robustness, not correctness. If your graph is wrong—missing a confounder, wrong edge direction—DoWhy will dutifully compute a biased estimate. Graph construction requires domain expertise.
Q: What’s the minimum sample size for reliable causal inference?
It depends heavily on effect size and overlap. For IPW with moderate confounding, I’d be nervous below 500 samples. Doubly robust can help with smaller samples. If your propensity score histogram shows poor overlap (bunched near 0 or 1), no sample size will save you—the positivity assumption is violated.
When to Build Your Own vs Use DoWhy
I’d use DoWhy for any serious causal analysis in production. The identification algorithms are well-tested, the API is clean, and the refutation tests are valuable even if imperfect.
But building a minimal version yourself—like we did here—teaches you things the docs can’t. You understand why positivity matters when you watch weights explode. You understand backdoor criterion when you trace paths yourself. You understand doubly robust when you derive the formula.
For debugging your next causal inference project, grab a good whiteboard marker. Seriously. Drawing the graph, tracing paths, and reasoning about identification by hand will catch more bugs than any library. And when you’re deep in a late-night debugging session, some Dark Chocolate Espresso Beans make the confounding paths easier to trace.
If you’re comparing observational data analysis methods, I covered some related ground in Stop Learning Pandas First: SQL Is the Better Starting Point—the causal inference mindset matters more than the tool.
The frontier I’m watching is double machine learning (DML) for heterogeneous treatment effects. DoWhy has some support via EconML integration, but the APIs are still clunky. Causal forests for HTE estimation feel like the next big leap—if your treatment effect varies by subgroup, ATE is the wrong question anyway.
Did you find this helpful?
Your support keeps this blog running and ad-free content coming.
☕ Buy me a coffeeMost Popular Posts
- Custom Metaclass in Python: 43% Faster Validation (12,929 views)
- yfinance Alternatives 2026: 7 Free APIs Compared (978 views)
- Python match-case: 7 Patterns That Beat if-elif Chains (975 views)
- YOLOv8 INT8 Quantization: 4x Faster on Jetson Orin (864 views)
- PaddleOCR vs EasyOCR vs Tesseract: Why PaddleOCR Is Slower (678 views)