# Load necessary library
library(tidyverse)
library(broom)

# Set seed and generate data
set.seed(1000)
n <- 10000

# Create the data
data <- tibble(
  background = rnorm(n),
  pe = 1 + 1 * background + rnorm(n),
  income = 1 + 5 * pe + rnorm(n),
  college = 1 + 1 * background + 1 * pe + 1 * income + rnorm(n),
  earnings = 1 + 2 * college + 1 * income + rnorm(n)
)

# Treatment effect of college on earnings is 2

# Naive model is biased
naive_model <- lm(earnings ~ college, data = data)
summary(naive_model)

# Adjusted model controlling for income to satisfy the backdoor criterion
adjusted_model <- lm(earnings ~ college + income, data = data)
summary(adjusted_model)