Causal Inference

Dr. Lucy D’Agostino McGowan

What is Causal Inference?

Statistics tells us: Variables are associated

Causal inference asks: Does changing one variable cause changes in another?

Moving from “X is correlated with Y” to “X causes Y”

Causal DAGs

What is a DAG?

DAG: Directed Acyclic Graph

A visual tool for representing our assumptions about the data generating mechanism

A Simple DAG

This DAG says: X has a causal effect on Y

Example

Example: Ice cream sales and drowning deaths are correlated

Does ice cream cause drowning?

No! Temperature causes both:

DAGs help us think clearly about what causes what

The Fundamental Problem of Causal Inference

Potential Outcomes

For each person \(i\), there are two potential outcomes:

  • \(Y_i(1)\) = outcome if person \(i\) receives an exposure
  • \(Y_i(0)\) = outcome if person \(i\) does not receive exposure

Individual treatment effect: \[\tau_i = Y_i(1) - Y_i(0)\]

The Fundamental Problem

We can only observe ONE potential outcome for each person!

If person \(i\) gets exposed, we see \(Y_i(1)\) but not \(Y_i(0)\)

If person \(i\) doesn’t get exposed, we see \(Y_i(0)\) but not \(Y_i(1)\)

The missing outcome is called the counterfactual

Average Treatment Effect (ATE)

Since we can’t measure individual effects, we focus on averages:

\[\text{ATE} = E[Y(1) - Y(0)] = E[Y(1)] - E[Y(0)]\]

In words: The average difference in outcomes if everyone was exposed vs. if no one was exposed

Naive Comparison Doesn’t Work

Tempting approach: Compare exposed vs. unexposed groups \[E[Y|X=1] - E[Y|X=0]\]

Problem: This is often NOT the same as the ATE because exposed and unexposed groups might differ in other ways

Example: The Drug Study

Imagine a drug study where:

  • Older people are more likely to receive the drug
  • Older people tend to have worse health outcomes
  • The drug actually helps (true effect = +10 points)
# Naive comparison
naive_estimate <- mean(data$outcome[data$treatment == 1]) - 
                  mean(data$outcome[data$treatment == 0])
naive_estimate
[1] 5.212386

The naive estimate is wrong! It’s confounded by age

Confounders

What is a Confounder?

A confounder is a variable that:

  1. Affects the treatment/exposure
  2. Affects the outcome
  3. Is not caused by the treatment

Why Confounders Matter

Confounders create spurious associations

The correlation between X and Y is not purely causal, some of it is due to C

Solution: We need to “control for” or “adjust for” the confounder

Multiple Confounders

Real problems often have multiple confounders

We need to adjust for all of them to get the causal effect

Colliders

What is a Collider?

A collider is a variable that is caused by two other variables

Key point: Arrows point into the collider

Why Colliders Are Dangerous

Do NOT control for colliders!

Controlling for a collider creates a spurious association between its causes

This is called collider bias or selection bias

Example: The Talent Paradox

Among people who were hired (conditioning on the collider):

  • High talent people may have fewer connections
  • Low talent people may have more connections
  • This creates a negative correlation that doesn’t exist in the full population!

Visualizing Collider Bias

This is collider bias in action!

Confounder vs. Collider

Rule of thumb:

  • Arrows out → confounder → control for it
  • Arrows in → collider → don’t control for it

Paths and Backdoor Paths

Paths in a DAG

A path is any route from X to Y following the arrows (in any direction)

Path 1: X → Y (direct causal path)

Path 2: X ← C → Y (path through confounder)

Open vs. Closed Paths

Open path: Association can flow through it

Closed path: Association is blocked

A path is closed if:

  1. There’s a collider on the path (and we don’t condition on it), OR
  2. We condition on a non-collider on the path

Backdoor Paths

A backdoor path is a path from X to Y that starts with an arrow into X

Why Backdoor Paths Matter

Goal of causal inference: Isolate the causal effect of X on Y

Problem: Open backdoor paths allow non-causal association to flow

Solution: Close all backdoor paths by conditioning on the right variables

Finding All Paths

library(ggdag)

# Define the DAG
dag <- dagify(
  C1 ~ C2,
  X ~ C1 + C2,
  M ~ C2 + X,
  Y ~ C1 + X + M,
  exposure = "X",
  outcome = "Y"
)

# Find all paths from X to Y
ggdag_paths(dag, from = "X", to = "Y") +
  theme_dag()

Adjustment Sets

What is an Adjustment Set?

An adjustment set is a set of variables that, if we condition on them (control for them), will:

  1. Close all backdoor paths from X to Y
  2. Not open any new biasing paths

In other words: Variables that give us the causal effect when we include them in a regression

Finding Adjustment Sets

For this simple DAG, we need to adjust for C

Using ggdag + dagitty to Find Adjustment Sets

library(dagitty)

# Simple confounder DAG
simple_dag <- dagify(
  X ~ C,
  Y ~ C + X,
  exposure = "X",
  outcome = "Y"
)
# Find adjustment sets
adjustmentSets(simple_dag, exposure = "X", outcome = "Y")
{ C }

Result: We need to adjust for C

The Backdoor Criterion

A set of variables satisfies the backdoor criterion if:

  1. No variable in the set is a descendant of X (not caused by X)
  2. The set blocks all backdoor paths from X to Y

Practical meaning: These are valid adjustment sets

Example: What NOT to Control For

Don’t control for mediators! They’re part of the causal effect we want to measure

Multiple Valid Adjustment Sets

# More complex DAG
complex_dag <- dagify(
  X ~ C1 + C2,
  Y ~ C1 + C2 + C3 + X,
  C3 ~ C2,
  exposure = "X",
  outcome = "Y"
)

adjustmentSets(complex_dag, exposure = "X", outcome = "Y", type = "all")
{ C1, C2 }
{ C1, C2, C3 }

Both work! Usually we prefer smaller sets (fewer variables to measure)

The Minimally Sufficient Adjustment Set

# Same complex DAG
adjustmentSets(complex_dag, exposure = "X", outcome = "Y", 
               type = "minimal")
{ C1, C2 }

Minimal set: The smallest set that closes all backdoor paths

Visualizing Adjustment Sets

From DAGs to Regression

The Adjustment Formula

Once we identify an adjustment set Z, we can estimate the ATE:

\[E[Y(1) - Y(0)] =\\ E_Z[E[Y|X=1, Z] - E[Y|X=0, Z]]\]

Proving the Adjustment Formula

What We Want to Prove

Goal: Show that adjusting for confounders gives us the ATE

\[E[Y(1) - Y(0)] = E_Z[E[Y|X=1, Z] - E[Y|X=0, Z]]\]

In words: The average treatment effect equals the average difference in outcomes when we condition on the confounders

Key Assumptions

1. Conditional Exchangeability (No unmeasured confounding)

\[Y(1), Y(0) \perp\!\!\!\perp X \mid Z\]

Given Z, exposure assignment is “as if random”

Key Assumptions

2. Positivity

\[0 < P(X=1|Z) < 1 \text{ for all } Z\]

Everyone has some chance of being exposed and unexposed

Key Assumptions

3. Consistency

\[Y = Y(1) \cdot X + Y(0) \cdot (1-X)\]

The outcome we observe equals the potential outcome under the exposure we received

The Proof: Step 1

Start with the definition of ATE:

\[E[Y(1) - Y(0)] = E[Y(1)] - E[Y(0)]\]

Focus on \(E[Y(1)]\) first (same logic applies to \(E[Y(0)]\))

The Proof: Step 2

Use the law of iterated expectations:

\[E[Y(1)] = E_Z[E[Y(1)|Z]]\]

The Proof: Step 3

Now use conditional exchangeability:

\[E[Y(1)|Z] = E[Y(1)|X=1, Z]\]

Why? Because \(Y(1) \perp\!\!\!\perp X \mid Z\)

Given Z, the potential outcome \(Y(1)\) is independent of whether you actually were exposed

Intuition: Among people with the same Z, those who were exposed are no different (in terms of potential outcomes) from those who weren’t

The Proof: Step 4

Use consistency:

\[E[Y(1)|X=1, Z] = E[Y|X=1, Z]\]

Why? For people who actually were exposed (\(X=1\)), their observed outcome \(Y\) equals their potential outcome under exposure \(Y(1)\)

This is the crucial step: We’ve gone from unobservable potential outcomes to observable actual outcomes!

The Proof: Putting It Together

Combining Steps 2-4:

\[E[Y(1)] = E_Z[E[Y(1)|Z]]\] \[= E_Z[E[Y(1)|X=1, Z]]\] \[= E_Z[E[Y|X=1, Z]]\]

By the same logic:

\[E[Y(0)] = E_Z[E[Y|X=0, Z]]\]

The Proof: Final Result

Therefore:

\[E[Y(1) - Y(0)] = E[Y(1)] - E[Y(0)]\] \[= E_Z[E[Y|X=1, Z]] - E_Z[E[Y|X=0, Z]]\] \[= E_Z[E[Y|X=1, Z] - E[Y|X=0, Z]]\]

From Theory to Practice

The Practical Question

How do we actually compute \(E[Y|X, Z]\)?

If Z is continuous or multidimensional, we can’t just stratify and take group means

Some Approaches

Parametric modeling: Assume \(E[Y|X, Z] = \beta_0 + \beta_1 X + \boldsymbol\beta^T \mathbf{Z}\) - Fit with linear regression - Use model to predict under \(X=1\) and \(X=0\)

Matching: Pair treated and control units with similar Z values

Weighting: Reweight the sample to balance Z across treatment groups

Simple Linear Regression Approach

If we assume a linear model: \[Y = \beta_0 + \beta_1 X + \boldsymbol\beta^T \mathbf{Z} + \varepsilon\]

Under (strong!) conditions: \[\hat{\beta}_1 = \text{Average Treatment Effect}\]

Assumptions for Linear Regression = ATE

For \(\hat{\beta}_1\) to equal the ATE, we need:

  1. Correct adjustment set: Z satisfies the backdoor criterion
  2. No unmeasured confounding: We measured all confounders
  3. Linearity: Effects are additive (no interactions)
  4. Positivity: Every unit has non-zero probability of both exposures

Example: Simulated Data

set.seed(1)
n <- 1000

# Confounder
C <- rnorm(n)

# Exposure (affected by confounder)
X <- rbinom(n, 1, plogis(0.5 * C))

# Outcome (true ATE = 2)
Y <- 2 * X + 3 * C + rnorm(n)

Naive Regression (Confounded)

# Without adjusting for C
naive_model <- lm(Y ~ X)
coef(naive_model)["X"]
      X 
3.70334

Wrong! This is biased because C affects both X and Y

Adjusted Regression (Correct)

# Adjusting for C
adjusted_model <- lm(Y ~ X + C)
coef(adjusted_model)["X"]
       X 
2.029299

Correct! This recovers the true ATE of 2 (plus sampling error)

Why This Works: FWL Perspective

Remember the Frisch-Waugh-Lovell theorem?

When we regress Y ~ X + C, the coefficient on X is:

The relationship between the part of Y that C can’t explain and the part of X that C can’t explain

This is exactly what we want! We’ve removed C’s confounding influence

Visualizing the Adjustment

G-Computation for Linear Regression

What is G-Computation?

G-computation is a method to estimate causal effects by:

  1. Fitting a model for the outcome
  2. Predicting outcomes under different treatment values
  3. Averaging these predictions

The G-Computation Formula

Recall the adjustment formula: \[E[Y(1)] = E_Z[E[Y|X=1, Z]]\]

The G-Computation Formula

G-computation procedure:

  1. Fit \(E[Y|X, Z]\) (our outcome model)
  2. For each unit, predict \(\hat{Y}_i(1)\) setting \(X=1\), using their \(Z_i\)
  3. Average: \(\frac{1}{n}\sum_i \hat{Y}_i(1)\)
  4. Repeat for \(X=0\) to get \(\hat{E}[Y(0)]\)
  5. ATE = \(\hat{E}[Y(1)] - \hat{E}[Y(0)]\)

Why Bother with G-Computation?

For simple linear models with no interactions: The regression coefficient equals the ATE

G-computation becomes useful when:

  • You have interactions between treatment and covariates
  • You have non-linear models
  • You want to estimate effects on different scales

G-Computation: Step by Step

# Use our previous example data
set.seed(1)
n <- 1000
C <- rnorm(n, 1)
X <- rbinom(n, 1, plogis(0.5 * C))
Y <- 2 * X + 3 * C + rnorm(n)
data_all <- data.frame(X = X, C = C, Y = Y)

# Step 1: Fit the outcome model
outcome_model <- lm(Y ~ X + C, data_all)

G-Computation: Step 2-3

# Step 2: Create datasets with everyone treated and untreated

# Everyone treated (X = 1)
data_treated <- data_all
data_treated$X <- 1

# Everyone untreated (X = 0)
data_control <- data_all
data_control$X <- 0

G-Computation: Step 4-5

# Step 3: Predict outcomes under both scenarios
Y1_pred <- predict(outcome_model, newdata = data_treated)
Y0_pred <- predict(outcome_model, newdata = data_control)

# Step 4: Average predictions
E_Y1 <- mean(Y1_pred)
E_Y0 <- mean(Y0_pred)

# Step 5: Compute ATE
ATE_gcomp <- E_Y1 - E_Y0
ATE_gcomp
[1] 2.025261

Same as the regression coefficient! (for linear models with no interactions)

When G-Computation Differs: Interactions

set.seed(1)
n <- 1000
C <- rnorm(n, 1)
X <- rbinom(n, 1, plogis(0.25 * C))

# Treatment effect varies by C! (interaction)
Y_interact <- 2 * X + 3 * C + 1.5 * X * C + rnorm(n)

Now the treatment effect is not constant:

  • When C = 0, effect is 2
  • When C = 1, effect is 2 + 1.5 = 3.5

Interaction Model: Regression Approach

# Fit model with interaction
interact_model <- lm(Y_interact ~ X + C + X:C)
coef(interact_model)
(Intercept)           X           C         X:C 
-0.09771085  2.14169352  3.04054859  1.44840785

Question: What’s the average treatment effect?

Problem: The coefficient on X (2.14) is the effect when C = 0, not the average!

Interaction Model: G-Computation

# Create counterfactual datasets
data_all_int <- data.frame(X = X, C = C, Y = Y_interact)
data_treated_int <- data_all_int
data_treated_int$X <- 1
data_control_int <- data_all_int
data_control_int$X <- 0

# Predict and average
Y1_pred_int <- predict(interact_model, newdata = data_treated_int)
Y0_pred_int <- predict(interact_model, newdata = data_control_int)

ATE_gcomp_int <- mean(Y1_pred_int) - mean(Y0_pred_int)
ATE_gcomp_int
[1] 3.57323

G-computation gives us the average effect across the distribution of C!

Why G-Computation Works Here

The average treatment effect is: \(\text{ATE} = E[2 + 1.5 \times C]\)

Since \(E[C] \approx 1\), the ATE \(\approx 3.5\)

G-computation automatically:

  1. Computes effect for each value of C
  2. Averages using the actual distribution of C

Regression coefficient: Only gives effect at C = 0

Visualizing Heterogeneous Effects

G-computation averages the blue line using the distribution on the right

G-Computation Algorithm Summary

For any outcome model:

  1. Fit model: \(E[Y|X, Z]\)
  2. Create \(n\) copies of data with \(X = 1\)
  3. Predict: \(\hat{Y}_i(1)\) for each unit
  4. Create \(n\) copies of data with \(X = 0\)
  5. Predict: \(\hat{Y}_i(0)\) for each unit
  6. Compute: \(\hat{\text{ATE}} = \frac{1}{n}\sum_i[\hat{Y}_i(1) - \hat{Y}_i(0)]\)

Putting It All Together

The Complete Workflow

Step 1: Draw the DAG

  • List all relevant variables
  • Draw arrows based on domain knowledge
  • Identify treatment (X) and outcome (Y)

The Complete Workflow

Step 2: Identify adjustment set

  • Find all backdoor paths
  • Determine minimal sufficient adjustment set
  • Check for colliders (don’t adjust!)

The Complete Workflow

Step 3: Collect data

  • Measure treatment, outcome, and adjustment set
  • More data is better (reduces sampling variability)

Step 4: Fit the model

  • Start with linear regression: Y ~ X + Z
  • Consider interactions if theory suggests heterogeneity

The Complete Workflow

Step 5: Estimate the ATE

  • Simple model: Use coefficient
  • Complex model: Use G-computation

Key Assumptions to Remember

1. No unmeasured confounding

  • We measured all confounders
  • Requires domain knowledge

2. Positivity

  • Everyone has some chance of treatment and control
  • Check for overlap in covariate distributions

3. Correct model specification

  • For linear regression: linearity, include important interactions

Common Pitfalls

❌ Controlling for colliders

Creates bias where there was none!

❌ Not controlling for confounders

Leads to biased estimates

❌ Controlling for mediators

Blocks the causal effect you want to measure

❌ Interpreting regression coefficients as causal without thinking about the DAG

Correlation ≠ causation!

Tips for Success

✓ Start with the DAG

Think causally before looking at data

✓ Use domain knowledge

Statistics can’t tell you what causes what

✓ Be transparent about assumptions

No unmeasured confounding is a strong assumption

✓ Consider sensitivity analyses

What if there’s an unmeasured confounder?

Sensitivity to Unmeasured Confounding

The Problem

What if we missed a confounder?

Our estimate will be biased if U affects both X and Y

Omitted Variable Bias Formula

What we estimate when U is unmeasured:

\[\hat{\beta}_X = \beta_X + \beta_U \times \delta_{UX}\]

Example: Exercise and Heart Health

set.seed(1)
n <- 1000

# Unmeasured confounder: genetic fitness
genetic_fitness <- rnorm(n)

# Measured confounder: age
age <- rnorm(n, mean = 50, sd = 10)

# Exercise depends on both age and genetics
exercise_hours <- 5 + (-0.08 * age) + 0.6 * genetic_fitness + rnorm(n, sd = 1)
exercise_hours <- pmax(exercise_hours, 0)  # Can't be negative

# Heart health depends on exercise AND genetics
# True effect of exercise = -2 (negative = better)
heart_disease_risk <- 60 + (-2 * exercise_hours) + 
                      0.3 * age + (-3 * genetic_fitness) + rnorm(n, sd = 3)

Naive Estimate (No Controls)

# No controls
naive <- lm(heart_disease_risk ~ exercise_hours)
coef(naive)["exercise_hours"]
exercise_hours 
     -4.503693

Biased toward zero (or even positive): People with good genetics exercise more AND have lower risk

Controlling for Measured Confounder

# Control for age only
partial_control <- lm(heart_disease_risk ~ exercise_hours + age)
coef(partial_control)["exercise_hours"]
exercise_hours 
     -3.518357

Still biased! We adjusted for age, but genetic fitness remains unmeasured

Oracle Estimate (If We Could Measure Everything)

# Control for both (hypothetical - we can't measure genetics)
full_control <- lm(heart_disease_risk ~ exercise_hours + age + genetic_fitness)
coef(full_control)["exercise_hours"]
exercise_hours 
     -2.008628

Unbiased: Recovers the true effect of -2 (plus sampling error)

Computing the Bias Using FWL

# Our estimate from partial control
beta_x_biased <- coef(partial_control)["exercise_hours"]

# Regress outcome on unmeasured U (controlling for measured)
beta_u <- coef(lm(heart_disease_risk ~ genetic_fitness + age))["genetic_fitness"]

# Regress treatment on unmeasured U (controlling for measured)  
delta_ux <- coef(lm(exercise_hours ~ genetic_fitness + age))["genetic_fitness"]

# Predicted bias
predicted_bias <- beta_u * delta_ux
predicted_bias
genetic_fitness 
      -1.884238

Predicted vs. Actual Bias

# Actual bias
true_effect <- -2
actual_bias <- beta_x_biased - true_effect

data.frame(
  actual_bias = actual_bias,
  predicted_bias = predicted_bias
)
               actual_bias predicted_bias
exercise_hours   -1.518357      -1.884238

They match! The OVB formula works