# install.packages(c("tidyverse", "dagitty", "ggdag"))
library(tidyverse)
library(dagitty)
library(ggdag)

The big idea

A DAG (directed acyclic graph) is a picture of our causal assumptions. Each arrow says “this variable directly influences that one.”

To check that the DAG does what we claim, we do this each time:

  1. Draw the DAG.
  2. Simulate data from it. Each variable is drawn with rnorm(), and its mean depends on the variables pointing into it through coefficients we choose (the “true parameters”).
  3. Fit a regression and see whether the association between X and Y shows up or disappears when we condition on (add to the model) a third variable Z.

Because we picked the true coefficients, we know the right answer, so we can see exactly when a model is biased.

set.seed(2026)
n <- 2000   # number of simulated observations

1. The Fork: X <- Z -> Y (common cause / confounder)

Z causes both X and Y. X does not cause Y, so the true effect of X on Y is 0.

fork_dag <- dagify(
  X ~ Z,
  Y ~ Z,
  coords = list(x = c(X = 0, Z = 1, Y = 2),
                y = c(X = 0, Z = 1, Y = 0))
)

ggdag(fork_dag, text_size = 5) + theme_dag()

Simulate (true parameters: Z -> X is 1.5, Z -> Y is 2, X -> Y is 0):

fork_data <- tibble(
  Z = rnorm(n, mean = 0, sd = 1),
  X = rnorm(n, mean = 1.5 * Z, sd = 1),
  Y = rnorm(n, mean = 2 * Z, sd = 1)
)
# Not conditioning on Z: a spurious association appears
fork_no_z <- lm(Y ~ X, data = fork_data)
# Conditioning on Z: the association should vanish
fork_with_z <- lm(Y ~ X + Z, data = fork_data)

coef(summary(fork_no_z))["X", c("Estimate", "Std. Error")]
##   Estimate Std. Error 
## 0.89335482 0.01967857
coef(summary(fork_with_z))["X", c("Estimate", "Std. Error")]
##    Estimate  Std. Error 
## -0.04744988  0.02277931

Without Z, X looks like it affects Y (the estimate is far from 0), even though the true effect is 0. With Z in the model, the estimate is about 0.

ggplot(fork_data, aes(x = X, y = Y)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Fork: X and Y are associated (through Z)")

D-separation check. dagitty can tell us this without any simulation:

dseparated(fork_dag, "X", "Y")                        # FALSE: path is open
## [1] FALSE
dseparated(fork_dag, "X", "Y", "Z")                   # TRUE: conditioning on Z closes it
## [1] TRUE
impliedConditionalIndependencies(fork_dag)            # X _||_ Y | Z
## X _||_ Y | Z
ggdag_dseparated(fork_dag, from = "X", to = "Y", controlling_for = "Z", text_size = 5) +
  theme_dag()

2. The Pipe (chain): X -> Z -> Y (mediator)

X causes Z, and Z causes Y. Z carries the effect of X to Y.

pipe_dag <- dagify(
  Z ~ X,
  Y ~ Z,
  coords = list(x = c(X = 0, Z = 1, Y = 2),
                y = c(X = 0, Z = 0, Y = 0))
)

ggdag(pipe_dag, text_size = 5) + theme_dag()

Simulate (X -> Z is 1.5, Z -> Y is 2, so the total effect of X on Y is 1.5 x 2 = 3):

pipe_data <- tibble(
  X = rnorm(n, mean = 0, sd = 1),
  Z = rnorm(n, mean = 1.5 * X, sd = 1),
  Y = rnorm(n, mean = 2 * Z, sd = 1)
)
# Not conditioning on Z: we recover the TOTAL effect (about 3)
pipe_no_z <- lm(Y ~ X, data = pipe_data)
# Conditioning on Z: we block the path, effect goes to about 0
pipe_with_z <- lm(Y ~ X + Z, data = pipe_data)

coef(summary(pipe_no_z))["X", c("Estimate", "Std. Error")]
##   Estimate Std. Error 
## 3.07211226 0.04952265
coef(summary(pipe_with_z))["X", c("Estimate", "Std. Error")]
##   Estimate Std. Error 
## 0.06037674 0.04093205

Conditioning on a mediator removes the association, but this is not always what we want. If we care about the total effect of X on Y, adding Z is a mistake (“overcontrol bias” or “post-treatment bias”).

dseparated(pipe_dag, "X", "Y")
## [1] FALSE
dseparated(pipe_dag, "X", "Y", "Z")
## [1] TRUE

3. The Collider: X -> Z <- Y

Both X and Y cause Z. X and Y are independent of each other. Here the path is closed by default, and conditioning on Z opens it.

collider_dag <- dagify(
  Z ~ X + Y,
  coords = list(x = c(X = 0, Z = 1, Y = 2),
                y = c(X = 1, Z = 0, Y = 1))
)

ggdag(collider_dag, text_size = 5) + theme_dag()

Simulate (X and Y are unrelated; Z = 1.5 X + 2 Y + noise):

collider_data <- tibble(
  X = rnorm(n, mean = 0, sd = 1),
  Y = rnorm(n, mean = 0, sd = 1),
  Z = rnorm(n, mean = 1.5 * X + 2 * Y, sd = 1)
)
# Not conditioning on Z: no association (correct!)
collider_no_z <- lm(Y ~ X, data = collider_data)
# Conditioning on Z: a FAKE association appears (collider bias)
collider_with_z <- lm(Y ~ X + Z, data = collider_data)

coef(summary(collider_no_z))["X", c("Estimate", "Std. Error")]
##    Estimate  Std. Error 
## 0.009100393 0.022700026
coef(summary(collider_with_z))["X", c("Estimate", "Std. Error")]
##    Estimate  Std. Error 
## -0.58903175  0.01197979

The estimate with Z in the model is negative even though the true effect is 0. A common intuition: if Z is high and X is low, Y must be high to make up for it.

# Visualize: slice Z into groups. Within each group, X and Y become negatively related.
collider_data %>%
  mutate(Z_group = cut(Z, breaks = quantile(Z, c(0, 1/3, 2/3, 1)),
                       labels = c("Low Z", "Mid Z", "High Z"),
                       include.lowest = TRUE)) %>%
  ggplot(aes(x = X, y = Y, color = Z_group)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Collider: conditioning on Z creates an X-Y association")

dseparated(collider_dag, "X", "Y")        # TRUE: closed
## [1] TRUE
dseparated(collider_dag, "X", "Y", "Z")   # FALSE: conditioning opened it
## [1] FALSE

4. The Descendant: conditioning on a child of a collider

If D is caused by the collider Z, then conditioning on D partly conditions on Z, and opens the path too.

descendant_dag <- dagify(
  Z ~ X + Y,
  D ~ Z,
  coords = list(x = c(X = 0, Z = 1, Y = 2, D = 1),
                y = c(X = 1, Z = 0, Y = 1, D = -1))
)

ggdag(descendant_dag, text_size = 5) + theme_dag()

descendant_data <- tibble(
  X = rnorm(n, mean = 0, sd = 1),
  Y = rnorm(n, mean = 0, sd = 1),
  Z = rnorm(n, mean = 1.5 * X + 2 * Y, sd = 1),
  D = rnorm(n, mean = 1 * Z, sd = 1)
)
descendant_no_d <- lm(Y ~ X, data = descendant_data)
descendant_with_d <- lm(Y ~ X + D, data = descendant_data)

coef(summary(descendant_no_d))["X", c("Estimate", "Std. Error")]
##   Estimate Std. Error 
## 0.01928239 0.02241204
coef(summary(descendant_with_d))["X", c("Estimate", "Std. Error")]
##    Estimate  Std. Error 
## -0.48247287  0.01533924
dseparated(descendant_dag, "X", "Y")
## [1] TRUE
dseparated(descendant_dag, "X", "Y", "D")
## [1] FALSE

5. Putting it together: a bias example with a real effect

Now let’s use a DAG with a true effect of X on Y and a confounder Z.

  • Z -> X (0.8), Z -> Y (1.2), and X -> Y (0.5, the true causal effect)
bias_dag <- dagify(
  X ~ Z,
  Y ~ X + Z,
  exposure = "X",
  outcome = "Y",
  coords = list(x = c(X = 0, Z = 1, Y = 2),
                y = c(X = 0, Z = 1, Y = 0))
)

ggdag(bias_dag, text_size = 5) + theme_dag()

dagitty can tell us what to adjust for to estimate the effect of X on Y:

adjustmentSets(bias_dag)     # { Z }
## { Z }
ggdag_adjustment_set(bias_dag, text_size = 4) + theme_dag()

bias_data <- tibble(
  Z = rnorm(n, mean = 0, sd = 1),
  X = rnorm(n, mean = 0.8 * Z, sd = 1),
  Y = rnorm(n, mean = 0.5 * X + 1.2 * Z, sd = 1)
)
bias_naive <- lm(Y ~ X, data = bias_data)       # biased: confounded by Z
bias_adjusted <- lm(Y ~ X + Z, data = bias_data) # adjusted: should be about 0.5

results <- tibble(
  model = c("Naive: Y ~ X", "Adjusted: Y ~ X + Z", "Truth"),
  estimate = c(coef(bias_naive)["X"], coef(bias_adjusted)["X"], 0.5)
)
results
ggplot(results, aes(x = model, y = estimate, fill = model)) +
  geom_col(show.legend = FALSE) +
  geom_hline(yintercept = 0.5, linetype = "dashed") +
  labs(title = "Confounding bias: effect of X on Y",
       x = NULL, y = "Estimated effect of X")

The naive model overstates the effect of X, because part of the X-Y association is just Z driving both. Conditioning on Z fixes it.

Summary

Structure DAG X and Y associated (unconditioned)? Condition on Z What to do
Fork X <- Z -> Y Yes (spurious) Removes association Condition on Z (confounder)
Pipe X -> Z -> Y Yes (causal, via Z) Removes association Don’t condition if you want the total effect
Collider X -> Z <- Y No Creates an association Don’t condition on Z
Descendant collider with D below Z No Creates a (weaker) association Don’t condition on D

Rule of thumb: whether to add a variable to a model depends on where it sits in the DAG, not on whether it improves fit or is “significant.”