library(tidyverse); library(MatchIt); library(cobalt); library(estimatr); library(causalsens)
data("lalonde.psid", package = "causalsens")
data("lalonde.exp", package = "causalsens")
# Experimental benchmark from lalonde.exp (the gold standard we want to recover)
fit_exp <- # ___ YOUR CODE: lm_robust(re78 ~ treat, data = lalonde.exp, se_type = "HC1")
# Persist a CSV so the Python tab can read the same sample
write.csv(lalonde.psid, "lalonde.csv", row.names = FALSE)
covars <- c("age", "education", "black", "hispanic", "married", "nodegree",
"re74", "re75", "u74", "u75")
formula_match <- as.formula(paste("treat ~", paste(covars, collapse = " + ")))
# (a) Naive OLS — observational diff-in-means with HC1 SE
fit_naive <- # ___ YOUR CODE
# (b) Adjusted OLS — include all 10 covariates linearly, HC1 SE
fit_adj <- # ___ YOUR CODE
# (c) Exact matching on the four binary covariates only (method = "exact")
m_exact <- # ___ YOUR CODE: matchit(treat ~ black + hispanic + married + nodegree, ...)
# (d) Mahalanobis 1:1 matching with replacement
m_maha <- # ___ YOUR CODE: matchit(formula_match, method = "nearest",
# distance = "mahalanobis", replace = TRUE)
# (e) Coarsened Exact Matching with default coarsening
m_cem <- # ___ YOUR CODE: matchit(formula_match, method = "cem")
# (f) PSM — 1:1, logistic distance, caliper 0.2 * SD(logit PS), with replacement
m_psm <- # ___ YOUR CODE
# For each matched object, estimate the ATT as the weighted mean-difference on the
# matched sample (re78 ~ treat, weights = weights). For the SE use matched-set
# cluster-robust SEs — cluster by `subclass` (without replacement), or by `subclass`
# and the reused `id` via get_matches() (with replacement) — NOT plain HC1. Assemble
# all six estimates into a single results tibble with method, ATT, SE, and 95% CI.
# Bonus — Route 1 (Abadie-Imbens): also compute the AI analytic SE via
# Matching::Match(..., BiasAdjust = TRUE) for Mahalanobis and PSM, and compare
# the two routes. (CEM is not in the `Matching` package.)
# Produce three Love plots (Mahalanobis, CEM, PSM) via cobalt::love.plot,
# annotated with the headline result for each method.9 Session 6: Observed Confounding
9.1 Lecture Slides
9.2 Overview
This session is the foundation of Part III: Identification Under (Un-)Observed Confounding. We move from S04–S05’s experimental setting — where the analyst controls treatment assignment — to observational data, where treatment \(D\) is correlated with everything that predicts \(Y\). The pivot rests on a single assumption: conditional independence, \((Y(0), Y(1)) \perp D \mid \mathbf{X}\), given an adjustment set \(\mathbf{X}\) that closes every backdoor path.
The session walks through three estimator families that exploit this assumption — matching, inverse probability weighting (IPW), and augmented IPW (AIPW) — culminating in the doubly robust property: AIPW is consistent if either the outcome model or the propensity model is correctly specified. We then pivot to the question every honest observational analyst must answer: how strong would an unobserved confounder have to be to overturn the estimate? This motivates the modern toolkit of Rosenbaum bounds, E-values (VanderWeele and Ding 2017), and sensemakr (Cinelli and Hazlett 2020). The PBL completes the loop with a six-step workflow on the Lalonde NSW dataset and a controlled DGP where the truth is known.
Three forward-pointers: S07 (instrumental variables, when conditional independence fails), S08 (double machine learning, when \(\mathbf{X}\) is high-dimensional), S09 (treatment effect heterogeneity, when \(\tau\) varies with \(\mathbf{X}\)).
9.3 Lecture Notes
9.3.1 Identification under observed confounding
9.3.1.1 The conditional independence assumption — five names for one idea
The defining assumption of this session is
\[(Y(0), Y(1)) \perp D \mid \mathbf{X}\]
read as “potential outcomes are independent of treatment, given the adjustment set \(\mathbf{X}\)”. It is known by five names across disciplines: conditional independence (statistics), unconfoundedness (econometrics), ignorability (Bayesian inference), conditional exchangeability (epidemiology), and selection-on-observables (labour economics / sociology). All five describe the same condition: every variable that drives both \(D\) and \(Y\) is in \(\mathbf{X}\).
In S04, randomization made this hold unconditionally — the experimenter picked \(D\) uniformly at random, so \(D \perp (Y(0), Y(1))\) by construction. In S06, the conditional version places the burden on the analyst’s domain knowledge: the DAG that defines “the right \(\mathbf{X}\)” is a theoretical claim, not a data product.
9.3.1.2 Backdoor adjustment and the “good and bad controls” taxonomy
By S02’s backdoor criterion, an adjustment set \(\mathbf{X}\) is sufficient if it blocks every backdoor path from \(D\) to \(Y\) and contains no descendants of \(D\). Cinelli, Forney, and Pearl (Cinelli, Forney, and Pearl 2024) enumerate eight canonical configurations, sorted into three families:
- Good controls (M1–M3): classic confounders, indirect ancestors, and precision variables. Including them reduces bias or variance.
- Bad controls (M4–M6): post-treatment variables, colliders on \(D\) and \(Y\), and pre-treatment colliders of two unmeasured confounders (“M-bias”). Including them creates bias.
- Neutral / risky (M7–M8): pure instruments, which amplify any remaining unobserved-confounder bias, and pre-treatment variables of unknown structure, where the default is to include unless M6/M7 is suspected.
The pedagogical point: “throw everything at the model” is wrong. The analyst’s job is to draw a DAG, derive an adjustment set, and justify the variables included.
9.3.1.3 Positivity — every unit must have a chance of either treatment
The second identification assumption is positivity (also called overlap or common support):
\[0 < P(D = 1 \mid \mathbf{X} = \mathbf{x}) < 1 \quad \text{for every } \mathbf{x} \in \text{supp}(\mathbf{X}).\]
If positivity fails for some \(\mathbf{x}\), any estimator extrapolates beyond the data at that point. Two responses: trim the sample to the overlap region (e.g., propensity score \(\in [0.05, 0.95]\), Crump et al. 2009), or redefine the estimand via overlap weights (Li, Morgan, and Zaslavsky 2018). Both are now standard.
9.3.1.4 Identification — three roads to the same ATE
Under conditional independence and positivity, the ATE is identified from observed data in three algebraically equivalent ways:
- Adjustment formula (a.k.a. the g-formula): \(\tau = \mathbb{E}_X\!\left[\mathbb{E}(Y \mid D=1, X) - \mathbb{E}(Y \mid D=0, X)\right]\) — motivates the outcome-regression and matching families.
- IPW: \(\tau = \mathbb{E}\!\left[\dfrac{D \, Y}{e(X)} - \dfrac{(1-D) Y}{1 - e(X)}\right]\) — motivates inverse-probability weighting.
- AIPW: the two combined with a residual correction — yields a doubly robust estimator (AIPW).
In finite samples, the three expressions diverge — and the divergence is informative about which assumption is wrong.
9.3.2 Matching
9.3.2.1 The matching idea — build a synthetic RCT
Matching reads the adjustment formula literally: for each treated unit at \(\mathbf{X}_i\), find its nearest untreated twin under some distance metric, then average the paired differences. The target is the ATT, \(\tau_{\text{ATT}} = \mathbb{E}\!\left[\mathbb{E}(Y \mid D=1, X) - \mathbb{E}(Y \mid D=0, X) \mid D=1\right]\) — the conditional contrast averaged over the treated distribution of \(X\). The sample estimator (1:1, no replacement) is
\[\hat\tau_{\text{ATT}} = \frac{1}{n_1}\sum_{i: D_i = 1} \bigl(Y_i - Y_{m(i)}\bigr)\]
where \(m(i)\) is the index of the matched untreated unit. Three design knobs define the estimator: the distance (exact / Mahalanobis / propensity / CEM-coarsened), the matching ratio (1:1 / 1:M / variable-M with caliper), and replacement (yes / no). For the ATE rather than the ATT, also match each untreated unit to its nearest treated twin and combine by arm shares: \(\hat\tau_{\text{ATE}} = (n_1/n)\hat\tau_{\text{ATT}} + (n_0/n)\hat\tau_{\text{ATU}}\). The bias–variance direction of each knob: a caliper (drop matches with distance \(> c\); rule of thumb \(c = 0.2 \times\) SD of the logit propensity score) and 1:M matching trade bias against variance, and matching with replacement lowers bias at the cost of variance. The modern default is 1:M (\(M = 4\)) with replacement and a caliper of \(0.2\), reporting the number of unique controls used and the share of treated units matched.
9.3.2.2 Exact, Mahalanobis, and CEM
- Exact matching pairs \(i, j\) only if \(\mathbf{X}_i = \mathbf{X}_j\). Zero bias from the matching step; impossible in practice for continuous or high-dimensional \(\mathbf{X}\) (curse of dimensionality).
- Mahalanobis distance \(d_M(\mathbf{X}_i, \mathbf{X}_j) = \sqrt{(\mathbf{X}_i - \mathbf{X}_j)^\top \hat\Sigma_X^{-1} (\mathbf{X}_i - \mathbf{X}_j)}\) is scale-invariant and decorrelates. Pairs are not exact but are “balanced”. Default for continuous and ordinal \(\mathbf{X}\).
- Coarsened Exact Matching (CEM) (Iacus, King, and Porro 2012) bins each continuous covariate into a small number of strata, then does exact matching on the bins. Treated units in bins without untreated units (and vice versa) are pruned. The Monotonic Imbalance Bounding (MIB) property: bias is bounded before the data are touched, by the coarsening width — the analyst sets the maximum imbalance up front rather than diagnosing it ex post.
CEM is the modern recommendation for matching when sample size permits the pruning. Its price is sample loss: aggressive coarsening prunes more units; the estimand effectively becomes the ATT on the overlap region.
9.3.2.3 Propensity score matching and its discontents
The propensity score \(e(\mathbf{X}) = P(D = 1 \mid \mathbf{X})\) is a one-dimensional summary of the adjustment set. By the Rosenbaum–Rubin (1983) balancing-score theorem (Rosenbaum and Rubin 1983), conditioning on \(e(\mathbf{X})\) achieves the same identification as conditioning on the full vector \(\mathbf{X}\). Pre-2019 textbooks present propensity score matching (PSM) as the canonical matching method: estimate \(\hat e\) by logistic regression, match treated and untreated on \(\hat e\), estimate the ATT.
King & Nielsen (2019) (King and Nielsen 2019) published the modern critique: PSM can increase imbalance relative to no matching, because matching on \(\hat e\) approximates random sampling within strata of \(\hat e\) — analogous to a completely-randomized RCT rather than a blocked one. CEM and Mahalanobis approximate a fully blocked design and produce sample-by-sample covariate balance. The operational consensus today: prefer covariate-space matching (CEM or Mahalanobis) over PSM, choosing between them by checking both balance and retained-sample representativeness — CEM can over-prune (as the Ex 1 Lalonde lab shows, where it retains only 51 of 185 treated and lands far from the experimental benchmark), in which case Mahalanobis is preferred; reserve the propensity score for weighting (IPW) and AIPW (AIPW).
9.3.2.4 Balance diagnostics and the Love plot
After matching, the analyst must verify that the matched sample is balanced on every covariate in \(\mathbf{X}\). The standard metric is the standardized mean difference (SMD): \(\text{SMD}_k = (\bar X_{k,1} - \bar X_{k,0}) / \sqrt{(s_{k,1}^2 + s_{k,0}^2)/2}\) for covariate \(k\). Rule of thumb: \(|\text{SMD}_k| < 0.10\) means balance achieved (Stuart 2010). The Love plot visualizes raw vs. matched SMDs for every covariate side by side — cobalt::love.plot() in R.
Additional metrics for the same plot: variance ratio \(s_{k,1}^2 / s_{k,0}^2 \in [0.5, 2]\); K-S statistic for distribution-level checks; L1 distance for multidimensional histogram comparison (CEM’s native metric).
9.3.2.5 Inference after matching
The point estimate is mechanical (average paired difference). Standard errors are not:
- The naive t-test treats matched pairs as i.i.d. — wrong (matching introduces dependence through shared controls).
- The naive bootstrap of the full matching procedure is not generally consistent for matching SEs [Abadie & Imbens 2008].
- Abadie–Imbens (2006) (Abadie and Imbens 2006) analytic SEs account for the matching procedure and handle matching with replacement; implemented in the
Matchingpackage (Matching::Match()).MatchItitself returns no SEs. - Regression on the matched sample with
weights, plus matched-set cluster-robust SEs — not plain HC1, which treats the matched sample as i.i.d. and is inconsistent for matching SEs (errs in either direction) [Abadie & Spiess 2022]. Cluster bysubclasswithout replacement; bysubclassand the reused-unitid(viaget_matches()) with replacement. (Adding the covariates to the regression additionally regresses out residual imbalance — a related, not identical, alternative to the Abadie–Imbens estimator.)
The two routes target the same estimand and a careful analyst reports both as a robustness check: they agree on the qualitative story but differ in the exact finite-sample numbers, with the Abadie–Imbens variance the more conservative (it carries the full matching + conditional-variance uncertainty). The Ex 1 lab computes both for the with-replacement methods (Mahalanobis, PSM).
9.3.3 Inverse Probability Weighting (IPW)
9.3.3.1 The Horvitz–Thompson theorem
Under conditional independence and positivity, the ATE satisfies
\[\tau = \mathbb{E}\!\left[\frac{D \, Y}{e(\mathbf{X})} - \frac{(1 - D) \, Y}{1 - e(\mathbf{X})}\right]\]
— the Horvitz–Thompson (HT) estimator (Hirano, Imbens, and Ridder 2003). The proof is two lines of iterated expectation + conditional independence. The sample estimator plugs in \(\hat e(\mathbf{X}_i)\):
\[\hat\tau^{\text{HT}} = \frac{1}{n}\sum_{i=1}^n \left[\frac{D_i Y_i}{\hat e(\mathbf{X}_i)} - \frac{(1-D_i) Y_i}{1 - \hat e(\mathbf{X}_i)}\right].\]
9.3.3.2 Hájek normalization and diagnostics
The HT estimator is unbiased but high-variance. The Hájek (normalized) IPW divides by the sum of weights within each arm:
\[\hat\tau^{\text{Hájek}} = \frac{\sum_i D_i Y_i / \hat e_i}{\sum_i D_i / \hat e_i} - \frac{\sum_i (1-D_i) Y_i / (1-\hat e_i)}{\sum_i (1-D_i) / (1-\hat e_i)}.\]
Hájek is biased \(O(1/n)\) but lower-variance and bounded by \([\min Y, \max Y]\). It is the default in every IPW software package.
Diagnostics: report the effective sample size \(\text{ESS} = (\sum w_i)^2 / \sum w_i^2\), the maximum weight, and the propensity-score histogram for the two arms. Trim units with \(\hat e \notin [\alpha, 1-\alpha]\) (Crump et al. 2009) or use overlap weights \(\propto e(1-e)\) (Li, Morgan, Zaslavsky 2018) for a smoother alternative.
9.3.3.3 Modern balancing weights — enforce balance, don’t diagnose it
The classical workflow fits a propensity model, then checks balance ex post. The modern alternative enforces balance as a constraint on the weights:
- Entropy balancing (Hainmueller 2012): minimum-KL-divergence weights subject to exact moment-balance constraints.
- Covariate balancing propensity score (CBPS) (Imai and Ratkovic 2014): fit \(\hat e\) under a balancing objective rather than maximum likelihood.
- Stable balancing weights (SBW) (Zubizarreta 2015): minimum-variance weights subject to user-tolerated imbalance.
Athey & Imbens (2017) (Athey and Imbens 2017) recommend these as the modern default for applied causal inference.
9.3.4 Augmented IPW (AIPW)
9.3.4.1 The AIPW estimator
The two ingredients are the outcome regression (OR) \(\mu_d(\mathbf{x}) = \mathbb{E}(Y \mid D = d, \mathbf{X} = \mathbf{x})\) and the propensity score \(e(\mathbf{x})\). The AIPW estimator (Robins, Rotnitzky, and Zhao 1994) is
\[\hat\tau_{\text{AIPW}} = \frac{1}{n}\sum_{i=1}^n \!\left[\hat\mu_1(\mathbf{X}_i) - \hat\mu_0(\mathbf{X}_i) + \frac{D_i (Y_i - \hat\mu_1(\mathbf{X}_i))}{\hat e(\mathbf{X}_i)} - \frac{(1-D_i)(Y_i - \hat\mu_0(\mathbf{X}_i))}{1 - \hat e(\mathbf{X}_i)}\right].\]
Two readings of the same formula: (a) the OR difference, augmented by an IPW correction on the residuals; (b) the IPW estimator, with a correction that removes its dependence on \(\hat e\) where the OR model fits.
9.3.4.2 Double robustness
\(\hat\tau_{\text{AIPW}} \xrightarrow{p} \tau\) if either \(\hat\mu_d \xrightarrow{p} \mu_d\) (OR correct) or \(\hat e \xrightarrow{p} e\) (PS correct). The bias decomposes as the product of the two model errors — small if either factor is small. AIPW is also semiparametrically efficient (Hahn 1998) when both are correctly specified, and asymptotically normal: a standard error follows from the influence-function variance \(\widehat{\text{Var}}(\hat\tau_{\text{AIPW}}) = \frac{1}{n^2}\sum_i (\hat\psi_i - \hat\tau_{\text{AIPW}})^2\), where \(\hat\psi_i\) is the per-unit AIPW score (the bracketed term in the estimator), or equivalently from the bootstrap — this is what the Exercise 3 confidence intervals use.
The failure mode: if both models are wrong and not coincidentally offsetting, AIPW is biased. Exercise 2 demonstrates this on a controlled DGP (true ATE = 2.0): double robustness rescues the estimate when only the OR (scenario B) or only the PS (scenario C) is misspecified, and fails when both are wrong (scenario D).
9.3.4.3 Forward-pointer to S08 — Double Machine Learning
Classical AIPW uses parametric nuisances (logistic PS + linear OR). Three modern moves: replace each with ML estimators, use cross-fitting (split-sample fit/evaluate to remove overfit bias), and use Neyman-orthogonal score functions for valid \(\sqrt{n}\)-CIs. Together = Double Machine Learning [Chernozhukov et al. 2018], the topic of S08.
9.3.5 Sensitivity analysis
Conditional independence is untestable from the data: by definition, the unobserved confounder is absent from the dataset. Sensitivity analysis flips the question from “is the assumption right?” to “how strong would the violation have to be?”
9.3.5.1 Rosenbaum bounds
For matching estimators, Rosenbaum (2002) (Rosenbaum 2002) parametrizes hidden bias on the odds-ratio scale: within a matched pair \((i, j)\) with \(\mathbf{X}_i = \mathbf{X}_j\), the odds ratio of treatment due to an unobserved \(U\) lies in \([1/\Gamma, \Gamma]\). The sensitivity output is the largest \(\Gamma\) at which the matched-pair test (Wilcoxon signed-rank, say) still rejects. \(\Gamma = 1\) is no hidden bias; \(\Gamma \geq 2\) is strong robustness; \(\Gamma \geq 3\) is highly robust.
9.3.5.2 E-values
VanderWeele & Ding (2017) (VanderWeele and Ding 2017) define the E-value as the minimum risk-ratio strength that an unobserved \(U\) would need to have with both \(D\) and \(Y\) to fully explain away an observed effect:
\[E = \widehat{\text{RR}} + \sqrt{\widehat{\text{RR}}\bigl(\widehat{\text{RR}} - 1\bigr)} \quad \text{for } \widehat{\text{RR}} > 1.\]
For continuous outcomes, convert via \(\widehat{\text{RR}} \approx \exp(0.91 \cdot d)\) where \(d\) is the standardized mean difference. The E-value’s appeal: no model required, computable in one line, interpretable on the risk-ratio scale.
9.3.5.3 sensemakr
Cinelli & Hazlett (2020) (Cinelli and Hazlett 2020) build a sensitivity framework for linear regression (and AIPW with linear OR) parametrized by partial \(R^2\). Two summary statistics:
- Robustness Value (RV): the value of \(R^2_{Y \sim U \mid D, X} = R^2_{D \sim U \mid X}\) at which the point estimate would be reduced to a threshold (default: zero).
- RV-α: the value at which the \(t\)-statistic no longer rejects at level \(\alpha\).
The headline innovation is benchmarking: compare RV against the partial \(R^2\) of named observed covariates (“\(U\) would have to be at least \(k\) times stronger than age”). This makes the abstract sensitivity number actionable.
9.3.5.4 The reporting standard
For any observational estimate, report: (1) point estimate + CI; (2) robustness value; (3) at least one benchmark covariate with bound multiplier; (4) ovb_contour_plot; (5) one paragraph on the assumed lower bound on \(U\)’s strength and the operational conclusion.
9.3.6 The observed-confounding workflow
The six-step pipeline:
- Articulate the DAG — derive an adjustment set; flag good/bad controls.
- Check overlap — propensity histogram; trim or use overlap weights.
- Choose estimator — CEM / Mahalanobis / IPW / AIPW.
- Verify balance — Love plot + SMD table for every covariate.
- Stress-test — Rosenbaum bounds + E-value + sensemakr.
- Report — DAG + estimand + estimator + balance + RV.
Lock steps 1–3 before unblinding outcome data. This is the operational discipline that distinguishes a defensible analysis from a fishing expedition.
Why locking matters — the garden of forking paths. An observational analysis has many defensible choices (adjustment set, matching method, caliper width, trimming threshold, outcome model, subsample); with \(K\) binary choices there are \(2^K\) analytical paths, and the best-looking one will be biased (Gelman & Loken 2014). Pre-registration lists each decision and its value before \(Y\) is unblinded. A minimal pre-analysis plan for S06 work: (1) the DAG with the adjustment set highlighted; (2) the estimator and nuisance models with explicit functional form; (3) the trimming/caliper rule; (4) the inference method (sandwich / bootstrap); (5) the sensitivity benchmarks (which observed covariates). Pre-registration is now required in many medical and social-science journals, becoming standard in tech-company observational analyses, and supported for observational designs on OSF and EGAP (the AEA RCT Registry remains RCT-only).
9.4 Coding Lab
Four exercises walk through the observed-confounding pipeline end-to-end. Total budget ≈ 135 minutes:
| # | Topic | Slide stack | Time |
|---|---|---|---|
| 1 | Matching family on Lalonde (six estimators, three Love plots) | Matching | ~40 min |
| 2 | IPW: estimation, diagnostics, and balancing weights | IPW | ~40 min |
| 3 | AIPW and double robustness (four-scenario 2×2) | AIPW | ~30 min |
| 4 | Sensitivity analysis on the AIPW estimate | Sensitivity | ~25 min |
Exercise 1 runs six matching estimators side-by-side on Lalonde NSW, reproduces King–Nielsen (2019) in miniature, and asks which estimator recovers the experimental benchmark. Exercise 2 walks the full IPW workflow on a controlled DGP (known true ATE = 2.0): estimate (Horvitz–Thompson and Hájek), check balance on the weighted sample, then stress the DGP until positivity nearly fails and apply three remedies (trimming, overlap weights, entropy balancing) — contrasting fit-then-check IPW against enforce-by-design balancing weights. Exercise 3 reuses that DGP to demonstrate AIPW and the doubly-robust property across a four-scenario (OR × PS correctness) design, and shows how a misspecified propensity model is caught by the balance check. Exercise 4 runs sensitivity analysis on the Exercise 3 AIPW estimate after withholding the true confounder, using both E-values and sensemakr.
Each exercise has solutions in collapsed code blocks (both R and Python). Open the relevant panel-tabset, work on the exercise, then expand the solution to compare.
9.4.1 Exercise 1 — Matching family on Lalonde (~40 min)
The Lalonde National Supported Work (NSW) demonstration (1976–1977) is the canonical matching benchmark (Dehejia and Wahba 2002). The experimental treatment was a job-training program; the outcome is real earnings in 1978. We use two paired datasets shipped together in the causalsens R package:
lalonde.exp— the experimental sample (185 treated + 260 NSW controls). The unbiased ATT here is the gold-standard benchmark: \(\hat\tau_{\text{exp}} = \$1{,}794\) (HC1 SE \(\approx \$671\)).lalonde.psid— the observational variant. The 185 NSW treated units are paired with 2,490 full PSID controls — a sample chosen specifically because it does not match the treated population on covariates. PSID respondents are a representative slice of the US workforce; NSW participants are disadvantaged job-training enrollees. The naive earnings gap is therefore enormous and negative.
The full-PSID control sample is what makes the demonstration vivid: matching methods have $17,000 of selection bias to close, and the failure modes of each method become clearly visible.
Task: estimate the ATT on the observational lalonde.psid data using six methods, compute balance diagnostics on each, and reflect on which method converges to the experimental benchmark of $1,794.
9.4.1.1 Scaffold
import numpy as np, pandas as pd
import statsmodels.formula.api as smf
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from scipy.spatial import cKDTree
# Read the same lalonde.psid data the R tab wrote out
lalonde = pd.read_csv("lalonde.csv")
X_cols = ["age", "education", "black", "hispanic", "married", "nodegree",
"re74", "re75", "u74", "u75"]
# (a) Naive OLS
m_naive = # ___ YOUR CODE: smf.ols("re78 ~ treat", data=lalonde).fit(cov_type="HC1")
# (b) Adjusted OLS — include all 10 covariates linearly
m_adj = # ___ YOUR CODE
# (f) Hand-rolled 1:1 PSM:
# 1. StandardScaler on X (binary + dollar-scale features mix badly)
# 2. LogisticRegression(treat ~ X) → predict propensity ps
# 3. logit_ps = log(ps / (1 - ps))
# 4. scipy.spatial.cKDTree on untreated logit_ps; query for each treated unit's nearest neighbour
# 5. ATT = mean(Y[treated] - Y[matched_untreated])
# ___ YOUR CODE: implement steps 1–5; assemble a 3-row results DataFrame
# (Naive, Adjusted, hand-rolled PSM) with method, ATT, SE.
TipShow solution
Code
library(tidyverse)
library(MatchIt)
library(cobalt)
library(estimatr)
library(causalsens)
# Two paired Dehejia-Wahba datasets:
# - lalonde.psid : 185 NSW treated + 2,490 PSID controls (observational)
# - lalonde.exp : 185 NSW treated + 260 NSW experimental controls
data("lalonde.psid", package = "causalsens")
data("lalonde.exp", package = "causalsens")
# The experimental benchmark — what the matching methods are trying to recover.
fit_exp <- lm_robust(re78 ~ treat, data = lalonde.exp, se_type = "HC1")
cat(sprintf("Experimental benchmark ATT: $%0.1f (HC1 SE = $%0.1f)\n",
coef(fit_exp)["treat"], fit_exp$std.error["treat"]))Experimental benchmark ATT: $1794.3 (HC1 SE = $670.8)Code
# Persist a CSV copy so the Python tab can load the same sample without rpy2.
write.csv(lalonde.psid, "lalonde.csv", row.names = FALSE)
covars <- c("age", "education", "black", "hispanic", "married", "nodegree",
"re74", "re75", "u74", "u75")
formula_match <- as.formula(paste("treat ~", paste(covars, collapse = " + ")))
# (a) Naive OLS — treats observational sample as if randomized
fit_naive <- lm_robust(re78 ~ treat, data = lalonde.psid, se_type = "HC1")
# (b) Adjusted OLS — includes the covariates linearly
fit_adj <- lm_robust(reformulate(c("treat", covars), response = "re78"),
data = lalonde.psid, se_type = "HC1")
# (c) Exact matching on the four binary covariates only — `method = "exact"`
# uses every variable in the formula, so we pass a binary-only formula here.
m_exact <- matchit(treat ~ black + hispanic + married + nodegree,
data = lalonde.psid, method = "exact")
# (d) Mahalanobis distance matching, 1:1 with replacement
m_maha <- matchit(formula_match, data = lalonde.psid,
method = "nearest", distance = "mahalanobis", replace = TRUE)
# (e) Coarsened Exact Matching (CEM)
m_cem <- matchit(formula_match, data = lalonde.psid, method = "cem")
# (f) Propensity Score Matching, 1:1, logistic, caliper 0.2 * SD(logit PS)
m_psm <- matchit(formula_match, data = lalonde.psid,
method = "nearest", distance = "glm",
replace = TRUE, caliper = 0.2)
# Estimate the ATT from each matched sample as the weighted mean-difference
# (`re78 ~ treat`) — the matching does the balancing. Standard errors must
# respect the matched-set dependence (pairs share controls), so plain HC1 is
# invalid (Abadie & Spiess 2022). The valid cluster-robust SE differs by design:
# - WITHOUT replacement (exact, CEM): cluster by the matched set (`subclass`).
# - WITH replacement (Mahalanobis, PSM): a control can serve many treated units,
# so use `get_matches()` and cluster by BOTH the matched set and the reused
# unit `id` (two-way). The point estimate is identical either way.
estimate_att <- function(m, name, replace = FALSE) {
if (replace) {
gm <- get_matches(m)
fit <- lm(re78 ~ treat, data = gm, weights = weights)
V <- sandwich::vcovCL(fit, cluster = gm[c("subclass", "id")])
tibble(method = name, att = coef(fit)[["treat"]], se = sqrt(V["treat", "treat"]))
} else {
md <- match.data(m)
fit <- lm_robust(re78 ~ treat, data = md, weights = weights,
clusters = subclass, se_type = "CR2")
tibble(method = name, att = coef(fit)[["treat"]], se = fit$std.error[["treat"]])
}
}
results <- bind_rows(
tibble(method = "(a) Naive OLS", att = coef(fit_naive)["treat"], se = fit_naive$std.error["treat"]),
tibble(method = "(b) Adjusted OLS", att = coef(fit_adj)["treat"], se = fit_adj$std.error["treat"]),
estimate_att(m_exact, "(c) Exact", replace = FALSE),
estimate_att(m_maha, "(d) Mahalanobis", replace = TRUE),
estimate_att(m_cem, "(e) CEM", replace = FALSE),
estimate_att(m_psm, "(f) PSM", replace = TRUE)
) |>
mutate(ci_lo = att - 1.96 * se, ci_hi = att + 1.96 * se)
knitr::kable(results, digits = 0,
caption = "ATT estimates on `lalonde.psid` (NSW treated + full PSID controls). Experimental benchmark from `lalonde.exp` ≈ \\$1,794.")| method | att | se | ci_lo | ci_hi |
|---|---|---|---|---|
| (a) Naive OLS | -15205 | 656 | -16490 | -13919 |
| (b) Adjusted OLS | 4 | 854 | -1670 | 1678 |
| (c) Exact | -5605 | 1740 | -9016 | -2194 |
| (d) Mahalanobis | 1941 | 1162 | -337 | 4219 |
| (e) CEM | -3013 | 1451 | -5857 | -169 |
| (f) PSM | 970 | 1074 | -1135 | 3075 |
# Two routes to the SE, shown on the two WITH-replacement methods.
# Route 2 (the table above): a weighted regression on the matched sample plus a
# matched-set cluster-robust SE.
# Route 1: the Abadie-Imbens analytic variance via the `Matching` package, with
# the regression bias correction (BiasAdjust = TRUE).
# CEM is not implemented in `Matching`, so Route 1 here covers Mahalanobis and PSM.
# (Namespaced `Matching::Match` avoids attaching Matching/MASS in the shared session.)
Xm <- as.matrix(lalonde.psid[, covars])
ps <- glm(formula_match, data = lalonde.psid, family = binomial)$fitted
ai_maha <- Matching::Match(Y = lalonde.psid$re78, Tr = lalonde.psid$treat, X = Xm,
estimand = "ATT", Weight = 2, replace = TRUE, BiasAdjust = TRUE)
ai_psm <- Matching::Match(Y = lalonde.psid$re78, Tr = lalonde.psid$treat, X = ps,
estimand = "ATT", replace = TRUE, BiasAdjust = TRUE)
r2 <- results |> filter(method %in% c("(d) Mahalanobis", "(f) PSM"))
two_routes <- tibble(
method = r2$method,
`Route 2 — ATT` = r2$att,
`Route 2 — cluster SE` = r2$se,
`Route 1 — ATT (bias-adj)` = c(ai_maha$est, ai_psm$est),
`Route 1 — AI SE` = c(ai_maha$se, ai_psm$se)
)
knitr::kable(two_routes, digits = 0,
caption = "Two routes to inference (with-replacement methods). Route 2 = matched-sample regression + matched-set cluster-robust SE; Route 1 = Abadie-Imbens analytic SE via `Matching::Match(BiasAdjust = TRUE)`. Same estimand, different finite-sample numbers; the AI SE is the more conservative.")| method | Route 2 — ATT | Route 2 — cluster SE | Route 1 — ATT (bias-adj) | Route 1 — AI SE |
|---|---|---|---|---|
| (d) Mahalanobis | 1941 | 1162 | 2520 | 1705 |
| (f) PSM | 970 | 1074 | 762 | 1756 |
# Love plots: standardized mean differences before vs. after matching.
# Showing all three relevant methods so the reader can see which method
# achieves balance — and on which covariates.
love.plot(bal.tab(m_maha, un = TRUE),
stats = "mean.diffs",
thresholds = c(m = 0.1),
var.order = "unadjusted",
title = "Balance — Mahalanobis (winner: ATT $1,941, n_treated = 185)")Warning: Large mean differences detected; you may not be using standardized mean
differences for continuous variables.Warning: Standardized mean differences and raw mean differences are present in the same
plot. Use the `stars` argument to distinguish between them and appropriately
label the x-axis. See `love.plot()` for details.
love.plot(bal.tab(m_cem, un = TRUE),
stats = "mean.diffs",
thresholds = c(m = 0.1),
var.order = "unadjusted",
title = "Balance — CEM (over-pruned: only 51 of 185 treated retained)")Warning: Large mean differences detected; you may not be using standardized mean
differences for continuous variables.Warning: Standardized mean differences and raw mean differences are present in the same
plot. Use the `stars` argument to distinguish between them and appropriately
label the x-axis. See `love.plot()` for details.
love.plot(bal.tab(m_psm, un = TRUE),
stats = "mean.diffs",
thresholds = c(m = 0.1),
var.order = "unadjusted",
title = "Balance — PSM (King-Nielsen: residual imbalance on 6 of 10 covariates)")Warning: Large mean differences detected; you may not be using standardized mean
differences for continuous variables.Warning: Standardized mean differences and raw mean differences are present in the same
plot. Use the `stars` argument to distinguish between them and appropriately
label the x-axis. See `love.plot()` for details.
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
# Read the same lalonde.psid data the R tab wrote out (R writes CSV, Python reads CSV).
lalonde = pd.read_csv("lalonde.csv")
# (a) Naive OLS
m_naive = smf.ols("re78 ~ treat", data=lalonde).fit(cov_type="HC1")
# (b) Adjusted OLS
m_adj = smf.ols(
"re78 ~ treat + age + education + black + hispanic + married + nodegree"
" + re74 + re75 + u74 + u75",
data=lalonde,
).fit(cov_type="HC1")
# (c) Hand-rolled propensity score matching with scikit-learn + scipy.
# We estimate the propensity, do 1:1 nearest-neighbour matching on the logit
# of the propensity with replacement, then take the matched-sample mean
# difference as the ATT.
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from scipy.spatial import cKDTree
X_cols = ["age", "education", "black", "hispanic", "married", "nodegree",
"re74", "re75", "u74", "u75"]
X_raw = lalonde[X_cols].astype(float)
D = lalonde["treat"].values.astype(int)
Y = lalonde["re78"].values
# Standardize features before fitting (the dollar-scale re74/re75 are 4 orders
# of magnitude larger than the binary indicators; lbfgs would otherwise
# converge slowly or not at all).
X = StandardScaler().fit_transform(X_raw)
logit = LogisticRegression(max_iter=1000).fit(X, D)
ps = logit.predict_proba(X)[:, 1]
logit_ps = np.log(ps / (1 - ps))
# 1:1 nearest-neighbour matching on logit(PS), with replacement
treated_idx = np.where(D == 1)[0]
untreated_idx = np.where(D == 0)[0]
tree = cKDTree(logit_ps[untreated_idx].reshape(-1, 1))
dist, neighbor = tree.query(logit_ps[treated_idx].reshape(-1, 1), k=1)
matched_untreated = untreated_idx[neighbor]
att_psm_py = (Y[treated_idx] - Y[matched_untreated]).mean()
results = pd.DataFrame({
"method": ["(a) Naive OLS", "(b) Adjusted OLS", "(f) PSM (1:1, hand-rolled)"],
"att": [m_naive.params["treat"], m_adj.params["treat"], att_psm_py],
"se": [m_naive.bse["treat"], m_adj.bse["treat"], np.nan],
})
print(results.round(0).to_string(index=False)) method att se
(a) Naive OLS -15205.0 656.0
(b) Adjusted OLS 4.0 854.0
(f) PSM (1:1, hand-rolled) 2318.0 NaNprint("\nExperimental benchmark (lalonde.exp) ≈ $1,794")
Experimental benchmark (lalonde.exp) ≈ $1,794The Python lab implements (a) Naive OLS, (b) Adjusted OLS, and (f) a hand-rolled 1:1 PSM via sklearn + scipy.spatial. For CEM and Mahalanobis, the cleanest production tool is R’s MatchIt; the R tab above is the reference implementation. Note on the hand-rolled PSM: it is deliberately simplified — it standardizes the covariates, fits an L2-regularized sklearn logistic propensity, and matches every treated unit to its nearest control on logit(\(\hat e\)) with no caliper — so it retains all 185 treated and lands near $2,318, not R’s $970 (which uses an unregularized glm distance and a 0.2-SD caliper that drops poorly-matched treated units). The two illustrate the same idea, not a like-for-like replication, and the gap is itself the lesson: PSM is highly implementation-dependent. Bonus exercise: extend the Python code with a Mahalanobis matcher via scipy.spatial.cKDTree and the inverse covariance matrix. Inference note: the Route-1 Abadie–Imbens analytic SEs (R’s Matching package) have no off-the-shelf Python equivalent — the Python tab reports point estimates, with matching-aware inference deferred to the R reference implementation above.
NoteSolution sketch — Exercise 1
The six estimators on lalonde.psid (185 NSW treated + 2,490 PSID controls) tell a story in three acts.
Act 1 — selection bias is enormous. Naive OLS gives −$15,205 (HC1 SE = $656). PSID respondents earn an average of ~$21,500 in 1978 while NSW treated earn ~$6,300; the unconditional comparison loads the entire ~$17k gap onto “treatment”. The selection-on-observables assumption is wildly violated when we ignore covariates.
Act 2 — pure linear adjustment fails. Adjusted OLS lands at $4 — essentially zero, $1,790 below the experimental benchmark. The linear functional form simply cannot extrapolate across two populations whose covariate distributions barely overlap; OLS fits a global hyperplane that is dominated by the 2,490-row PSID region and learns a covariate-treatment relationship from sample units the treated group does not resemble.
Act 3 — matching unlocks the comparison. Now the methods diverge:
- (d) Mahalanobis matching wins decisively: ATT = $1,941 (matched-set cluster-robust SE = $1,162), $147 from the $1,794 experimental benchmark and within its 95% CI. All 185 treated units kept; max post-match |SMD| = 0.25 across 10 covariates.
- (f) PSM is biased low at $970, with the worst remaining balance (6 of 10 covariates with |SMD|>0.10, max 0.36). The propensity-score collapse onto a single dimension allows opposite-direction imbalances to cancel inside \(\hat e(\mathbf{X})\) — exactly the King–Nielsen (2019) mechanism.
- (e) CEM is the most surprising result: ATT = −$3,013 (negative!), with only 51 of 185 treated units surviving the default coarsening. The remaining matched sample is small, badly balanced (max |SMD| = 0.63), and unrepresentative — over-aggressive coarsening prunes the treated units that lack PSID neighbours, leaving a non-random sliver. CEM’s monotonic-imbalance-bounding property is about the coarsening grid, not about the retained sample’s representativeness; the analyst must inspect both.
- (c) Exact (binary covariates only) lands at −$5,605 — exact matching on the four binaries leaves continuous covariates (age, education, re74, re75) wildly unbalanced. Pure exact matching is rarely useful for continuous covariates.
The headline: Mahalanobis is the workhorse for this sample. CEM with default settings over-prunes — the analyst must reduce coarsening or trade off against PSM-style “everyone retained but residual imbalance”. PSM’s residual imbalance and biased point estimate illustrate exactly the King–Nielsen critique. The Adjusted-OLS = $4 result is a sharp pedagogical reminder that “include the covariates linearly” is not a substitute for explicit balance.
Act 4 — two routes to the SE. The headline table uses Route 2 (a weighted regression on the matched sample with matched-set cluster-robust SEs). Route 1 — the Abadie–Imbens analytic variance via Matching::Match(BiasAdjust = TRUE) — tells the same story with different exact numbers: Mahalanobis moves from $1,941 (cluster SE $1,162) to a bias-corrected $2,520 (AI SE $1,705), and PSM from $970 ($1,074) to $762 ($1,756). The two routes share a probability limit and agree on the verdict (Mahalanobis ≈ the benchmark, PSM lower), but the AI SE is the more conservative — it carries the full matching and conditional-variance uncertainty. Report both as a robustness check; never fall back to plain HC1 on the matched sample (Abadie and Spiess 2022).
9.4.2 Exercise 2 — IPW: estimation, diagnostics, and balancing weights (~40 min)
The propensity-score weighting workflow has three moves, in order: estimate the effect (Horvitz–Thompson and Hájek IPW), check balance on the weighted sample (did the weights actually equalize the covariates?), and — when standard IPW turns out fragile — reach for the modern balancing-weight alternatives. We run all three on a controlled DGP with a known true ATE = 2.0, then stress the DGP until positivity nearly fails and compare three remedies.
9.4.2.1 Scaffold
set.seed(20260601L)
n <- 5000
# Controlled DGP: two true confounders (X1, X2), noise covariate (X3),
# post-treatment M (used in Ex 3). True ATE = 2.0.
df <- tibble(
X1 = rnorm(n), X2 = rnorm(n), X3 = rnorm(n),
D = rbinom(n, 1, plogis(-0.5 + 0.8 * X1 - 0.6 * X2)),
Y = 1.0 + 2.0 * D + 1.5 * X1 - 1.0 * X2 + 0.5 * D * X1 + rnorm(n, sd = 1.2),
M = D + 0.3 * X1 + rnorm(n, sd = 0.3)
)
# (1) ESTIMATE — fit the correct PS, then Horvitz-Thompson and Hájek IPW.
e_hat <- # ___ YOUR CODE: predict(glm(D ~ X1 + X2, df, family = binomial()), type="response")
ht <- # ___ YOUR CODE: mean(D*Y/e_hat) - mean((1-D)*Y/(1-e_hat))
hajek <- # ___ YOUR CODE: self-normalized Hájek form
# (2) CHECK BALANCE — did the weights equalize X1, X2? The weighted Love plot.
# w <- WeightIt::weightit(D ~ X1 + X2 + X3, data = df, estimand = "ATE")
# cobalt::love.plot(w, stats = "mean.diffs", abs = TRUE, thresholds = c(m = 0.1))
# (3) DIAGNOSE under stress — rebuild with steep selection so positivity nearly fails.
# df_stress: D ~ plogis(-0.5 + 2.4*X1 - 1.8*X2). Compute ESS/n, max(w),
# out-of-[0.05,0.95] count; plot the mirrored PS histogram.
# (4) REMEDIES — standard Hájek, trim, overlap weights (ATO), entropy balancing.
# Then a weighted Love plot contrasting raw vs IPW vs entropy-balanced.import numpy as np, pandas as pd
import statsmodels.api as sm, statsmodels.formula.api as smf
from scipy.optimize import minimize
from scipy.special import logsumexp
rng = np.random.default_rng(20260601)
n = 5000
X1, X2, X3 = rng.standard_normal(n), rng.standard_normal(n), rng.standard_normal(n)
D = rng.binomial(1, 1 / (1 + np.exp(-(-0.5 + 0.8 * X1 - 0.6 * X2))))
Y = (1.0 + 2.0 * D + 1.5 * X1 - 1.0 * X2 + 0.5 * D * X1
+ rng.standard_normal(n) * 1.2)
M = D + 0.3 * X1 + rng.standard_normal(n) * 0.3
df = pd.DataFrame({"X1": X1, "X2": X2, "X3": X3, "D": D, "Y": Y, "M": M})
# (1) ESTIMATE — HT and Hájek IPW with the correct PS.
# (2) CHECK BALANCE — weighted SMD of X1, X2, X3 (treated vs control) under IPW weights.
# (3) DIAGNOSE under stress — ESS/n, max(w), out-of-range count.
# (4) REMEDIES — Hájek, trim, overlap weights, hand-rolled entropy balancing; then
# compare weighted SMDs raw vs IPW vs entropy-balanced.
TipShow solution
Code
set.seed(20260601L)
n <- 5000
# Two true confounders (X1, X2), one noise covariate (X3), one post-treatment
# variable (M, used in Ex 3 — students should NOT include it in any adjustment set).
# True ATE = 2.0.
df <- tibble(
X1 = rnorm(n),
X2 = rnorm(n),
X3 = rnorm(n),
D = rbinom(n, 1, plogis(-0.5 + 0.8 * X1 - 0.6 * X2)),
Y = 1.0 + 2.0 * D + 1.5 * X1 - 1.0 * X2 + 0.5 * D * X1 + rnorm(n, sd = 1.2),
M = D + 0.3 * X1 + rnorm(n, sd = 0.3)
)
# True ATE:
true_ate <- mean((1.0 + 2.0 * 1 + 1.5 * df$X1 - 1.0 * df$X2 + 0.5 * 1 * df$X1) -
(1.0 + 2.0 * 0 + 1.5 * df$X1 - 1.0 * df$X2 + 0.5 * 0 * df$X1))
cat("True ATE (population): ", round(true_ate, 3), "\n")True ATE (population): 2.008 # (1) ESTIMATE — fit the correct PS, then Horvitz-Thompson and Hájek IPW.
e_hat <- predict(glm(D ~ X1 + X2, data = df, family = binomial()), type = "response")
ht <- with(df, mean(D * Y / e_hat) - mean((1 - D) * Y / (1 - e_hat)))
hajek <- with(df, sum(D * Y / e_hat) / sum(D / e_hat) -
sum((1 - D) * Y / (1 - e_hat)) / sum((1 - D) / (1 - e_hat)))
cat(sprintf("Horvitz-Thompson IPW: %.3f | Hájek IPW: %.3f (true ATE = 2.0)\n",
ht, hajek))Horvitz-Thompson IPW: 2.036 | Hájek IPW: 2.084 (true ATE = 2.0)# (2) CHECK BALANCE — did the weights equalize the covariates? The weighted Love plot.
library(WeightIt); library(cobalt)
w_ate <- weightit(D ~ X1 + X2 + X3, data = df, method = "glm", estimand = "ATE")
love.plot(w_ate, stats = "mean.diffs", abs = TRUE, thresholds = c(m = 0.1),
colors = c("#FF7E15", "#005e73"), drop.distance = TRUE,
sample.names = c("Unweighted", "IPW-weighted"),
title = "IPW balance check — clean DGP")
Code
# (3) DIAGNOSE — stress the DGP: bigger treatment-selection coefficients on
# X1, X2 push propensities toward 0 and 1.
set.seed(20260601L)
n <- 5000
df_stress <- tibble(
X1 = rnorm(n),
X2 = rnorm(n),
D = rbinom(n, 1, plogis(-0.5 + 2.4 * X1 - 1.8 * X2)), # was 0.8, -0.6
Y = 1.0 + 2.0 * D + 1.5 * X1 - 1.0 * X2 + 0.5 * D * X1 + rnorm(n, sd = 1.2)
)
# Correctly-specified PS — the issue here is not misspecification, it's overlap.
ps_fit <- glm(D ~ X1 + X2, data = df_stress, family = binomial())
e_hat <- predict(ps_fit, type = "response")
w_std <- with(df_stress, ifelse(D == 1, 1 / e_hat, 1 / (1 - e_hat)))
ess_pct <- (sum(w_std))^2 / sum(w_std^2) / n * 100
cat(sprintf("Standard IPW diagnostics on the stressed DGP (true ATE = 2.0):\n"))Standard IPW diagnostics on the stressed DGP (true ATE = 2.0):Code
cat(sprintf(" ESS / n = %.1f%% (rule of thumb: < 50%% is heavy concentration)\n", ess_pct)) ESS / n = 14.8% (rule of thumb: < 50% is heavy concentration)Code
cat(sprintf(" max(w) = %.1f (rule of thumb: > 10 deserves investigation)\n", max(w_std))) max(w) = 172.1 (rule of thumb: > 10 deserves investigation)Code
cat(sprintf(" units with e_hat < 0.05 or > 0.95: %d (out of %d)\n",
sum(e_hat < 0.05 | e_hat > 0.95), n)) units with e_hat < 0.05 or > 0.95: 1643 (out of 5000)# Mirrored PS histogram — the visual that goes with the numbers above.
ymax <- max(table(cut(e_hat, breaks = seq(0, 1, length.out = 41))))
ggplot() +
annotate("rect", xmin = 0.05, xmax = 0.95,
ymin = -ymax * 1.1, ymax = ymax * 1.1,
fill = "grey70", alpha = 0.10) +
geom_histogram(data = tibble(e_hat = e_hat[df_stress$D == 1]),
aes(x = e_hat, y = after_stat(count)),
fill = "#FF7E15", color = NA, alpha = 0.85,
breaks = seq(0, 1, length.out = 41)) +
geom_histogram(data = tibble(e_hat = e_hat[df_stress$D == 0]),
aes(x = e_hat, y = -after_stat(count)),
fill = "#005e73", color = NA, alpha = 0.85,
breaks = seq(0, 1, length.out = 41)) +
geom_vline(xintercept = c(0.05, 0.95), linetype = "dashed",
color = "grey40", linewidth = 0.4) +
geom_hline(yintercept = 0, color = "white", linewidth = 0.6) +
scale_y_continuous(labels = abs) +
labs(x = expression(hat(e)(X)),
y = "count (D=1 above, D=0 below)",
title = "Stressed DGP — propensities pile up at the boundaries") +
theme_minimal(base_size = 12)
# (4) REMEDIES for marginal positivity.
# (i) Standard IPW (Hájek)
hj_std <- with(df_stress,
sum(D * Y / e_hat) / sum(D / e_hat) -
sum((1 - D) * Y / (1 - e_hat)) / sum((1 - D) / (1 - e_hat))
)
# (ii) Trim at [0.05, 0.95] and re-estimate Hájek IPW on the overlap region
keep <- e_hat > 0.05 & e_hat < 0.95
d_t <- df_stress[keep, ]; e_t <- e_hat[keep]
hj_trim <- with(d_t,
sum(D * Y / e_t) / sum(D / e_t) -
sum((1 - D) * Y / (1 - e_t)) / sum((1 - D) / (1 - e_t))
)
# (iii) Overlap weights (Li, Morgan, Zaslavsky 2018) — estimand: ATO
w_ow <- e_hat * (1 - e_hat)
ato_ow <- with(df_stress,
sum(D * Y * w_ow / e_hat) / sum(D * w_ow / e_hat) -
sum((1 - D) * Y * w_ow / (1 - e_hat)) / sum((1 - D) * w_ow / (1 - e_hat))
)
# (iv) Entropy balancing (Hainmueller 2012) — exact moment balance, no PS fit
ebw <- weightit(D ~ X1 + X2, data = df_stress, method = "ebal", estimand = "ATE")
fit_ebw <- lm_robust(Y ~ D, data = df_stress, weights = ebw$weights, se_type = "HC1")
ate_ebw <- coef(fit_ebw)["D"]
tibble(
estimator = c("Hájek IPW (standard)",
"Hájek IPW (trim at [0.05, 0.95])",
"Overlap-weighted (Li-Morgan-Zaslavsky)",
"Entropy balancing (Hainmueller)"),
estimand = c("ATE", "ATE on overlap region", "ATO (equipoise pop.)", "ATE"),
estimate = round(c(hj_std, hj_trim, ato_ow, ate_ebw), 3)
) |>
knitr::kable(caption = "Three remedies for marginal positivity. True ATE = 2.0; the ATO need not equal the ATE but should land in the same neighborhood under conditional ignorability.")| estimator | estimand | estimate |
|---|---|---|
| Hájek IPW (standard) | ATE | 2.289 |
| Hájek IPW (trim at [0.05, 0.95]) | ATE on overlap region | 2.043 |
| Overlap-weighted (Li-Morgan-Zaslavsky) | ATO (equipoise pop.) | 2.065 |
| Entropy balancing (Hainmueller) | ATE | 2.028 |
# (5) BALANCE under stress — fit-then-check (IPW) vs enforce-by-design (entropy balancing).
w_ipw_s <- weightit(D ~ X1 + X2, data = df_stress, method = "glm", estimand = "ATE")
w_eb_s <- weightit(D ~ X1 + X2, data = df_stress, method = "ebal", estimand = "ATE")
love.plot(D ~ X1 + X2, data = df_stress,
weights = list(IPW = w_ipw_s$weights, `Entropy bal.` = w_eb_s$weights),
stats = "mean.diffs", abs = TRUE, thresholds = c(m = 0.1),
colors = c("grey65", "#FF7E15", "#005e73"),
sample.names = c("Unweighted", "IPW", "Entropy bal."),
title = "Stressed DGP: IPW vs entropy balancing")
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from scipy.optimize import minimize
from scipy.special import logsumexp
rng = np.random.default_rng(20260601)
n = 5000
X1 = rng.standard_normal(n)
X2 = rng.standard_normal(n)
X3 = rng.standard_normal(n)
D = rng.binomial(1, 1 / (1 + np.exp(-(-0.5 + 0.8 * X1 - 0.6 * X2))))
Y = (1.0 + 2.0 * D + 1.5 * X1 - 1.0 * X2 + 0.5 * D * X1
+ rng.standard_normal(n) * 1.2)
M = D + 0.3 * X1 + rng.standard_normal(n) * 0.3
df = pd.DataFrame({"X1": X1, "X2": X2, "X3": X3, "D": D, "Y": Y, "M": M})
true_ate = float(((1.0 + 2.0 * 1 + 1.5 * X1 - 1.0 * X2 + 0.5 * 1 * X1)
- (1.0 + 2.0 * 0 + 1.5 * X1 - 1.0 * X2 + 0.5 * 0 * X1)).mean())
print(f"True ATE (population): {true_ate:.3f}")True ATE (population): 2.006# (1) ESTIMATE — HT and Hájek IPW with the correct PS.
e_hat = smf.glm("D ~ X1 + X2", data=df, family=sm.families.Binomial()).fit().predict(df).values
ht = np.mean(D * Y / e_hat) - np.mean((1 - D) * Y / (1 - e_hat))
hj = (np.sum(D * Y / e_hat) / np.sum(D / e_hat)
- np.sum((1 - D) * Y / (1 - e_hat)) / np.sum((1 - D) / (1 - e_hat)))
print(f"Horvitz-Thompson IPW: {ht:.3f} | Hajek IPW: {hj:.3f} (true ATE = 2.0)")Horvitz-Thompson IPW: 2.053 | Hajek IPW: 2.038 (true ATE = 2.0)# (2) CHECK BALANCE — weighted standardized mean differences (treated vs control).
# Denominator = unweighted pooled SD (matches cobalt's default), so unadjusted and
# weighted SMDs are comparable.
def smd_w(x, D, w):
mt = np.average(x[D == 1], weights=w[D == 1])
mc = np.average(x[D == 0], weights=w[D == 0])
s = np.sqrt((np.var(x[D == 1], ddof=1) + np.var(x[D == 0], ddof=1)) / 2)
return (mt - mc) / s
w_ate = np.where(D == 1, 1 / e_hat, 1 / (1 - e_hat))
print("Balance check — IPW weights on the clean DGP (|SMD| < 0.1 = balanced):")Balance check — IPW weights on the clean DGP (|SMD| < 0.1 = balanced):print(pd.DataFrame({
"covariate": ["X1", "X2", "X3"],
"SMD_unweighted": [smd_w(v, D, np.ones(n)) for v in (X1, X2, X3)],
"SMD_IPW": [smd_w(v, D, w_ate) for v in (X1, X2, X3)],
}).round(3).to_string(index=False))covariate SMD_unweighted SMD_IPW
X1 0.719 -0.023
X2 -0.471 0.005
X3 0.009 0.003# (3) DIAGNOSE — stress the DGP.
rng = np.random.default_rng(20260601)
n = 5000
X1 = rng.standard_normal(n)
X2 = rng.standard_normal(n)
D = rng.binomial(1, 1 / (1 + np.exp(-(-0.5 + 2.4 * X1 - 1.8 * X2))))
Y = (1.0 + 2.0 * D + 1.5 * X1 - 1.0 * X2 + 0.5 * D * X1
+ rng.standard_normal(n) * 1.2)
df_stress = pd.DataFrame({"X1": X1, "X2": X2, "D": D, "Y": Y})
ps_fit = smf.glm("D ~ X1 + X2", data=df_stress,
family=sm.families.Binomial()).fit()
e_hat = ps_fit.predict(df_stress).values
w_std = np.where(D == 1, 1 / e_hat, 1 / (1 - e_hat))
ess_pct = (w_std.sum() ** 2) / (w_std ** 2).sum() / n * 100
print(f"ESS / n = {ess_pct:.1f}%")ESS / n = 22.1%print(f"max(w) = {w_std.max():.1f}")max(w) = 95.9print(f"units with e_hat < 0.05 or > 0.95: "
f"{int(((e_hat < 0.05) | (e_hat > 0.95)).sum())} of {n}")units with e_hat < 0.05 or > 0.95: 1679 of 5000# (4) REMEDIES
def hajek(D, Y, e):
return (np.sum(D * Y / e) / np.sum(D / e)
- np.sum((1 - D) * Y / (1 - e)) / np.sum((1 - D) / (1 - e)))
hj_std = hajek(D, Y, e_hat)
keep = (e_hat > 0.05) & (e_hat < 0.95)
hj_trim = hajek(D[keep], Y[keep], e_hat[keep])
w_ow = e_hat * (1 - e_hat)
ato_ow = (np.sum(D * Y * w_ow / e_hat) / np.sum(D * w_ow / e_hat)
- np.sum((1 - D) * Y * w_ow / (1 - e_hat))
/ np.sum((1 - D) * w_ow / (1 - e_hat)))
# Hand-rolled entropy balancing via the dual (see Hainmueller 2012). CENTER the
# covariates on the target moments, then minimize logsumexp(-(X_c @ lambda)); at
# the optimum the weights satisfy E_w[X_c] = mu_target by construction.
def ebal_weights(X_c, mu_target):
Xc = X_c - mu_target
def loss(lam):
return logsumexp(-(Xc @ lam))
def grad(lam):
w = np.exp(-(Xc @ lam) - logsumexp(-(Xc @ lam)))
return -(w[:, None] * Xc).sum(axis=0)
res = minimize(loss, x0=np.zeros(Xc.shape[1]), jac=grad, method="BFGS")
w = np.exp(-(Xc @ res.x) - logsumexp(-(Xc @ res.x)))
return w / w.sum()
X_all = np.column_stack([X1, X2])
mu_all = X_all.mean(axis=0)
w_t = ebal_weights(X_all[D == 1], mu_all)
w_c = ebal_weights(X_all[D == 0], mu_all)
ate_ebw = (Y[D == 1] @ w_t) - (Y[D == 0] @ w_c)
print(pd.DataFrame({
"estimator": ["Hájek IPW (standard)",
"Hájek IPW (trim at [0.05, 0.95])",
"Overlap-weighted",
"Entropy balancing"],
"estimand": ["ATE", "ATE on overlap", "ATO", "ATE"],
"estimate": np.round([hj_std, hj_trim, ato_ow, ate_ebw], 3),
}).to_string(index=False)) estimator estimand estimate
Hájek IPW (standard) ATE 2.564
Hájek IPW (trim at [0.05, 0.95]) ATE on overlap 2.042
Overlap-weighted ATO 2.069
Entropy balancing ATE 2.019# (5) BALANCE under stress — IPW (approximate) vs entropy balancing (exact).
w_eb = np.empty(n); w_eb[D == 1] = w_t; w_eb[D == 0] = w_c
print("\nWeighted |SMD| under stress (denominator = unweighted pooled SD):")
Weighted |SMD| under stress (denominator = unweighted pooled SD):print(pd.DataFrame({
"covariate": ["X1", "X2"],
"unweighted": [abs(smd_w(v, D, np.ones(n))) for v in (X1, X2)],
"IPW": [abs(smd_w(v, D, w_std)) for v in (X1, X2)],
"ebal": [abs(smd_w(v, D, w_eb)) for v in (X1, X2)],
}).round(4).to_string(index=False))covariate unweighted IPW ebal
X1 1.3046 0.2617 0.0
X2 0.9209 0.1451 0.09.4.3 Exercise 3 — AIPW and double robustness (~30 min)
We reuse the controlled DGP from Exercise 2 (true ATE = 2.0) to see the doubly-robust property cleanly: AIPW augments the IPW estimator with an outcome-regression correction, and is consistent if either the propensity model or the outcome model is right. We compare three estimators (IPW, OR-alone, AIPW) across four DGP-specified scenarios for model correctness.
9.4.3.1 Scaffold
# Reuses `df` from Exercise 2. Helper: fit PS via glm, fit two outcome
# regressions (D=1, D=0), and return Horvitz-Thompson, Hájek, and AIPW.
ipw_aipw <- function(df, ps_form, or_form) {
ps_fit <- # ___ YOUR CODE: glm(ps_form, data = df, family = binomial())
e_hat <- # ___ YOUR CODE: predict(ps_fit, type = "response")
keep <- # ___ YOUR CODE: trim e_hat to [0.02, 0.98]
d <- df[keep, ]; e_hat <- e_hat[keep]; D <- d$D; Y <- d$Y
ht <- # ___ YOUR CODE
hj <- # ___ YOUR CODE: self-normalized Hájek form
or1 <- # ___ YOUR CODE: lm(or_form, data = d[d$D == 1, ])
or0 <- # ___ YOUR CODE
mu1 <- predict(or1, newdata = d); mu0 <- predict(or0, newdata = d)
aipw <- # ___ YOUR CODE: mean(mu1 - mu0 + D*(Y-mu1)/e_hat - (1-D)*(Y-mu0)/(1-e_hat))
tibble(HT = ht, Hajek = hj, AIPW = aipw, n_kept = sum(keep))
}
# Run the four scenarios with:
# ps_correct = D ~ X1 + X2; ps_wrong = D ~ X3;
# or_correct = Y ~ X1 + X2 + I(X1^2); or_wrong = Y ~ X3.
# Then: a per-scenario balance check (weighted SMD) — confirm balance FAILS when
# the PS is misspecified; and a PS-overlap histogram on the correct specification.# Reuses `df` from Exercise 2.
def ipw_aipw(df, ps_form, or_form):
ps_fit = # ___ YOUR CODE
e_hat = # ___ YOUR CODE
keep = # ___ YOUR CODE: trim at (0.02, 0.98)
d, e_hat = df[keep].reset_index(drop=True), e_hat[keep].reset_index(drop=True)
D, Y = d["D"].values, d["Y"].values
ht = # ___ YOUR CODE
hj = # ___ YOUR CODE
or1 = # ___ YOUR CODE: smf.ols(or_form, data=d[d.D == 1]).fit()
or0 = # ___ YOUR CODE
mu1, mu0 = or1.predict(d), or0.predict(d)
aipw = # ___ YOUR CODE
return {"HT": ht, "Hajek": hj, "AIPW": aipw, "n_kept": int(keep.sum())}
# Build the 4-scenario table; then the per-scenario weighted-SMD balance check.
TipShow solution
# Helper: compute Horvitz-Thompson IPW, Hájek IPW, and AIPW given PS and OR forms.
# Reuses `df` (the controlled DGP) from Exercise 2.
ipw_aipw <- function(df, ps_form, or_form) {
ps_fit <- glm(ps_form, data = df, family = binomial())
e_hat <- predict(ps_fit, type = "response")
keep <- e_hat > 0.02 & e_hat < 0.98 # trim to avoid divide-by-near-zero
d <- df[keep, ]; e_hat <- e_hat[keep]
D <- d$D; Y <- d$Y
ht <- mean(D * Y / e_hat) - mean((1 - D) * Y / (1 - e_hat))
hj <- sum(D * Y / e_hat) / sum(D / e_hat) -
sum((1 - D) * Y / (1 - e_hat)) / sum((1 - D) / (1 - e_hat))
or1 <- lm(or_form, data = d[d$D == 1, ])
or0 <- lm(or_form, data = d[d$D == 0, ])
mu1 <- predict(or1, newdata = d)
mu0 <- predict(or0, newdata = d)
aipw <- mean(mu1 - mu0 +
D * (Y - mu1) / e_hat -
(1 - D) * (Y - mu0) / (1 - e_hat))
tibble(HT = ht, Hajek = hj, AIPW = aipw, n_kept = sum(keep))
}ps_correct <- D ~ X1 + X2
ps_wrong <- D ~ X3 # Misspecified (omits true confounders)
or_correct <- Y ~ X1 + X2 + I(X1^2) # per-arm fit captures the D*X1 interaction
or_wrong <- Y ~ X3 # Misspecified
scenarios <- bind_rows(
cbind(scenario = "A: both right", ipw_aipw(df, ps_correct, or_correct)),
cbind(scenario = "B: OR wrong", ipw_aipw(df, ps_correct, or_wrong)),
cbind(scenario = "C: PS wrong", ipw_aipw(df, ps_wrong, or_correct)),
cbind(scenario = "D: both wrong", ipw_aipw(df, ps_wrong, or_wrong))
)
knitr::kable(scenarios, digits = 3,
caption = "IPW/AIPW point estimates by scenario. True ATE = 2.0.")| scenario | HT | Hajek | AIPW | n_kept |
|---|---|---|---|---|
| A: both right | 2.034 | 2.080 | 2.026 | 4997 |
| B: OR wrong | 2.034 | 2.080 | 2.091 | 4997 |
| C: PS wrong | 3.767 | 3.767 | 2.023 | 5000 |
| D: both wrong | 3.767 | 3.767 | 3.767 | 5000 |
# Q2: 95% bootstrap CIs for the AIPW estimate in each scenario, then a forest plot.
# Resample rows of `df` with replacement, recompute AIPW per scenario, B times.
set.seed(20260601L)
B <- 500
scen <- list(
"A: both right" = list(ps_correct, or_correct),
"B: OR wrong" = list(ps_correct, or_wrong),
"C: PS wrong" = list(ps_wrong, or_correct),
"D: both wrong" = list(ps_wrong, or_wrong)
)
aipw_ci <- function(ps_form, or_form) {
reps <- replicate(B, ipw_aipw(df[sample.int(nrow(df), replace = TRUE), ],
ps_form, or_form)$AIPW)
tibble(est = ipw_aipw(df, ps_form, or_form)$AIPW,
lo = quantile(reps, 0.025, names = FALSE),
hi = quantile(reps, 0.975, names = FALSE))
}
forest <- bind_rows(lapply(names(scen), function(nm) {
s <- scen[[nm]]
mutate(aipw_ci(s[[1]], s[[2]]), scenario = nm)
})) |>
mutate(scenario = factor(scenario, levels = rev(names(scen))),
covers = lo <= 2.0 & hi >= 2.0)
ggplot(forest, aes(x = est, y = scenario, colour = covers)) +
geom_vline(xintercept = 2.0, linetype = "dashed", colour = "grey35") +
geom_errorbarh(aes(xmin = lo, xmax = hi), height = 0.22, linewidth = 1) +
geom_point(size = 3) +
scale_colour_manual(values = c(`TRUE` = "#00C1D4", `FALSE` = "#C81E1E"),
labels = c(`TRUE` = "95% CI covers 2.0", `FALSE` = "CI misses"),
name = NULL) +
labs(x = "AIPW estimate with 95% bootstrap CI (true ATE = 2.0, dashed)",
y = NULL,
title = "Double robustness: AIPW brackets the truth when either model is right") +
theme_minimal()Warning: `geom_errorbarh()` was deprecated in ggplot2 4.0.0.
ℹ Please use the `orientation` argument of `geom_errorbar()` instead.`height` was translated to `width`.
# Per-scenario balance check: weighted SMD of X1, X2 under each scenario's PS
# weights. A failed weighted-balance check is the signature of a misspecified PS.
smd_w <- function(x, D, w) {
mt <- weighted.mean(x[D == 1], w[D == 1]); mc <- weighted.mean(x[D == 0], w[D == 0])
s <- sqrt((var(x[D == 1]) + var(x[D == 0])) / 2)
(mt - mc) / s
}
bal_scenario <- function(ps_form) {
e <- predict(glm(ps_form, data = df, family = binomial()), type = "response")
e <- pmin(pmax(e, 0.02), 0.98)
w <- ifelse(df$D == 1, 1 / e, 1 / (1 - e))
c(`|SMD X1|` = abs(smd_w(df$X1, df$D, w)), `|SMD X2|` = abs(smd_w(df$X2, df$D, w)))
}
tibble(
scenario = c("A/B: PS correct", "C/D: PS wrong"),
`PS model` = c("D ~ X1 + X2", "D ~ X3"),
bind_rows(bal_scenario(ps_correct), bal_scenario(ps_wrong))
) |>
knitr::kable(digits = 3,
caption = "Weighted covariate balance by PS specification (|SMD| < 0.10 = balanced). A wrong PS leaves the confounders grossly imbalanced — the diagnostic that tells you IPW is untrustworthy and the OR side must carry the AIPW estimate.")| scenario | PS model | |SMD X1| | |SMD X2| |
|---|---|---|---|
| A/B: PS correct | D ~ X1 + X2 | 0.018 | 0.023 |
| C/D: PS wrong | D ~ X3 | 0.699 | 0.584 |
# Propensity score overlap (with the correct PS specification)
ps_fit <- glm(D ~ X1 + X2, data = df, family = binomial())
df$ps <- predict(ps_fit, type = "response")
ggplot(df, aes(x = ps, fill = factor(D))) +
geom_histogram(alpha = 0.6, bins = 40, position = "identity") +
scale_fill_manual(values = c("0" = "#888", "1" = "#FF7E15"),
labels = c("Untreated", "Treated"),
name = NULL) +
labs(x = "Estimated propensity score",
y = "Count",
title = "Propensity overlap by treatment group") +
theme_minimal()
def ipw_aipw(df, ps_form, or_form):
ps_fit = smf.glm(ps_form, data=df,
family=sm.families.Binomial()).fit()
e_hat = ps_fit.predict(df)
keep = (e_hat > 0.02) & (e_hat < 0.98)
d, e_hat = df[keep].reset_index(drop=True), e_hat[keep].reset_index(drop=True)
D, Y = d["D"].values, d["Y"].values
ht = np.mean(D * Y / e_hat) - np.mean((1 - D) * Y / (1 - e_hat))
hj = (np.sum(D * Y / e_hat) / np.sum(D / e_hat)
- np.sum((1 - D) * Y / (1 - e_hat)) / np.sum((1 - D) / (1 - e_hat)))
or1 = smf.ols(or_form, data=d[d["D"] == 1]).fit()
or0 = smf.ols(or_form, data=d[d["D"] == 0]).fit()
mu1 = or1.predict(d); mu0 = or0.predict(d)
aipw = np.mean(mu1 - mu0 + D * (Y - mu1) / e_hat
- (1 - D) * (Y - mu0) / (1 - e_hat))
return {"HT": ht, "Hajek": hj, "AIPW": aipw, "n_kept": int(keep.sum())}
# `df` persists from the Ex 2 Python block (shared session).
scenarios = pd.DataFrame([
{"scenario": "A: both right",
**ipw_aipw(df, "D ~ X1 + X2", "Y ~ X1 + X2 + I(X1**2)")},
{"scenario": "B: OR wrong",
**ipw_aipw(df, "D ~ X1 + X2", "Y ~ X3")},
{"scenario": "C: PS wrong",
**ipw_aipw(df, "D ~ X3", "Y ~ X1 + X2 + I(X1**2)")},
{"scenario": "D: both wrong",
**ipw_aipw(df, "D ~ X3", "Y ~ X3")},
])
print(scenarios.round(3).to_string(index=False)) scenario HT Hajek AIPW n_kept
A: both right 2.053 2.038 2.077 5000
B: OR wrong 2.053 2.038 2.032 5000
C: PS wrong 3.728 3.728 2.078 5000
D: both wrong 3.728 3.728 3.728 5000# Q2: 95% bootstrap CIs for the AIPW estimate in each scenario, then a forest plot.
import matplotlib.pyplot as plt
rng_boot = np.random.default_rng(20260601)
B = 500
scen = {
"A: both right": ("D ~ X1 + X2", "Y ~ X1 + X2 + I(X1**2)"),
"B: OR wrong": ("D ~ X1 + X2", "Y ~ X3"),
"C: PS wrong": ("D ~ X3", "Y ~ X1 + X2 + I(X1**2)"),
"D: both wrong": ("D ~ X3", "Y ~ X3"),
}
def aipw_ci(ps_form, or_form):
reps = np.array([
ipw_aipw(df.iloc[rng_boot.integers(0, len(df), len(df))].reset_index(drop=True),
ps_form, or_form)["AIPW"]
for _ in range(B)
])
est = ipw_aipw(df, ps_form, or_form)["AIPW"]
return est, np.quantile(reps, 0.025), np.quantile(reps, 0.975)
rows = [(s, *aipw_ci(p, o)) for s, (p, o) in scen.items()]
y = np.arange(len(rows))[::-1] # scenario A at the top
fig, ax = plt.subplots(figsize=(7.5, 3.6))
ax.axvline(2.0, ls="--", color="#595959")
for yi, (s, est, lo, hi) in zip(y, rows):
col = "#00C1D4" if lo <= 2.0 <= hi else "#C81E1E"
ax.plot([lo, hi], [yi, yi], color=col, lw=2.5)
ax.plot(est, yi, "o", color=col, ms=8)
ax.set_yticks(y); ax.set_yticklabels([r[0] for r in rows])
ax.set_xlabel("AIPW estimate with 95% bootstrap CI (true ATE = 2.0, dashed)")
ax.set_title("Double robustness: AIPW brackets the truth when either model is right")
plt.tight_layout(); plt.show()
# Per-scenario weighted-SMD balance check (denominator = unweighted pooled SD).
Dv, X1v, X2v = df["D"].values, df["X1"].values, df["X2"].values
def smd_w(x, D, w):
mt = np.average(x[D == 1], weights=w[D == 1])
mc = np.average(x[D == 0], weights=w[D == 0])
s = np.sqrt((np.var(x[D == 1], ddof=1) + np.var(x[D == 0], ddof=1)) / 2)
return (mt - mc) / s
def bal_scenario(ps_form):
e = smf.glm(ps_form, data=df, family=sm.families.Binomial()).fit().predict(df).values
e = np.clip(e, 0.02, 0.98)
w = np.where(Dv == 1, 1 / e, 1 / (1 - e))
return abs(smd_w(X1v, Dv, w)), abs(smd_w(X2v, Dv, w))
rows = []
for label, form in [("A/B: PS correct", "D ~ X1 + X2"), ("C/D: PS wrong", "D ~ X3")]:
s1, s2 = bal_scenario(form)
rows.append({"scenario": label, "PS model": form,
"|SMD X1|": round(s1, 3), "|SMD X2|": round(s2, 3)})
print(pd.DataFrame(rows).to_string(index=False)) scenario PS model |SMD X1| |SMD X2|
A/B: PS correct D ~ X1 + X2 0.023 0.005
C/D: PS wrong D ~ X3 0.718 0.4719.4.4 Exercise 4 — Sensitivity analysis on the AIPW estimate (~25 min)
We now play the role of an analyst who built the AIPW estimate in Exercise 3 but did not have access to one of the true confounders (\(X_2\)). We compute two sensitivity quantities — the E-value and the sensemakr robustness value — and judge whether the analysis is robust. (Rosenbaum bounds apply to matched-pair estimators, not to this AIPW/OLS estimate.)
9.4.4.1 Scaffold
# (a) Restricted-adjustment AIPW on the controlled DGP (df, from Ex 2), withholding X2.
ps_restricted <- # ___ YOUR CODE: D ~ X1
or_restricted <- # ___ YOUR CODE: Y ~ X1
restricted <- # ___ YOUR CODE: ipw_aipw(df, ps_restricted, or_restricted)
# Print the restricted-AIPW point estimate.
# (b) Linear OLS for sensemakr (which needs an lm fit).
fit_for_sens <- # ___ YOUR CODE: lm(Y ~ D + X1, data = df)
# (c) sensemakr: robustness value + benchmark against X1 at kd = 1, 2, 3.
# Produce summary() and ovb_contour_plot().
library(sensemakr)
s <- # ___ YOUR CODE: sensemakr(fit_for_sens, treatment = "D", benchmark_covariates = "X1", kd = c(1, 2, 3))
# summary(s); ovb_contour_plot(s, sensitivity.of = "estimate")
# (d) E-value (VanderWeele-Ding 2017) via 0.91 * d_smd approximation.
library(EValue)
beta_hat <- coef(fit_for_sens)["D"]
sd_resid <- sigma(fit_for_sens)
ci <- confint(fit_for_sens)["D", ]
d_smd <- # ___ YOUR CODE: beta_hat / sd_resid
rr_hat <- # ___ YOUR CODE: exp(0.91 * d_smd)
# Compute lower / upper RR bounds analogously; call evalues.RR() and print.import numpy as np
import statsmodels.formula.api as smf
# Restricted-adjustment OLS (no X2)
fit_for_sens = # ___ YOUR CODE: smf.ols("Y ~ D + X1", data=df).fit()
# Restricted-adjustment AIPW (reuse the ipw_aipw helper from Ex 3 and df from Ex 2)
restricted = # ___ YOUR CODE: ipw_aipw(df, "D ~ X1", "Y ~ X1")
# E-value via the 0.91 * d approximation (no package needed)
beta_hat = # ___ YOUR CODE: fit_for_sens.params["D"]
sd_resid = # ___ YOUR CODE: np.sqrt(fit_for_sens.scale)
rr_hat = # ___ YOUR CODE: np.exp(0.91 * beta_hat / sd_resid)
# e_value(rr) = rr + sqrt(rr * (rr - 1)) for rr >= 1
# PySensemakr port for the robustness value + contour (pip install PySensemakr)
try:
import sensemakr as smkr
s = # ___ YOUR CODE: smkr.Sensemakr(model=fit_for_sens, treatment="D",
# benchmark_covariates=["X1"], kd=[1, 2, 3])
# print(s.summary()); s.plot(plot_type="contour", sensitivity_of="estimate")
except ImportError:
print("[PySensemakr not available; read the R sensemakr output instead]")
TipShow solution
Code
# Analyst's restricted adjustment set: X1 only (X2 is "unobserved")
ps_restricted <- D ~ X1
or_restricted <- Y ~ X1
restricted <- ipw_aipw(df, ps_restricted, or_restricted)
cat("Restricted-AIPW point estimate (true ATE = 2.0):", round(restricted$AIPW, 3), "\n")Restricted-AIPW point estimate (true ATE = 2.0): 2.643 Code
# Run a linear AIPW-flavoured regression for sensemakr (which works on lm fits)
fit_for_sens <- lm(Y ~ D + X1, data = df)
summary(fit_for_sens)$coefficients["D", ] Estimate Std. Error t value Pr(>|t|)
2.70668216 0.04764959 56.80389580 0.00000000 library(sensemakr)See details in:Carlos Cinelli and Chad Hazlett (2020). Making Sense of Sensitivity: Extending Omitted Variable Bias. Journal of the Royal Statistical Society, Series B (Statistical Methodology).s <- sensemakr(model = fit_for_sens, treatment = "D",
benchmark_covariates = "X1",
kd = c(1, 2, 3))Warning in ovb_partial_r2_bound.numeric(r2dxj.x = r2dxj.x[i], r2yxj.dx =
r2yxj.dx[i], : Implied bound on r2yz.dx greater than 1, try lower kd and/or ky.
Setting r2yz.dx to 1.summary(s)Sensitivity Analysis to Unobserved Confounding
Model Formula: Y ~ D + X1
Null hypothesis: q = 1 and reduce = TRUE
-- This means we are considering biases that reduce the absolute value of the current estimate.
-- The null hypothesis deemed problematic is H0:tau = 0
Unadjusted Estimates of 'D':
Coef. estimate: 2.7067
Standard Error: 0.0476
t-value (H0:tau = 0): 56.8039
Sensitivity Statistics:
Partial R2 of treatment with outcome: 0.3924
Robustness Value, q = 1: 0.5431
Robustness Value, q = 1, alpha = 0.05: 0.5312
Verbal interpretation of sensitivity statistics:
-- Partial R2 of the treatment with the outcome: an extreme confounder (orthogonal to the covariates) that explains 100% of the residual variance of the outcome, would need to explain at least 39.24% of the residual variance of the treatment to fully account for the observed estimated effect.
-- Robustness Value, q = 1: unobserved confounders (orthogonal to the covariates) that explain more than 54.31% of the residual variance of both the treatment and the outcome are strong enough to bring the point estimate to 0 (a bias of 100% of the original estimate). Conversely, unobserved confounders that do not explain more than 54.31% of the residual variance of both the treatment and the outcome are not strong enough to bring the point estimate to 0.
-- Robustness Value, q = 1, alpha = 0.05: unobserved confounders (orthogonal to the covariates) that explain more than 53.12% of the residual variance of both the treatment and the outcome are strong enough to bring the estimate to a range where it is no longer 'statistically different' from 0 (a bias of 100% of the original estimate), at the significance level of alpha = 0.05. Conversely, unobserved confounders that do not explain more than 53.12% of the residual variance of both the treatment and the outcome are not strong enough to bring the estimate to a range where it is no longer 'statistically different' from 0, at the significance level of alpha = 0.05.
Bounds on omitted variable bias:
--The table below shows the maximum strength of unobserved confounders with association with the treatment and the outcome bounded by a multiple of the observed explanatory power of the chosen benchmark covariate(s).
Bound Label R2dz.x R2yz.dx Treatment Adjusted Estimate Adjusted Se Adjusted T
1x X1 0.1175 1 D 1.4779 0 Inf
2x X1 0.2349 1 D 0.8402 0 Inf
3x X1 0.3524 1 D 0.2221 0 Inf
Adjusted Lower CI Adjusted Upper CI
1.4779 1.4779
0.8402 0.8402
0.2221 0.2221# X1 is so predictive of Y that a confounder "as strong as X1" implies a partial
# R^2 with the outcome > 1 (sensemakr caps it at 1), parking all three benchmark
# diamonds on the top edge — where sensemakr's own labels are clipped off the plot.
# Suppress them and label BELOW each diamond instead (carrying the adjusted D-effect).
ovb_contour_plot(s, sensitivity.of = "estimate", label.text = FALSE)
bnd <- s$bounds
text(bnd$r2dz.x, bnd$r2yz.dx,
labels = sprintf("%d× X1\n(%.2f)", seq_len(nrow(bnd)), bnd$adjusted_estimate),
pos = 1, offset = 0.5, cex = 0.8, font = 2)
library(EValue)
# Convert linear OLS estimate to risk-ratio scale via 0.91*d approximation
# (VanderWeele & Ding 2017): d = standardized mean difference; for a linear
# model, d ≈ beta / sigma_resid.
beta_hat <- coef(fit_for_sens)["D"]
sd_resid <- sigma(fit_for_sens)
ci_lo <- confint(fit_for_sens)["D", 1]
ci_hi <- confint(fit_for_sens)["D", 2]
d_smd <- beta_hat / sd_resid
d_lo <- ci_lo / sd_resid
d_hi <- ci_hi / sd_resid
rr_hat <- exp(0.91 * d_smd)
rr_lo <- exp(0.91 * d_lo)
rr_hi <- exp(0.91 * d_hi)
ev <- evalues.RR(est = rr_hat, lo = rr_lo, hi = rr_hi)
cat("Approximate risk-ratio (D->Y | X1):", round(rr_hat, 3), "\n")Approximate risk-ratio (D->Y | X1): 4.825 cat("E-value for point estimate: ", round(ev["E-values", "point"], 3), "\n")E-value for point estimate: 9.121 cat("E-value for the lower CI limit: ", round(ev["E-values", "lower"], 3), "\n")E-value for the lower CI limit: 8.609 import numpy as np
import statsmodels.formula.api as smf
# Restricted-adjustment OLS (no X2)
fit_for_sens = smf.ols("Y ~ D + X1", data=df).fit()
print(fit_for_sens.summary().tables[1])==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 0.7812 0.029 26.710 0.000 0.724 0.839
D 2.6664 0.047 56.470 0.000 2.574 2.759
X1 1.5632 0.023 67.558 0.000 1.518 1.609
==============================================================================# Restricted-adjustment AIPW — reuses the ipw_aipw helper from Ex 3 and df from Ex 2
# (both persist in Quarto's shared Python session). Mirrors the R tab's ~2.64.
restricted = ipw_aipw(df, "D ~ X1", "Y ~ X1")
print(f"Restricted-AIPW point estimate (true ATE = 2.0): {restricted['AIPW']:.3f}")Restricted-AIPW point estimate (true ATE = 2.0): 2.583# E-value (VanderWeele-Ding 2017) on the risk-ratio scale — language-agnostic,
# no package needed: convert the OLS estimate via the 0.91 * d approximation.
beta_hat = fit_for_sens.params["D"]
sd_resid = np.sqrt(fit_for_sens.scale)
ci_lo = fit_for_sens.conf_int().loc["D"].iloc[0]
rr_hat = np.exp(0.91 * beta_hat / sd_resid)
rr_lo = np.exp(0.91 * ci_lo / sd_resid)
def e_value(rr):
rr = rr if rr >= 1 else 1 / rr
return rr + np.sqrt(rr * (rr - 1))
print(f"Approximate risk-ratio (D->Y | X1): {rr_hat:.3f}")Approximate risk-ratio (D->Y | X1): 4.786print(f"E-value (point): {e_value(rr_hat):.3f} E-value (lower CI): {e_value(rr_lo):.3f}")E-value (point): 9.043 E-value (lower CI): 8.535# Robustness value + contour plot via the sensemakr port (pip install PySensemakr).
# If the package is unavailable, read the RV and contour plot from the R tab.
try:
import sensemakr as smkr
s = smkr.Sensemakr(
model=fit_for_sens, treatment="D",
benchmark_covariates=["X1"], kd=[1, 2, 3]
)
print(s.summary())
s.plot(plot_type="contour", sensitivity_of="estimate")
except Exception as e:
print(f"[PySensemakr not installed: {e}]")
print("[The robustness value and contour plot are in the R tab; the AIPW "
"estimate and E-value above already reproduce R's headline numbers.]")[PySensemakr not installed: No module named 'sensemakr']
[The robustness value and contour plot are in the R tab; the AIPW estimate and E-value above already reproduce R's headline numbers.]9.5 References & Resources
9.5.1 Academic literature
The full curated literature menu — 24 academic refs + 9 practitioner resources, every entry DOI-verified on 2026-05-20 — lives at literature/discovered/06_observed_confounding.md. Highlights below.
Foundational
- Rosenbaum & Rubin (1983), Biometrika — the central role of the propensity score (Rosenbaum and Rubin 1983)
- Dehejia & Wahba (2002), REStat — the canonical Lalonde reanalysis (Dehejia and Wahba 2002)
- Imbens (2004), REStat — review of nonparametric ATE under exogeneity
Matching
- Abadie & Imbens (2006), Econometrica — matching SEs (Abadie and Imbens 2006)
- Stuart (2010), Statistical Science — assignable single-paper overview (Stuart 2010)
- Iacus, King, Porro (2012), Political Analysis — CEM (Iacus, King, and Porro 2012)
- King & Nielsen (2019), Political Analysis — the PSM critique (King and Nielsen 2019)
Weighting
- Hirano, Imbens, Ridder (2003), Econometrica — IPW efficiency (Hirano, Imbens, and Ridder 2003)
- Hainmueller (2012), Political Analysis — entropy balancing (Hainmueller 2012)
- Imai & Ratkovic (2014), JRSS-B — CBPS (Imai and Ratkovic 2014)
- Zubizarreta (2015), JASA — stable balancing weights (Zubizarreta 2015)
AIPW and double robustness
- Robins, Rotnitzky, Zhao (1994), JASA — AIPW origin (Robins, Rotnitzky, and Zhao 1994)
Sensitivity
- Rosenbaum (2002), Springer — Observational Studies (2e), Rosenbaum bounds (Rosenbaum 2002)
- VanderWeele & Ding (2017), Annals of Internal Medicine — E-values (VanderWeele and Ding 2017)
- Cinelli & Hazlett (2020), JRSS-B — sensemakr (Cinelli and Hazlett 2020)
Reviews and overlap
- Athey & Imbens (2017), JEP — modern applied causal inference (Athey and Imbens 2017)
- Crump, Hotz, Imbens, Mitnik (2009), Biometrika — overlap trimming (Crump et al. 2009)
- Li, Morgan, Zaslavsky (2018), JASA — overlap weights (Li, Morgan, and Zaslavsky 2018)
- Cinelli, Forney, Pearl (2024), Sociological Methods & Research — good and bad controls (Cinelli, Forney, and Pearl 2024)
9.5.2 Practitioner resources
- Hernán & Robins, Causal Inference: What If — free PDF, chapters 12–15 cover observed confounding and IPW.
- Facure, Causal Inference for the Brave and True — practical PSM/IPW chapters with Python code.
MatchItvignette: https://kosukeimai.github.io/MatchIt/articles/MatchIt.htmlWeightItvignette: https://ngreifer.github.io/WeightIt/cobaltvignette: https://ngreifer.github.io/cobalt/sensemakrR vignette: https://carloscinelli.com/sensemakr/articles/sensemakr.html- PySensemakr quickstart: https://pysensemakr.readthedocs.io/en/latest/notebooks/quickstart.html
- DoWhy end-to-end pipeline: https://www.pywhy.org/dowhy/
- EconML observational estimators: https://econml.azurewebsites.net/
- DoubleML observational + DML: https://docs.doubleml.org/stable/index.html
9.5.3 Bilingual package summary
| Topic | R packages | Python packages |
|---|---|---|
| Matching | MatchIt, Matching, optmatch, cem |
psmpy, causalinference, dowhy |
| Balance diagnostics | cobalt, MatchIt |
dowhy refute, causallib |
| IPW / weighting | WeightIt, CBPS, ebal, survey, sbw |
dowhy, causallib, econml |
| AIPW / DR | drgee, AIPW, WeightIt (with OLS) |
econml.dr.LinearDRLearner, doubleml |
| Sensitivity | sensemakr, EValue, rbounds, sensitivitymv |
PySensemakr |
Abadie, Alberto. 2005. “Semiparametric Difference-in-Differences
Estimators.” Review of Economic Studies 72 (1): 1–19. https://doi.org/10.1111/0034-6527.00321.
———. 2021. “Using Synthetic Controls: Feasibility, Data
Requirements, and Methodological Aspects.” Journal of
Economic Literature 59 (2): 391–425. https://doi.org/10.1257/jel.20191450.
Abadie, Alberto, Anish Agarwal, Raaz Dwivedi, and Abhin Shah. 2024.
“Doubly Robust Inference in Causal Latent Factor Models.”
arXiv Preprint arXiv:2402.11652. https://doi.org/10.48550/arXiv.2402.11652.
Abadie, Alberto, Alexis Diamond, and Jens Hainmueller. 2010.
“Synthetic Control Methods for Comparative Case Studies:
Estimating the Effect of California’s Tobacco Control Program.”
Journal of the American Statistical Association 105 (490):
493–505. https://doi.org/10.1198/jasa.2009.ap08746.
———. 2015. “Comparative Politics and the Synthetic Control
Method.” American Journal of Political Science 59 (2):
495–510. https://doi.org/10.1111/ajps.12116.
Abadie, Alberto, and Guido W. Imbens. 2006. “Large Sample
Properties of Matching Estimators for Average Treatment Effects.”
Econometrica 74 (1): 235–67. https://doi.org/10.1111/j.1468-0262.2006.00655.x.
Abadie, Alberto, and Jann Spiess. 2022. “Robust Post-Matching
Inference.” Journal of the American Statistical
Association 117 (538): 983–95. https://doi.org/10.1080/01621459.2020.1840383.
Acemoglu, Daron, and Pascual Restrepo. 2020. “Robots and Jobs:
Evidence from US Labor Markets.” Journal of
Political Economy 128 (6): 2188–2244. https://doi.org/10.1086/705716.
Acharki, Naoufal, Ramiro Lugo, Antoine Bertoncello, and Josselin
Garnier. 2023. “Comparison of Meta-Learners for Estimating
Multi-Valued Treatment Heterogeneous Effects.” In Proceedings
of the 40th International Conference on Machine Learning
(ICML), 202:91–132. Proceedings of Machine Learning
Research.
Adão, Rodrigo, Michal Kolesár, and Eduardo Morales. 2019.
“Shift-Share Designs: Theory and Inference.” The
Quarterly Journal of Economics 134 (4): 1949–2010. https://doi.org/10.1093/qje/qjz025.
Ahrens, Achim, Christian B. Hansen, Mark E. Schaffer, and Thomas
Wiemann. 2025. “Model Averaging and Double Machine
Learning.” Journal of Applied Econometrics 40 (3):
249–69. https://doi.org/10.1002/jae.3103.
Anderson, T. W., and Herman Rubin. 1949. “Estimation of the
Parameters of a Single Equation in a Complete System of Stochastic
Equations.” The Annals of Mathematical Statistics 20
(1): 46–63. https://doi.org/10.1214/aoms/1177730090.
Andrews, Isaiah, James H. Stock, and Liyang Sun. 2019. “Weak
Instruments in Instrumental Variables Regression: Theory and
Practice.” Annual Review of Economics 11: 727–53. https://doi.org/10.1146/annurev-economics-080218-025643.
Angrist, Joshua D., Guido W. Imbens, and Donald B. Rubin. 1996.
“Identification of Causal Effects Using Instrumental
Variables.” Journal of the American Statistical
Association 91 (434): 444–55. https://doi.org/10.1080/01621459.1996.10476902.
Angrist, Joshua D., and Jörn-Steffen Pischke. 2015. Mastering
’Metrics: The Path from Cause to Effect. Princeton, NJ: Princeton
University Press.
Arkhangelsky, Dmitry, Susan Athey, David A. Hirshberg, Guido W. Imbens,
and Stefan Wager. 2021. “Synthetic
Difference-in-Differences.” American Economic Review 111
(12): 4088–118. https://doi.org/10.1257/aer.20190159.
Arkhangelsky, Dmitry, and Guido Imbens. 2024. “Causal Models for
Longitudinal and Panel Data: A Survey.” The
Econometrics Journal 27 (3): C1–61. https://doi.org/10.1093/ectj/utae014.
Arnold, David, Will Dobbie, and Crystal S. Yang. 2018. “Racial
Bias in Bail Decisions.” Quarterly Journal of Economics
133 (4): 1885–1932. https://doi.org/10.1093/qje/qjy012.
Athey, Susan, Mohsen Bayati, Nikolay Doudchenko, Guido Imbens, and
Khashayar Khosravi. 2021. “Matrix Completion Methods for Causal
Panel Data Models.” Journal of the American Statistical
Association 116 (536): 1716–30. https://doi.org/10.1080/01621459.2021.1891924.
Athey, Susan, and Guido Imbens. 2016. “Recursive Partitioning for
Heterogeneous Causal Effects.” Proceedings of the National
Academy of Sciences 113 (27): 7353–60. https://doi.org/10.1073/pnas.1510489113.
Athey, Susan, and Guido W. Imbens. 2017b. “The State of Applied
Econometrics: Causality and Policy Evaluation.” Journal of
Economic Perspectives 31 (2): 3–32. https://doi.org/10.1257/jep.31.2.3.
———. 2017a. “The State of Applied Econometrics:
Causality and Policy Evaluation.” Journal of
Economic Perspectives 31 (2): 3–32. https://doi.org/10.1257/jep.31.2.3.
Athey, Susan, Julie Tibshirani, and Stefan Wager. 2019.
“Generalized Random Forests.” The Annals of
Statistics 47 (2): 1148–78. https://doi.org/10.1214/18-AOS1709.
Athey, Susan, and Stefan Wager. 2021. “Policy Learning with
Observational Data.” Econometrica 89 (1): 133–61. https://doi.org/10.3982/ECTA15732.
Bach, Philipp, Victor Chernozhukov, Carlos Cinelli, Lin Jia, Sven
Klaassen, Nils Skotara, and Martin Spindler. 2025b. “Sensitivity
Analysis for Causal ML: A Use Case at Booking.com.” https://arxiv.org/abs/2510.09109.
———. 2025a. “Sensitivity Analysis for Causal ML:
A Use Case at Booking.com.” arXiv Preprint
arXiv:2510.09109. https://doi.org/10.48550/arXiv.2510.09109.
Bach, Philipp, Victor Chernozhukov, Malte S. Kurz, and Martin Spindler.
2022. “DoubleML – an Object-Oriented Implementation
of Double Machine Learning in Python.” Journal
of Machine Learning Research 23 (53): 1–6. https://jmlr.org/papers/v23/21-0862.html.
Bach, Philipp, Malte S. Kurz, Victor Chernozhukov, Martin Spindler, and
Sven Klaassen. 2024. “DoubleML — An
Object-Oriented Implementation of Double Machine Learning in
R.” Journal of Statistical Software 108
(3): 1–56. https://doi.org/10.18637/jss.v108.i03.
Baker, Andrew C., Brantly Callaway, Scott Cunningham, Andrew
Goodman-Bacon, and Pedro H. C. Sant’Anna. 2025.
“Difference-in-Differences: A Practitioner’s
Guide.” arXiv Preprint arXiv:2503.13323.
Ballinari, Daniele, and Nora Bearth. 2024. “Improving the Finite
Sample Estimation of Average Treatment Effects Using
Double/Debiased Machine Learning with Propensity Score
Calibration.” https://arxiv.org/abs/2409.04874.
Bang, Heejung, and James M. Robins. 2005. “Doubly Robust
Estimation in Missing Data and Causal Inference Models.”
Biometrics 61 (4): 962–73. https://doi.org/10.1111/j.1541-0420.2005.00377.x.
Bareinboim, Elias, Juan D. Correa, Duligur Ibeling, and Thomas Icard.
2022. “On Pearl’s Hierarchy and the
Foundations of Causal Inference.” Probabilistic and Causal
Inference: The Works of Judea Pearl, 507–56. https://doi.org/10.1145/3501714.3501743.
Bartik, Timothy J. 1991. Who Benefits from State and Local Economic
Development Policies? Kalamazoo, MI: W.E. Upjohn Institute for
Employment Research. https://doi.org/10.17848/9780585223940.
Ben-Michael, Eli, Avi Feller, and Jesse Rothstein. 2021. “The
Augmented Synthetic Control Method.” Journal of the American
Statistical Association 116 (536): 1789–1803. https://doi.org/10.1080/01621459.2021.1929245.
Borusyak, Kirill, Peter Hull, and Xavier Jaravel. 2022.
“Quasi-Experimental Shift-Share Research Designs.” The
Review of Economic Studies 89 (1): 181–213. https://doi.org/10.1093/restud/rdab030.
Borusyak, Kirill, Xavier Jaravel, and Jann Spiess. 2024.
“Revisiting Event-Study Designs: Robust and Efficient
Estimation.” Review of Economic Studies 91 (6): 3253–85.
https://doi.org/10.1093/restud/rdae007.
Bound, John, David A. Jaeger, and Regina M. Baker. 1995. “Problems
with Instrumental Variables Estimation When the Correlation Between the
Instruments and the Endogenous Explanatory Variable Is Weak.”
Journal of the American Statistical Association 90 (430):
443–50. https://doi.org/10.1080/01621459.1995.10476536.
Caetano, Carolina, Brantly Callaway, Stroud Payne, and Hugo Sant’Anna
Rodrigues. 2024. “Difference in Differences with Time-Varying
Covariates.” https://doi.org/10.48550/arXiv.2202.02903.
Callaway, Brantly, and Pedro H. C. Sant’Anna. 2021.
“Difference-in-Differences with Multiple Time Periods.”
Journal of Econometrics 225 (2): 200–230. https://doi.org/10.1016/j.jeconom.2020.12.001.
Calonico, Sebastian, Matias D. Cattaneo, and Rocío Titiunik. 2014.
“Robust Nonparametric Confidence Intervals for
Regression-Discontinuity Designs.” Econometrica 82 (6):
2295–2326. https://doi.org/10.3982/ECTA11757.
Cattaneo, Matias D., Yingjie Feng, Filippo Palomba, and Rocío Titiunik.
2025. “Scpi: Uncertainty Quantification for Synthetic Control
Methods.” Journal of Statistical Software 113 (1): 1–38.
https://doi.org/10.18637/jss.v113.i01.
Cattaneo, Matias D., Nicolás Idrobo, and Rocío Titiunik. 2020. A
Practical Introduction to Regression Discontinuity Designs:
Foundations. Cambridge Elements in Quantitative and Computational
Methods for the Social Sciences. Cambridge: Cambridge University Press.
https://doi.org/10.1017/9781108684606.
Cattaneo, Matias D., Michael Jansson, and Xinwei Ma. 2020. “Simple
Local Polynomial Density Estimators.” Journal of the American
Statistical Association 115 (531): 1449–55. https://doi.org/10.1080/01621459.2019.1635480.
Chaisemartin, Clément de, and Xavier D’Haultfœuille. 2020.
“Two-Way Fixed Effects Estimators with Heterogeneous Treatment
Effects.” American Economic Review 110 (9): 2964–96. https://doi.org/10.1257/aer.20181169.
———. 2023. “Two-Way Fixed Effects and Differences-in-Differences
with Heterogeneous Treatment Effects: A Survey.”
The Econometrics Journal 26 (3): C1–30. https://doi.org/10.1093/ectj/utac017.
Chang, Neng-Chieh. 2020. “Double/Debiased Machine Learning for
Difference-in-Differences Models.” The Econometrics
Journal 23 (2): 177–91. https://doi.org/10.1093/ectj/utaa001.
Chernozhukov, Victor, Denis Chetverikov, Mert Demirer, Esther Duflo,
Christian Hansen, Whitney Newey, and James Robins. 2018.
“Double/Debiased Machine Learning for Treatment and Structural
Parameters.” The Econometrics Journal 21 (1): C1–68. https://doi.org/10.1111/ectj.12097.
Chernozhukov, Victor, Carlos Cinelli, Whitney K. Newey, Amit Sharma, and
Vasilis Syrgkanis. 2026. “Long Story Short: Omitted
Variable Bias in Causal Machine Learning.” Review of
Economics and Statistics, 1–45. https://doi.org/10.1162/rest.a.1705.
Chernozhukov, Victor, Mert Demirer, Esther Duflo, and Iván
Fernández-Val. 2018. “Fisher-Schultz Lecture: Generic
Machine Learning Inference on Heterogenous Treatment Effects in
Randomized Experiments, with an Application to Immunization in
India.” https://arxiv.org/abs/1712.04802.
Chernozhukov, Victor, Whitney K. Newey, Víctor Quintas-Martínez, and
Vasilis Syrgkanis. 2021. “Automatic Debiased Machine Learning via
Riesz Regression.” https://arxiv.org/abs/2104.14737.
Chernozhukov, Victor, Whitney K. Newey, and Rahul Singh. 2022.
“Automatic Debiased Machine Learning of Causal and Structural
Effects.” Econometrica 90 (3): 967–1027. https://doi.org/10.3982/ECTA18515.
Chernozhukov, Victor, Kaspar Wüthrich, and Yinchu Zhu. 2021. “An
Exact and Robust Conformal Inference Method for Counterfactual and
Synthetic Controls.” Journal of the American Statistical
Association 116 (536): 1849–64. https://doi.org/10.1080/01621459.2021.1920957.
Chickering, David Maxwell. 2002. “Optimal Structure Identification
with Greedy Search.” Journal of Machine Learning
Research 3: 507–54.
Chiu, Albert, Xingchen Lan, Ziyi Liu, and Yiqing Xu. 2026. “Causal
Panel Analysis Under Parallel Trends: Lessons from a Large Reanalysis
Study.” American Political Science Review 120 (1):
245–66. https://doi.org/10.1017/S0003055425000243.
Choi, Jungjun, Hyukjun Kwon, and Yuan Liao. 2024. “Inference for
Low-Rank Completion Without Sample Splitting with Application to
Treatment Effect Estimation.” Journal of Econometrics
240 (1): 105682. https://doi.org/10.1016/j.jeconom.2024.105682.
Choi, Jungjun, and Ming Yuan. 2024. “Matrix Completion When
Missing Is Not at Random and Its Applications in Causal Panel Data
Models.” Journal of the American Statistical
Association. https://doi.org/10.1080/01621459.2024.2380105.
Chyn, Eric, Brigham R. Frandsen, and Emily C. Leslie. 2025.
“Examiner and Judge Designs in Economics: A Practitioner’s
Guide.” Journal of Economic Literature 63 (2): 401–39.
https://doi.org/10.1257/jel.20241719.
Cinelli, Carlos, Andrew Forney, and Judea Pearl. 2024. “A Crash
Course in Good and Bad Controls.” Sociological Methods &
Research 53 (3): 1071–1104. https://doi.org/10.1177/00491241221099552.
Cinelli, Carlos, and Chad Hazlett. 2020. “Making Sense of
Sensitivity: Extending Omitted Variable Bias.” Journal of the
Royal Statistical Society: Series B 82 (1): 39–67. https://doi.org/10.1111/rssb.12348.
Crump, Richard K., V. Joseph Hotz, Guido W. Imbens, and Oscar A. Mitnik.
2009. “Dealing with Limited Overlap in Estimation of Average
Treatment Effects.” Biometrika 96 (1): 187–99. https://doi.org/10.1093/biomet/asn055.
Cunningham, Scott. 2021. Causal Inference: The Mixtape. Yale
University Press. https://mixtape.scunning.com/.
Curth, Alicia, and Mihaela van der Schaar. 2023. “In Search of
Insights, Not Magic Bullets: Towards Demystification of the
Model Selection Dilemma in Heterogeneous Treatment Effect
Estimation.” In Proceedings of the 40th International
Conference on Machine Learning (ICML), 202:6623–42.
Proceedings of Machine Learning Research.
Dandl, Susanne, Christian Haslinger, Torsten Hothorn, Heidi Seibold,
Erik Sverdrup, Stefan Wager, and Achim Zeileis. 2024. “What Makes
Forest-Based Heterogeneous Treatment Effect Estimators Work?”
The Annals of Applied Statistics 18 (1): 506–28. https://doi.org/10.1214/23-AOAS1799.
Dehejia, Rajeev H., and Sadek Wahba. 2002. “Propensity
Score-Matching Methods for Nonexperimental Causal Studies.”
Review of Economics and Statistics 84 (1): 151–61. https://doi.org/10.1162/003465302317331982.
Deng, Alex, Ya Xu, Ron Kohavi, and Toby Walker. 2013. “Improving
the Sensitivity of Online Controlled Experiments by Utilizing
Pre-Experiment Data.” In Proceedings of the Sixth ACM
International Conference on Web Search and Data Mining (WSDM),
123–32. https://doi.org/10.1145/2433396.2433413.
Dobbie, Will, Andres Liberman, Daniel Paravisini, and Vikram Pathania.
2021. “Measuring Bias in Consumer Lending.” The Review
of Economic Studies 88 (6): 2799–2832. https://doi.org/10.1093/restud/rdaa057.
Elbers, Benjamin. 2023. “Encouragement Designs and Instrumental
Variables for a/b Testing.” Spotify Engineering. 2023. https://engineering.atspotify.com/2023/08/encouragement-designs-and-instrumental-variables-for-a-b-testing.
Fabijan, Aleksander, Jayant Gupchup, Somit Gupta, Jeff Omhover, Wen Qin,
Lukas Vermeer, and Pavel Dmitriev. 2019. “Diagnosing Sample Ratio
Mismatch in Online Controlled Experiments: A Taxonomy and
Rules of Thumb for Practitioners.” In Proceedings of the 25th
ACM SIGKDD International Conference on Knowledge Discovery and Data
Mining (KDD), 2156–64. https://doi.org/10.1145/3292500.3330722.
Ferman, Bruno, and Cristine Pinto. 2021. “Synthetic Controls with
Imperfect Pretreatment Fit.” Quantitative Economics 12
(4): 1197–1221. https://doi.org/10.3982/QE1596.
Fernandes, Ana, Martin Huber, and Giannina Vaccaro. 2020. “Gender
Differences in Wage Expectations.” https://arxiv.org/abs/2003.11496.
Fernández-Loría, Carlos, and Foster Provost. 2022. “Causal
Decision Making and Causal Effect Estimation Are Not the Same…and Why It
Matters.” INFORMS Journal on Data Science 1 (1): 4–16.
https://doi.org/10.1287/ijds.2021.0006.
Feuerriegel, Stefan, Dennis Frauen, Valentyn Melnychuk, Jonas
Schweisthal, Konstantin Hess, Alicia Curth, Stefan Bauer, Niki
Kilbertus, Isaac S. Kohane, and Mihaela van der Schaar. 2024.
“Causal Machine Learning for Predicting Treatment
Outcomes.” Nature Medicine 30: 958–68. https://doi.org/10.1038/s41591-024-02902-1.
Frandsen, Brigham R., Lars J. Lefgren, and Emily C. Leslie. 2023.
“Judging Judge Fixed Effects.” American Economic
Review 113 (1): 253–77. https://doi.org/10.1257/aer.20201860.
Frisch, Ragnar, and Frederick V. Waugh. 1933. “Partial Time
Regressions as Compared with Individual Trends.”
Econometrica 1 (4): 387–401. https://doi.org/10.2307/1907330.
Frölich, Markus. 2007. “Nonparametric IV Estimation
of Local Average Treatment Effects with Covariates.” Journal
of Econometrics 139 (1): 35–75. https://doi.org/10.1016/j.jeconom.2006.06.004.
Gardner, John. 2022. “Two-Stage Differences in
Differences.” https://doi.org/10.48550/arXiv.2207.05943.
Gelman, Andrew, and Guido Imbens. 2019. “Why High-Order
Polynomials Should Not Be Used in Regression Discontinuity
Designs.” Journal of Business & Economic Statistics
37 (3): 447–56. https://doi.org/10.1080/07350015.2017.1366909.
Goldin, Jacob, Julian Nyarko, and Justin Young. 2024. “Forecasting
Algorithms for Causal Inference with Panel Data.” arXiv
Preprint arXiv:2208.03489. https://doi.org/10.48550/arXiv.2208.03489.
Goldsmith-Pinkham, Paul, Isaac Sorkin, and Henry Swift. 2020.
“Bartik Instruments: What, When, Why, and How.”
American Economic Review 110 (8): 2586–2624. https://doi.org/10.1257/aer.20181047.
Goodman-Bacon, Andrew. 2021. “Difference-in-Differences with
Variation in Treatment Timing.” Journal of Econometrics
225 (2): 254–77. https://doi.org/10.1016/j.jeconom.2021.03.014.
Gordon, Brett R., Florian Zettelmeyer, Neha Bhargava, and Dan Chapsky.
2019. “A Comparison of Approaches to Advertising Measurement:
Evidence from Big Field Experiments at Facebook.” Marketing
Science 38 (2): 193–225. https://doi.org/10.1287/mksc.2018.1135.
Guo, Yongyi, Dominic Coey, Mikael Konutgan, Wenting Li, Chris Schoener,
and Matt Goldman. 2021. “Machine Learning for Variance Reduction
in Online Experiments.” In Advances in Neural Information
Processing Systems 34 (NeurIPS). https://proceedings.neurips.cc/paper/2021/hash/488b084119a1c7a4950f00706ec7ea16-Abstract.html.
Hahn, Jinyong, Petra Todd, and Wilbert Van der Klaauw. 2001.
“Identification and Estimation of Treatment Effects with a
Regression-Discontinuity Design.” Econometrica 69 (1):
201–9. https://doi.org/10.1111/1468-0262.00183.
Hainmueller, Jens. 2012. “Entropy Balancing for Causal Effects: A
Multivariate Reweighting Method to Produce Balanced Samples in
Observational Studies.” Political Analysis 20 (1):
25–46. https://doi.org/10.1093/pan/mpr025.
Heckman, James J., Hidehiko Ichimura, and Petra E. Todd. 1997.
“Matching as an Econometric Evaluation Estimator: Evidence from
Evaluating a Job Training Programme.” Review of Economic
Studies 64 (4): 605–54. https://doi.org/10.2307/2971733.
Hirano, Keisuke, Guido W. Imbens, and Geert Ridder. 2003.
“Efficient Estimation of Average Treatment Effects Using the
Estimated Propensity Score.” Econometrica 71 (4):
1161–89. https://doi.org/10.1111/1468-0262.00442.
Howard, Steven R., Aaditya Ramdas, Jon McAuliffe, and Jasjeet Sekhon.
2021. “Time-Uniform, Nonparametric, Nonasymptotic Confidence
Sequences.” The Annals of Statistics 49 (2): 1055–80. https://doi.org/10.1214/20-AOS1991.
Huber, Martin. 2023. Causal Analysis: Impact Evaluation and Causal
Machine Learning with Applications in R. Cambridge,
MA: MIT Press.
Hünermund, Paul, and Elias Bareinboim. 2023. “Causal Inference and
Data Fusion in Econometrics.” The Econometrics Journal
26 (1): 1–34. https://doi.org/10.1093/ectj/utad008.
Iacus, Stefano M., Gary King, and Giuseppe Porro. 2012. “Causal
Inference Without Balance Checking: Coarsened Exact Matching.”
Political Analysis 20 (1): 1–24. https://doi.org/10.1093/pan/mpr013.
Imai, Kosuke, and In Song Kim. 2019. “When Should We Use Unit
Fixed Effects Regression Models for Causal Inference with Longitudinal
Data?” American Journal of Political Science 63 (2):
467–90. https://doi.org/10.1111/ajps.12417.
Imai, Kosuke, and Marc Ratkovic. 2014. “Covariate Balancing
Propensity Score.” Journal of the Royal Statistical Society:
Series B 76 (1): 243–63. https://doi.org/10.1111/rssb.12027.
Imbens, Guido W., and Joshua D. Angrist. 1994. “Identification and
Estimation of Local Average Treatment Effects.”
Econometrica 62 (2): 467–75. https://doi.org/10.2307/2951620.
Imbens, Guido W., and Thomas Lemieux. 2008. “Regression
Discontinuity Designs: A Guide to Practice.” Journal of
Econometrics 142 (2): 615–35. https://doi.org/10.1016/j.jeconom.2007.05.001.
Jacob, Daniel. 2021. “CATE Meets ML:
The Conditional Average Treatment Effect and Machine
Learning.” Digital Finance 3 (2): 99–148. https://doi.org/10.1007/s42521-021-00033-7.
Kennedy, Edward H. 2023. “Towards Optimal Doubly Robust Estimation
of Heterogeneous Causal Effects.” Electronic Journal of
Statistics 17 (2): 3008–49. https://doi.org/10.1214/23-EJS2157.
King, Gary, and Richard Nielsen. 2019. “Why Propensity Scores
Should Not Be Used for Matching.” Political Analysis 27
(4): 435–54. https://doi.org/10.1017/pan.2019.11.
Kleibergen, Frank, and Richard Paap. 2006. “Generalized Reduced
Rank Tests Using the Singular Value Decomposition.” Journal
of Econometrics 133 (1): 97–126. https://doi.org/10.1016/j.jeconom.2005.02.011.
Knaus, Michael C. 2022. “Double Machine Learning-Based Programme
Evaluation Under Unconfoundedness.” The Econometrics
Journal 25 (3): 602–27. https://doi.org/10.1093/ectj/utac015.
Kohavi, Ron, Diane Tang, and Ya Xu. 2020. Trustworthy Online
Controlled Experiments: A Practical Guide to a/b Testing. Cambridge
University Press. https://doi.org/10.1017/9781108653985.
Künzel, Sören R., Jasjeet S. Sekhon, Peter J. Bickel, and Bin Yu. 2019.
“Metalearners for Estimating Heterogeneous Treatment Effects Using
Machine Learning.” Proceedings of the National Academy of
Sciences 116 (10): 4156–65. https://doi.org/10.1073/pnas.1804597116.
Lechner, Michael. 2018. “Modified Causal Forests for Estimating
Heterogeneous Causal Effects.” https://arxiv.org/abs/1812.09487.
Lechner, Michael, and Jana Mareckova. 2024. “Comprehensive Causal
Machine Learning.” https://arxiv.org/abs/2405.10198.
Lee, David S., and Thomas Lemieux. 2010. “Regression Discontinuity
Designs in Economics.” Journal of Economic Literature 48
(2): 281–355. https://doi.org/10.1257/jel.48.2.281.
Lee, David S., Justin McCrary, Marcelo J. Moreira, and Jack Porter.
2022. “Valid t-Ratio Inference for IV.” American
Economic Review 112 (10): 3260–90. https://doi.org/10.1257/aer.20211063.
Li, Fan, Kari Lock Morgan, and Alan M. Zaslavsky. 2018. “Balancing
Covariates via Propensity Score Weighting.” Journal of the
American Statistical Association 113 (521): 390–400. https://doi.org/10.1080/01621459.2016.1260466.
Liu, Licheng, Ye Wang, and Yiqing Xu. 2024. “A Practical Guide to
Counterfactual Estimators for Causal Inference with Time-Series
Cross-Sectional Data.” American Journal of Political
Science 68 (1): 160–76. https://doi.org/10.1111/ajps.12723.
Lovell, Michael C. 1963. “Seasonal Adjustment of Economic Time
Series and Multiple Regression Analysis.” Journal of the
American Statistical Association 58 (304): 993–1010. https://doi.org/10.1080/01621459.1963.10480682.
Manski, Charles F. 2003. Partial Identification of Probability
Distributions. Springer Series in Statistics. New York:
Springer-Verlag. https://doi.org/10.1007/b97478.
McCrary, Justin. 2008. “Manipulation of the Running Variable in
the Regression Discontinuity Design: A Density Test.” Journal
of Econometrics 142 (2): 698–714. https://doi.org/10.1016/j.jeconom.2007.05.005.
Montiel Olea, José Luis, and Carolin Pflueger. 2013. “A Robust
Test for Weak Instruments.” Journal of Business &
Economic Statistics 31 (3): 358–69. https://doi.org/10.1080/00401706.2013.806694.
Nie, Xinkun, and Stefan Wager. 2021. “Quasi-Oracle Estimation of
Heterogeneous Treatment Effects.” Biometrika 108 (2):
299–319. https://doi.org/10.1093/biomet/asaa076.
Pearl, Judea. 2009. Causality: Models, Reasoning, and
Inference. 2nd ed. Cambridge University Press. https://doi.org/10.1017/CBO9780511803161.
Pocock, Stuart J. 1977. “Group Sequential Methods in the Design
and Analysis of Clinical Trials.” Biometrika 64 (2):
191–99. https://doi.org/10.1093/biomet/64.2.191.
Qian, Min, and Susan A. Murphy. 2011. “Performance Guarantees for
Individualized Treatment Rules.” The Annals of
Statistics 39 (2): 1180–210. https://doi.org/10.1214/10-AOS864.
Rambachan, Ashesh, and Jonathan Roth. 2023. “A More Credible
Approach to Parallel Trends.” Review of Economic Studies
90 (5): 2555–91. https://doi.org/10.1093/restud/rdad018.
Rehill, Patrick. 2025. “How Do Applied Researchers Use the Causal
Forest? A Methodological Review.” International
Statistical Review 93 (2): 288–316. https://doi.org/10.1111/insr.12610.
Robins, James M., Andrea Rotnitzky, and Lue Ping Zhao. 1994.
“Estimation of Regression Coefficients When Some Regressors Are
Not Always Observed.” Journal of the American Statistical
Association 89 (427): 846–66. https://doi.org/10.1080/01621459.1994.10476818.
Robinson, Peter M. 1988. “Root-N-Consistent
Semiparametric Regression.” Econometrica 56 (4): 931–54.
https://doi.org/10.2307/1912705.
Rosenbaum, Paul R. 2002. Observational Studies. 2nd ed.
Springer Series in Statistics. Springer. https://doi.org/10.1007/978-1-4757-3692-2.
Rosenbaum, Paul R., and Donald B. Rubin. 1983. “The Central Role
of the Propensity Score in Observational Studies for Causal
Effects.” Biometrika 70 (1): 41–55. https://doi.org/10.1093/biomet/70.1.41.
Roth, Jonathan. 2022. “Pretest with Caution:
Event-Study Estimates After Testing for Parallel
Trends.” American Economic Review: Insights 4 (3):
305–22. https://doi.org/10.1257/aeri.20210236.
———. 2024. “Interpreting Event-Studies from Recent
Difference-in-Differences Methods.” arXiv Preprint
arXiv:2401.12309.
Roth, Jonathan, and Pedro H. C. Sant’Anna. 2023. “When Is Parallel
Trends Sensitive to Functional Form?” Econometrica 91
(2): 737–47. https://doi.org/10.3982/ECTA19402.
Sant’Anna, Pedro H. C., and Jun Zhao. 2020. “Doubly Robust
Difference-in-Differences Estimators.” Journal of
Econometrics 219 (1): 101–22. https://doi.org/10.1016/j.jeconom.2020.06.003.
Semenova, Vira, and Victor Chernozhukov. 2021. “Debiased Machine
Learning of Conditional Average Treatment Effects and Other Causal
Functions.” The Econometrics Journal 24 (2): 264–89. https://doi.org/10.1093/ectj/utaa027.
Sharma, Amit, Jake M. Hofman, and Duncan J. Watts. 2015.
“Estimating the Causal Impact of Recommendation Systems from
Observational Data.” In Proceedings of the Sixteenth ACM
Conference on Economics and Computation (EC ’15), 453–70. https://doi.org/10.1145/2764468.2764488.
Shchetkina, Anya, and Ron Berman. 2024. “When Is Heterogeneity
Actionable for Personalization?” https://arxiv.org/abs/2411.16552.
Shimizu, Shohei, Patrik O. Hoyer, Aapo Hyvärinen, and Antti Kerminen.
2006. “A Linear Non-Gaussian Acyclic Model for Causal
Discovery.” Journal of Machine Learning Research 7:
2003–30.
Shpitser, Ilya, and Judea Pearl. 2006. “Identification of Joint
Interventional Distributions in Recursive Semi-Markovian
Causal Models.” Proceedings of the 21st National Conference
on Artificial Intelligence (AAAI), 1219–26.
Spirtes, Peter, Clark Glymour, and Richard Scheines. 2000.
Causation, Prediction, and Search. 2nd ed. Cambridge, MA: MIT
Press.
Staiger, Douglas, and James H. Stock. 1997. “Instrumental
Variables Regression with Weak Instruments.”
Econometrica 65 (3): 557–86. https://doi.org/10.2307/2171753.
Stock, James H., and Motohiro Yogo. 2005. “Testing for Weak
Instruments in Linear IV Regression.” In Identification and
Inference for Econometric Models: Essays in Honor of Thomas
Rothenberg, edited by Donald W. K. Andrews and James H. Stock,
80–108. Cambridge: Cambridge University Press. https://doi.org/10.1017/CBO9780511614491.006.
Stuart, Elizabeth A. 2010. “Matching Methods for Causal Inference:
A Review and a Look Forward.” Statistical Science 25
(1): 1–21. https://doi.org/10.1214/09-STS313.
Sun, Liyang, and Sarah Abraham. 2021. “Estimating Dynamic
Treatment Effects in Event Studies with Heterogeneous Treatment
Effects.” Journal of Econometrics 225 (2): 175–99. https://doi.org/10.1016/j.jeconom.2020.09.006.
Sverdrup, Erik, Han Wu, Susan Athey, and Stefan Wager. 2024. “Qini
Curves for Multi-Armed Treatment Rules.” Journal of
Computational and Graphical Statistics. https://doi.org/10.1080/10618600.2024.2418820.
VanderWeele, Tyler J., and Peng Ding. 2017. “Sensitivity Analysis
in Observational Research: Introducing the E-Value.”
Annals of Internal Medicine 167 (4): 268–74. https://doi.org/10.7326/M16-2607.
Veitch, Victor, and Anisha Zaveri. 2020. “Sense and Sensitivity
Analysis: Simple Post-Hoc Analysis of Bias Due to Unobserved
Confounding.” In Advances in Neural Information Processing
Systems (NeurIPS), 33:10999–1009.
Wager, Stefan, and Susan Athey. 2018. “Estimation and Inference of
Heterogeneous Treatment Effects Using Random Forests.”
Journal of the American Statistical Association 113 (523):
1228–42. https://doi.org/10.1080/01621459.2017.1319839.
Waudby-Smith, Ian, and Aaditya Ramdas. 2024. “Estimating Means of
Bounded Random Variables by Betting.” Journal of the Royal
Statistical Society: Series B 86 (1): 1–27. https://doi.org/10.1093/jrsssb/qkad009.
Wooldridge, Jeffrey M. 2010. Econometric Analysis of Cross Section
and Panel Data. 2nd ed. Cambridge, MA: MIT Press.
Xu, Yiqing. 2017. “Generalized Synthetic Control Method: Causal
Inference with Interactive Fixed Effects Models.” Political
Analysis 25 (1): 57–76. https://doi.org/10.1017/pan.2016.2.
Yadlowsky, Steve, Scott Fleming, Nigam Shah, Emma Brunskill, and Stefan
Wager. 2025. “Evaluating Treatment Prioritization Rules via
Rank-Weighted Average Treatment Effects.” Journal of the
American Statistical Association 120 (549): 38–51. https://doi.org/10.1080/01621459.2024.2393466.
Zeldow, Bret, and Laura A. Hatfield. 2021. “Confounding and
Regression Adjustment in Difference-in-Differences Studies.”
Health Services Research 56 (5): 932–41. https://doi.org/10.1111/1475-6773.13666.
Zhou, Zhengyuan, Susan Athey, and Stefan Wager. 2023. “Offline
Multi-Action Policy Learning: Generalization and
Optimization.” Operations Research 71 (1): 148–83. https://doi.org/10.1287/opre.2022.2271.
Zubizarreta, José R. 2015. “Stable Weights That Balance Covariates
for Estimation with Incomplete Outcome Data.” Journal of the
American Statistical Association 110 (511): 910–22. https://doi.org/10.1080/01621459.2015.1023805.