Advanced Causal Inference

The Potential Outcomes Framework defines causal effects through counterfactuals. For each unit i and treatment T ∈ {0,1}, Y_i(1) is the potential outcome under treatment and Y_i(0) under control. The individual treatment effect is τ_i = Y_i(1) - Y_i(0), but the fundamental problem of causal inference is that only one potential outcome is observed. Average Treatment Effect ATE = E[Y(1) - Y(0)]. ATT = E[Y(1) - Y(0) | T=1]. Identification requires assumptions (ignorability, SUTVA) to link observed data to potential outcomes. Randomization makes assignment independent of potential outcomes, enabling unbiased ATE estimation.

# ATE estimation from RCT
ate = y[treatment==1].mean() - y[treatment==0].mean()

A DAG represents causal structure as directed edges (X → Y means X causes Y) with no directed cycles. Nodes are variables; edges are causal relationships. DAGs encode conditional independence relationships via the d-separation criterion: two nodes are d-separated given a conditioning set Z if all paths between them are blocked. They identify adjustment sets needed for causal identification, expose collider bias, and guide instrument and mediator selection. The do-calculus (Pearl) uses DAGs to express interventions: P(Y | do(X)) vs P(Y | X). DAGs make assumptions explicit, enabling rigorous evaluation of whether causal effects are identifiable from observational data.

The backdoor criterion provides a sufficient condition for identifying the causal effect of X on Y from observational data. A set of variables Z satisfies the backdoor criterion for (X, Y) if: (1) Z blocks all backdoor paths from X to Y (paths with an arrow into X), and (2) Z contains no descendant of X. Conditioning on such a Z removes confounding: P(Y|do(X)) = Σ_Z P(Y|X,Z)P(Z). Finding a valid adjustment set requires the DAG structure. Not all confounders need to be included — only a valid adjustment set. The front-door criterion handles cases where no valid backdoor adjustment set exists but a mediator is available.

The front-door criterion identifies the causal effect of X on Y when no valid backdoor adjustment set exists (unobserved confounders between X and Y) but a mediator M is available such that: (1) M blocks all directed paths from X to Y (M fully mediates the effect), (2) there are no unblocked backdoor paths from X to M, and (3) all backdoor paths from M to Y are blocked by X. The front-door formula: P(Y|do(X)) = Σ_M P(M|X) * Σ_X' P(Y|M,X')P(X'). Classic example: smoking (X) → tar deposits (M) → cancer (Y) with unobserved genetic confounders between smoking and cancer. The front-door estimator exploits the mediation chain.

IV estimation identifies causal effects when treatment is endogenous (correlated with unobserved confounders). A valid instrument Z must satisfy: (1) relevance — Z is correlated with treatment X (testable); (2) exclusion restriction — Z affects outcome Y only through X, not directly (untestable, requires theory); (3) independence — Z is independent of unmeasured confounders. The IV estimator (2SLS): regress X on Z (first stage), use fitted values X̂ as exogenous variation, regress Y on X̂. IV estimates the Local Average Treatment Effect (LATE) for compliers — those whose treatment status is changed by Z. Examples: distance to college as IV for education, lottery as IV for military service.

from linearmodels.iv import IV2SLS
result = IV2SLS(dependent=df['y'],
                exog=df[['const']],
                endog=df['x'],
                instruments=df['z']).fit()

RDD identifies causal effects by exploiting a sharp cutoff in a continuous assignment variable (running variable) that determines treatment eligibility. Units just above the threshold are treated; those just below are not. Near the threshold, units are similar in expectation (local randomization), so comparing outcomes just above vs just below estimates the Local Average Treatment Effect at the cutoff. Sharp RDD: all units follow the threshold rule. Fuzzy RDD: some units don't comply, requiring IV estimation using threshold crossing as an instrument. Validity assumptions: no manipulation of the running variable (density test), no other discontinuities at the threshold, and continuity of potential outcomes in the running variable.

import rdd
rdd_result = rdd.rdd(df['outcome'], df['running_var'],
                      cutpoint=0, bw=5)
# McCrary density test for manipulation
rdd.density_test(df['running_var'], cutpoint=0)

DiD estimates the causal effect of a treatment that is applied to some units (treated group) at a specific time, by comparing the change in outcome over time for treated vs control groups. DiD estimate = (treated_after - treated_before) - (control_after - control_before). Key assumption: parallel trends — in the absence of treatment, the treated and control groups would have followed the same trend over time. This is untestable for the post-treatment period but testable in pre-treatment periods. Common violations: differential seasonality, treatment selection correlated with pre-trends, and anticipation effects. Modern extensions handle staggered adoption, heterogeneous treatment timing, and non-parallel pre-trends.

import statsmodels.formula.api as smf
result = smf.ols('outcome ~ treated + post + treated:post + controls',
                  data=df).fit()

Traditional DiD with a single treatment period is well-understood, but staggered adoption (units adopt treatment at different times) creates complications: the canonical TWFE (two-way fixed effects) estimator can produce negative weights on some treatment effects, potentially yielding a negative aggregate estimate even when all individual ATTs are positive (the Goodman-Bacon decomposition reveals this). Callaway and Sant'Anna (2021) propose group-time ATT(g,t) estimands for each cohort-time pair (never-treated or not-yet-treated as comparison), then aggregate them. Estimators: 'not yet treated' comparison group, 'never treated' comparison group. The did R/Python package implements these with influence-function standard errors.

# Python: csdid package
from csdid import ATTgt
att = ATTgt(y='outcome', t='time', id='unit', d='treatment_time', data=df)
att.fit()

Synthetic Control constructs a weighted combination of control units to serve as a counterfactual for a single treated unit, matching on pre-treatment outcome trends and covariates. Weights are chosen to minimize the discrepancy between the treated unit's pre-treatment outcomes and the weighted average of controls. The causal effect is the difference between actual and synthetic control outcomes post-treatment. It is ideal for case studies with one or few treated units (e.g., a policy change in one state, a product launch in one country). Inference uses permutation tests (applying the method to each control unit as placebo). The Synth R package and Python synth_runner implement it.

Propensity score matching (PSM) estimates the ATE or ATT by matching treated and control units with similar probabilities of treatment (propensity score), estimated via logistic regression on observed covariates. Matching removes confounding from observed covariates. The conditional independence assumption requires that assignment is ignorable given observed covariates — no unmeasured confounders. Limitations: PSM does not handle unobserved confounders; common support violations occur when treated and control don't overlap in propensity score (extrapolation); model misspecification in propensity estimation biases results; poor matches inflate variance. Doubly robust estimation (AIPW) combines propensity and outcome models for robustness to misspecification of one model.

from sklearn.linear_model import LogisticRegression
ps_model = LogisticRegression().fit(X_confounders, T)
propensity_scores = ps_model.predict_proba(X_confounders)[:,1]

The Augmented Inverse Probability Weighting (AIPW) estimator combines IPW and outcome regression, achieving double robustness: it is consistent if either (but not necessarily both) the propensity score model or the outcome model is correctly specified. The AIPW estimator: τ_AIPW = E[Y(1) - Y(0)] = E[(T/e(X)) * (Y - μ_1(X)) + μ_1(X)] - E[((1-T)/(1-e(X))) * (Y - μ_0(X)) + μ_0(X)], where e(X) is the propensity and μ_t(X) is the outcome model. Combined with cross-fitting (sample splitting) and machine learning for nuisance estimation, it achieves semiparametric efficiency bounds. Implemented in EconML and DoubleML packages.

from econml.dr import DRLearner
model = DRLearner()
model.fit(Y, T, X=X, W=W)
ate = model.ate(X)

IPW creates a pseudo-population in which treatment assignment is independent of observed confounders by weighting each unit by the inverse of its propensity score: w_i = T_i/e(X_i) + (1-T_i)/(1-e(X_i)). Treated units with low propensity (unlikely to be treated but were) get high weight; control units with high propensity (likely to be treated but weren't) get high weight. The weighted sample mimics randomization. IPW estimates are unbiased if the propensity model is correctly specified and there is adequate common support (no extreme propensity scores near 0 or 1). Trimming or truncating extreme weights stabilizes variance at the cost of some bias. Stabilized weights divide by the marginal treatment probability.

weights = T / propensity + (1-T) / (1-propensity)
# Weighted regression
import statsmodels.api as sm
wls = sm.WLS(Y, sm.add_constant(T), weights=weights).fit()

Mediation analysis decomposes the total causal effect of X on Y into: the Natural Direct Effect (NDE) — effect of X on Y not through the mediator M — and the Natural Indirect Effect (NIE) — effect operating through M. Total Effect = NDE + NIE. Baron and Kenny's three-step approach tests mediation via sequential regression, but it cannot handle non-linear models or confounders of the mediator-outcome relationship. The potential outcomes framework defines natural effects precisely: NDE = E[Y(1, M(0))] - E[Y(0, M(0))]; NIE = E[Y(1, M(1))] - E[Y(1, M(0))]. Identification requires no confounders of X→Y, X→M, and M→Y — including no mediator-outcome confounders affected by X.

from causalinference import CausalModel
# Or using Python's causalpy for mediation
import pingouin as pg
mediation = pg.mediation_analysis(data=df, x='X', m='M', y='Y', n_boot=1000)

IV estimation identifies the Local Average Treatment Effect (LATE) — the average causal effect for compliers: units whose treatment status changes because of the instrument (they comply with the instrument's push). It does not identify the ATE for the whole population or the ATT for those who chose treatment regardless of the instrument. Three subpopulations in an IV setting: always-takers (treated regardless of Z), never-takers (never treated regardless of Z), and compliers (treatment changes with Z). Defiers (treatment is opposite to Z) must not exist for monotonicity assumption. LATE = ITT_Y / ITT_X, where ITT_Y is the reduced-form effect of Z on Y and ITT_X is the first-stage effect of Z on X.

# LATE = Wald estimator
late = (y[z==1].mean() - y[z==0].mean()) / (d[z==1].mean() - d[z==0].mean())

A placebo test applies the causal method to a context where the treatment effect should be zero, and checks whether the estimated effect is indeed near zero. If the placebo test produces a spurious significant effect, the method's assumptions are likely violated. Examples: testing the DiD estimator in a pre-treatment period (should find no effect if parallel trends holds), applying IV estimation with a falsified instrument that shouldn't affect the outcome, applying RDD above/below a different threshold where no program exists, or testing synthetic control on a unit unaffected by the treatment. Placebo tests are essential for validating the credibility of quasi-experimental designs.

Sensitivity analysis quantifies how strong an unmeasured confounder would need to be to explain away an observed treatment effect. Rosenbaum's sensitivity analysis for matching studies asks: if units with identical observed covariates could differ in their odds of treatment by a factor Gamma, how large would Gamma need to be to make the treatment effect statistically insignificant? A study is robust if it requires an implausibly large Gamma (e.g., Gamma = 5, meaning 5x odds difference from unmeasured confounders) to overturn the conclusion. Frank's KonFound analysis and Cinelli-Hazlett's sensitivity framework provide modern alternatives expressing sensitivity in terms of partial R-squared.

# R: rbounds package
# Python equivalent
from pyrosenbaum import sensitivity
result = sensitivity(matched_pairs, gamma=2.0)

CATE = E[Y(1) - Y(0) | X = x] is the average treatment effect conditional on covariates X, capturing treatment effect heterogeneity across subgroups. ATE is the average of CATE over the covariate distribution. CATE estimation requires flexible machine learning methods: T-learner (separate outcome models for treated/control), S-learner (single model with treatment as a feature), X-learner (cross-fitting for imbalanced populations), and the R-learner (orthogonal ML based on the Robinson decomposition). The causalml and EconML Python libraries implement these. CATE is used for personalization (whom to treat?), optimal policy design, and heterogeneous effects analysis.

from econml.dml import CausalForestDML
cf = CausalForestDML()
cf.fit(Y, T, X=X, W=W)
cate_hat = cf.effect(X_test)

The R-learner (Robinson 1988, Chernozhukov et al. 2018 Double ML) estimates CATE by exploiting the partially linear model Y = θ(X)*T + g(X) + ε. Step 1: estimate residuals from nuisance models (ε̃ = Y - E[Y|X], T̃ = T - E[T|X]) using cross-fitting to avoid regularization bias. Step 2: regress ε̃ on T̃ using a flexible model with T̃ as weights, targeting θ(X). Cross-fitting (k-fold) prevents overfitting of nuisance models from contaminating the CATE estimate. The R-learner achieves the semiparametric efficiency bound when nuisance models converge at o(N^{-1/4}) rates — satisfied by many ML estimators. Implemented in EconML's DML and CausalForestDML.

from econml.dml import LinearDML
model = LinearDML(model_y='auto', model_t='auto', cv=5)
model.fit(Y, T, X=X, W=W)
print(model.ate(X))

Causal discovery aims to infer the causal graph (DAG) from observational data, not just associations. Key algorithms: PC algorithm (tests conditional independence to identify the skeleton, then orients edges using v-structures and orientation rules); FCI (Fast Causal Inference) extends PC to handle latent confounders, producing a PAG; GES (Greedy Equivalence Search) uses a score-based approach. LiNGAM (Linear Non-Gaussian Acyclic Models) exploits non-Gaussianity of noise to fully orient the DAG. Assumptions required: causal Markov condition, faithfulness (no cancellation of paths), and causal sufficiency (no unmeasured common causes, except FCI). Causal discovery is computationally expensive and sensitive to violations of faithfulness.

import causallearn.search.ConstraintBased.PC as pc_alg
cg = pc_alg.pc(data, alpha=0.05, indep_test='fisherz')

The Markov boundary of a variable Y is the minimal set of variables that renders Y conditionally independent of all other variables — it consists of Y's parents, children, and spouses (co-parents of Y's children) in the causal DAG. Conditioning on the Markov boundary is the minimal sufficient set for predicting Y. For causal feature selection, the Markov boundary provides the optimal prediction set: adding more variables does not improve prediction of Y given the boundary. Algorithms like IAMB (Incremental Association Markov Boundary) and HITON-MB discover the Markov boundary from data. This has applications in feature selection, variable importance, and causal structure learning.

The exclusion restriction requires that the instrument Z affects the outcome Y only through the treatment X, with no direct effect. It is fundamentally untestable from observational data (unlike relevance, which is testable via first-stage F-statistic). Partial testing strategies: (1) Apply IV to outcomes that should not be affected by the treatment — a placebo outcome test; (2) Check that the instrument has no effect in a subpopulation where treatment cannot change (never-takers); (3) Test the instrument's validity in a comparison period before treatment was possible; (4) Economic/theoretical reasoning about why the instrument would have no direct effect. Overidentification tests (Sargan-Hansen J test) test whether multiple instruments are mutually consistent but cannot prove exclusion.

# First-stage F-statistic (should be > 10)
import statsmodels.formula.api as smf
first_stage = smf.ols('X ~ Z + controls', data=df).fit()
print(f'F-stat: {first_stage.fvalue:.1f}')

A weak instrument has low correlation with the treatment (low first-stage F-statistic), causing the 2SLS estimator to be severely biased toward OLS and producing wide, unreliable confidence intervals. The conventional rule-of-thumb threshold is F > 10 (Stock-Yogo). With weak instruments, t-statistics for IV estimates are unreliable because the IV distribution becomes non-normal even in large samples. Solutions: use Anderson-Rubin (AR) test for inference, which is robust to weak instruments; Stock-Wright S-statistic; or conditional likelihood ratio (CLR) tests. Including multiple weak instruments does not solve the problem and can worsen the bias-variance trade-off. Finding a strong instrument is the primary solution.

# F > 10 rule of thumb
first_stage = smf.ols('X ~ Z + controls', data=df).fit()
print(f'First-stage F: {first_stage.fvalue:.1f}')  # needs > 10

ATT = E[Y(1) - Y(0) | T=1] is the average causal effect among those who actually received the treatment. ATE = E[Y(1) - Y(0)] averages over the full population. ATT is the appropriate estimand when: (1) the treatment cannot be assigned to the full population (e.g., surgical procedure — relevant only for those who qualify), (2) you want to evaluate a program by how well it worked for participants, or (3) there is limited overlap between treated and control distributions (common support issues make ATE estimation unreliable). Propensity score matching naturally estimates ATT by matching controls to treated units. IV estimates LATE for compliers, which is neither ATE nor ATT.

A collider is a variable C that is caused by both X and Y (both X → C ← Y). Conditioning on a collider — by selecting on it, controlling for it in regression, or studying a population defined by it — opens a spurious association between X and Y even if they are causally independent. Example: among hospitalized patients (selecting on hospitalization C), a spurious negative association appears between disease A and disease B because both cause hospitalization. In DAG terms, conditioning on C unblocks the path X → C ← Y, creating a backdoor path. Collider bias is a major source of spurious correlations in observational studies, clinical studies, and machine learning feature selection.

The sharp null hypothesis (Fisher) is that the treatment has absolutely no effect on any individual unit: Y_i(1) = Y_i(0) for all i. Under the sharp null, all potential outcomes are known from observed data, enabling exact permutation-based inference. The weak null hypothesis (Neyman) is that the average treatment effect is zero: E[Y(1) - Y(0)] = 0, but individual effects may be heterogeneous and cancel out. t-tests and z-tests test the weak null. A test can reject the weak null but not the sharp null (positive and negative effects that cancel), or reject the sharp null via distributional test while ATE = 0. These distinctions matter for interpreting experiment results under treatment heterogeneity.

The Angrist-Pischke (2009) first-stage check evaluates instrument relevance in IV regressions with multiple endogenous variables. For a single endogenous variable, the relevant statistic is the first-stage F-statistic for the excluded instruments. The Stock-Yogo critical values for 2SLS bias and size distortion depend on the number of instruments and the desired maximal bias/size. F < 10 indicates weak instruments. With multiple endogenous variables, the AP F-statistic applies separately to each endogenous variable — it is the F-statistic from regressing the residuals of one endogenous variable on the excluded instruments, conditioning on the other endogenous variables. Robust to heteroskedasticity using Kleibergen-Paap F-statistic.

# Kleibergen-Paap F-statistic in Python
from linearmodels.iv import IV2SLS
result = IV2SLS(y, X_exog, X_endog, Z).fit(cov_type='robust')
print(result.first_stage.diagnostics)

The Heckman two-step estimator corrects for sample selection bias when the outcome is only observed for a non-random subset of units. Step 1: probit model for selection into the observed sample, computing the inverse Mills ratio (IMR) = φ(Xδ)/Φ(Xδ) for each unit. Step 2: OLS regression of the outcome on treatment and covariates plus the IMR as a control variable. The IMR corrects for the correlation between unobservables in the selection and outcome equations. Requires an exclusion restriction: a variable in the selection equation but not the outcome equation (similar to IV). The exclusion restriction is often difficult to defend, making Heckit results sensitive to specification.

from heckman import Heckman
model = Heckman(outcome='y', selection='selected', exog=['x1','x2'],
                exog_select=['x1','x2','z'], data=df)
result = model.fit()

Consistency (also called SUTVA part 2) requires that the observed outcome equals the potential outcome under the received treatment: if T_i = t, then Y_i = Y_i(t). Violations occur when there are multiple versions of treatment — e.g., 'drug' could be delivered at different doses, by different clinicians, or in different formulations. Each version may have different effects, so Y_i(treatment) is not well-defined without specifying the version. In observational studies, treatment heterogeneity (different types of exposure) is a common violation. Solutions: define treatment precisely at the version level, or use principal stratification to handle multiple versions. In online experiments, consistency is often satisfied by design through precise feature specification.

A confounder causes both the treatment and the outcome, creating a spurious association (e.g., age causes both exercise and health). Conditioning on confounders removes this backdoor bias. A mediator lies on the causal path between treatment and outcome (X → M → Y), transmitting part of X's effect. Conditioning on a mediator blocks the indirect effect, biasing total effect estimates — mediators should not be included in regression for estimating total effects. A collider is caused by both X and Y; conditioning on it opens a spurious path. These distinctions are only possible to make through a causal DAG — statistical criteria alone (e.g., correlation with treatment) do not distinguish them.

Under the Neyman potential outcomes framework for experiments, the only source of randomness is the randomization mechanism itself — potential outcomes are treated as fixed constants, not random variables. The sampling distribution of estimators (e.g., difference in means) is derived from the distribution over all possible randomization realizations. This 'design-based' inference contrasts with 'model-based' inference where randomness comes from assumed population sampling. Design-based variance for the difference in means is Var(Y(1))/n_1 + Var(Y(0))/n_0 - Var(Y(1)-Y(0))/N. Since the covariance term is unobservable, the standard Neyman variance estimator is conservative (overestimates variance when individual treatment effects are homogeneous).

A natural experiment exploits exogenous variation in treatment assignment that occurs naturally (not controlled by the researcher) — policy changes, geographic boundaries, weather shocks, administrative cutoffs, lottery outcomes. Unlike designed experiments with explicit randomization, natural experiments require arguing that the variation is 'as-if' random conditional on some controls. Common natural experiments: draft lottery (as-if random military assignment), weather shocks affecting agricultural output, proximity to schools affecting educational outcomes. The credibility of the as-if randomization argument determines the validity of causal conclusions. Designed experiments have higher internal validity by construction; natural experiments may have higher external validity if the population is more representative.

Local Projections (LP, Jordà 2005) estimate impulse response functions (dynamic causal effects over multiple horizons) by running a sequence of regressions: for each horizon h, regress Y_{t+h} on treatment T_t and controls. Unlike VAR impulse responses, LP does not impose parametric restrictions across horizons — each regression is estimated independently, making LP more robust to model misspecification. LP standard errors are computed using Newey-West to account for serial correlation. In policy evaluation, LP estimates how a treatment at time t affects outcomes at t+1, t+2, ..., t+H, providing the full dynamic response path. Used in macroeconomics and increasingly in tech platform experiments for studying temporal treatment effects.

import numpy as np
# LP: regress Y_{t+h} on T_t for h = 1,...,H
for h in range(1, H+1):
    y_h = df['Y'].shift(-h)
    reg = smf.ols('y_h ~ T + controls', data=df).fit(cov_type='HAC', cov_kwds={'maxlags':h})

Angrist and Krueger (1991) used quarter of birth as an instrument for years of education to estimate the return to schooling. The argument: compulsory schooling laws require students to stay in school until a certain age; students born early in the year enter school younger and can drop out with less education. Quarter of birth (Q1 vs Q4) affects education but should not directly affect earnings — satisfying the exclusion restriction. IV estimates showed returns of ~8% per year of schooling. This study is canonical for: demonstrating a creative natural experiment IV strategy, illustrating the LATE interpretation (effect for compulsory-schooling-constrained students), and sparking debates about instrument validity and weak instruments (Staiger-Stock 1997).

ANCOVA (Analysis of Covariance) includes pre-experiment covariates as additional regressors in the treatment effect regression: Y = alpha + tau*T + beta*X + epsilon. This reduces residual variance attributable to X, tightening confidence intervals and increasing power without collecting more data. Unlike CUPED which adjusts the outcome variable, ANCOVA directly controls for covariates in the regression. Both are equivalent under linear relationships. The inclusion of covariates does not bias ATE estimates in randomized experiments (Lin 2013) — Lin's estimator interacts treatment with demeaned covariates for robustness. Variance reduction equals the R-squared from regressing Y on X in the control group: Var_reduction = R^2 * Var(Y).

import statsmodels.formula.api as smf
result = smf.ols('Y ~ T + X1 + X2 + T:X1 + T:X2', data=df).fit()
print(result.params['T'])  # Lin estimator for ATE

A DiD event study estimates treatment effects at each time period relative to treatment adoption, plotting ATT(t) for t = ..., -3, -2, -1, 0, +1, +2, ... where t=0 is treatment. Pre-treatment coefficients (t < 0) should be near zero if parallel trends holds — non-zero pre-trends indicate the assumption is violated. Post-treatment coefficients show the dynamic evolution of the effect. The regression includes unit and time fixed effects plus interactions of treatment indicator with time period dummies. This visualization is more informative than a single DiD estimate and is standard in economics journals. Dynamic effects (growing or decaying ATT over time) reveal treatment effect heterogeneity over the post-treatment horizon.

# Event study regression
result = smf.ols('Y ~ C(unit) + C(time) + C(relative_time):treated',
                  data=df).fit(cov_type='cluster', cov_kwds={'groups':df['unit']})

The fundamental problem of causal inference (Holland 1986) is that for any individual unit, only one potential outcome can be observed — the outcome under the treatment actually received. The counterfactual potential outcome (what would have happened under the alternative treatment) is fundamentally unobservable. Therefore, individual causal effects τ_i = Y_i(1) - Y_i(0) are never directly observable. All causal methods work around this by making assumptions that allow identification of average effects (ATE, ATT) from observable data: randomization, ignorability given observed covariates, IV assumptions, or structural assumptions from theory. This fundamentally distinguishes causal inference from prediction — prediction requires only P(Y|X), causation requires P(Y|do(X)).

Selection on observables (conditional ignorability) assumes that after conditioning on observed covariates X, treatment assignment is as good as random: T ⊥ {Y(0), Y(1)} | X. This justifies propensity score methods, regression adjustment, and matching. Selection on unobservables means assignment depends on unobserved factors correlated with potential outcomes — the 'endogeneity' problem that propensity score methods cannot handle. Methods for unobserved selection: IV (exploits exogenous variation), DiD (differences out time-invariant unobserved factors), RDD (local randomization at cutoff), fixed effects (removes unit-specific unobservables). The choice of method depends on the assumed selection mechanism, which must be justified by domain knowledge and theory.

When CATE(x) = E[Y(1)-Y(0)|X=x] varies across individuals, a uniform treatment policy is suboptimal. Optimal policy assigns treatment to units with positive expected benefit: π*(x) = 1 if CATE(x) > cost_of_treatment. Causal forests (Wager-Athey) estimate CATE with honest trees (splitting sample into discovery and estimation halves to prevent overfitting). Policy learning (Kitagawa-Tetenov, Athey-Wager) learns the best treatment assignment rule from data. Evaluation uses doubly robust welfare estimators. Applications: targeted interventions (who benefits from a drug?), personalized pricing, resource allocation (where to deploy social services?). The gain from targeting vs uniform policy equals the integral of positive CATE mass.

from econml.grf import CausalForest
cf = CausalForest(n_estimators=500, min_impurity_decrease=0.001)
cf.fit(X, T, Y)
cate = cf.predict(X_test)
policy = (cate > treatment_cost).astype(int)

ATU = E[Y(1) - Y(0) | T=0] is the average causal effect that untreated units would have experienced if they had received the treatment. ATE = ATT * P(T=1) + ATU * P(T=0). ATU is the relevant estimand when evaluating the effect of expanding a program to non-participants: would the job-training program benefit those who didn't enroll? Under selection on observables, ATU is identifiable. The difference ATT - ATU (treatment effect heterogeneity between choosers and non-choosers) reflects Roy model selection: individuals self-select into treatment partly based on their expected benefit, making ATT > ATU when selection is on gains.

Internal validity is the degree to which the study correctly identifies the causal effect within the study population and context, free from confounding, selection bias, and measurement error. A well-designed RCT has high internal validity. External validity (generalizability) is the degree to which findings from the study population transfer to other populations, settings, and time periods. An RCT run on volunteers in a clinical setting may have high internal validity but low external validity if patients in the general population differ systematically. Quasi-experimental designs often have higher external validity (natural populations) but lower internal validity. The LATE from IV estimates only the effect on compliers, limiting external validity to that subpopulation.

Principal stratification (Frangakis-Rubin 2002) defines subgroups based on the joint potential values of an intermediate variable — creating strata that are functions of potential outcomes and hence pre-treatment. The four principal strata for a binary treatment and binary intermediate: always-takers, never-takers, compliers, and defiers. Causal effects within strata are well-defined causal contrasts (principal causal effects) because strata are defined by potential outcomes, not observed post-treatment variables. It is used for: handling non-compliance in RCTs (CACE = effect for compliers), handling truncation by death (survivor average causal effect), and surrogate endpoint validation in clinical trials.

The Goodman-Bacon (2021) decomposition shows that the two-way fixed effects (TWFE) DiD estimator with staggered treatment adoption is a weighted average of all possible 2x2 DiD comparisons: early-treated vs control, late-treated vs control, and — problematically — early-treated vs late-treated (where the already-treated group serves as the 'control'). The last type uses units that are already treated as controls, producing negative weights when treatment effects are heterogeneous over time. This can make the TWFE estimate negative even if all individual treatment effects are positive. The decomposition motivates using staggered-design-robust estimators (Callaway-Sant'Anna, Sun-Abraham, Borusyak et al.).

Conformal prediction is a distribution-free framework for constructing prediction sets with guaranteed coverage: for a given confidence level 1-alpha, the prediction set covers the true outcome with probability >= 1-alpha without distributional assumptions, for any sample size. The split conformal method: fit a model on the training set, compute nonconformity scores (e.g., |y - y_hat|) on a calibration set, and set the prediction interval width at the (1-alpha)(1+1/n_cal) quantile of calibration scores. The coverage guarantee is marginal (averaged over test points). It is model-agnostic and works for regression, classification, and causal prediction problems. Widely used in ML deployment for risk-calibrated predictions.

from mapie.regression import MapieRegressor
from sklearn.linear_model import Ridge
mapie = MapieRegressor(Ridge(), method='base', cv='split')
mapie.fit(X_train, y_train)
y_pred, y_pi = mapie.predict(X_test, alpha=0.1)  # 90% coverage

In an RCT, randomization ensures treatment independence from potential outcomes: E[Y(0)|T=1] = E[Y(0)|T=0], so the difference in observed means is an unbiased ATE estimator. No covariate adjustment is required for consistency, though it improves efficiency. In observational studies, self-selection creates correlation between T and potential outcomes through confounders. ATE identification requires conditional ignorability (T ⊥ Y(0),Y(1) | X) plus common support. Estimators (IPW, AIPW, matching) must correct for confounding. Even with correct adjustment, observational ATE estimates may be biased by unmeasured confounders that RCTs automatically balance. RCTs are the gold standard for internal validity; observational studies compensate with larger, more representative samples.

Difference-in-discontinuities (DiD-RD) combines the RDD and DiD frameworks to identify causal effects when a policy changes the discontinuity at a known threshold across time. The design compares the discontinuity at the cutoff in the post-policy period to the discontinuity in the pre-policy period. The identifying assumption is that any pre-existing discontinuity at the threshold (not due to the policy) is stable over time. This relaxes the RDD continuity assumption by allowing for pre-existing jumps at the cutoff, as long as these jumps are constant over time (analogous to the parallel trends assumption in DiD). Used when a program rule applied at a cutoff changes its assignment formula at a specific point in time.

Angrist, Imbens, and Rubin (1996) extended IV analysis to settings with heterogeneous treatment effects, showing that 2SLS estimates the LATE (Local Average Treatment Effect) for compliers — units whose treatment status is determined by the instrument. Under monotonicity (no defiers: the instrument moves all units in the same direction), the LATE is identified by the Wald estimator: LATE = ITT_Y / ITT_D. The framework formalizes four assumptions: SUTVA, random assignment of instrument, exclusion restriction, and monotonicity. The LATE is policy-relevant when the instrument mimics the variation of interest (e.g., a draft lottery identifies the effect of military service for those whose service was determined by the lottery).

# Wald estimator for LATE
ITT_Y = y[z==1].mean() - y[z==0].mean()
ITT_D = d[z==1].mean() - d[z==0].mean()
LATE = ITT_Y / ITT_D
print(f'LATE: {LATE:.3f}')

TMLE (van der Laan & Rubin 2006) is a doubly robust semiparametric estimator that combines machine learning nuisance estimation with a targeting step that optimizes the bias-variance trade-off for a specific target parameter (e.g., ATE). The algorithm: (1) estimate the outcome model E[Y|T,X] and propensity score P(T=1|X) using flexible ML; (2) compute a clever covariate H and 'target' the outcome model by fitting epsilon in a one-step logistic regression update; (3) plug in the updated predictions to compute the ATE. TMLE achieves the semiparametric efficiency bound and is multiply robust (consistent if either the outcome or propensity model is correctly specified). Implemented in Python's causalml and R's tmle3.

from causalml.inference.meta import TMLELearner
tmle = TMLELearner()
tmle.fit(y, X, treatment)
ate, lb, ub = tmle.estimate_ate()

Matching estimators achieve identification by finding control units that are 'close' to each treated unit in covariate space. The bias of a matching estimator (nearest-neighbor match) is O(N^{-r/d}) where r is the order of the match and d is the number of covariates. As dimensionality d increases, the required sample size for bias control grows exponentially — the curse of dimensionality. With 10 binary covariates, the matching space has 2^10 = 1024 cells, making exact matching increasingly sparse. Solutions: propensity score matching reduces the matching problem to a one-dimensional scalar; exact matching on a subset of key variables with regression adjustment for the rest; or doubly robust estimators that use flexible ML for nuisance functions.

Interrupted Time Series (ITS) is a quasi-experimental design that estimates a treatment effect by modeling the pre-intervention trend and extrapolating it as the counterfactual, then comparing to the observed post-intervention series. It requires many pre- and post-intervention time points and assumes the pre-trend would have continued unchanged without intervention. The regression model includes time, treatment indicator, and a time-since-treatment slope change: Y_t = beta0 + beta1*t + beta2*D_t + beta3*(t-T0)*D_t + e_t, where D_t is post-intervention indicator. The coefficients capture immediate level change (beta2) and slope change (beta3). Controls for seasonality, autocorrelation (ARIMA errors), and concurrent interventions are critical for validity.

import statsmodels.formula.api as smf
df['post'] = (df['time'] >= intervention_time).astype(int)
df['time_since'] = (df['time'] - intervention_time) * df['post']
result = smf.ols('Y ~ time + post + time_since', data=df).fit(cov_type='HAC', cov_kwds={'maxlags':3})

The Frisch-Waugh-Lovell (FWL) theorem states that in a multiple regression Y = X1*beta1 + X2*beta2 + e, the OLS estimate of beta1 is identical to the OLS estimate from regressing the residuals of Y on X2 (M_X2 Y) onto the residuals of X1 on X2 (M_X2 X1), where M_X2 is the annihilator matrix. This is the 'partialing out' or 'residualing out' interpretation: the coefficient on X1 captures the variation in Y unexplained by X2, explained by variation in X1 unexplained by X2. The R-learner and Double ML are direct applications of FWL: residualizing both the outcome and treatment on confounders before regressing residualized outcome on residualized treatment.

import numpy as np
from sklearn.linear_model import LinearRegression
# FWL: partial out X2 from both Y and X1
resid_y = y - LinearRegression().fit(X2, y).predict(X2)
resid_x1 = x1 - LinearRegression().fit(X2, x1).predict(X2)
beta1_fwl = LinearRegression().fit(resid_x1.reshape(-1,1), resid_y).coef_[0]