# Load necessary libraries
library(tidyverse)
library(stargazer)

# Data-generating process (DGP)
set.seed(123)
tb <- tibble(
  female = ifelse(runif(10000) >= 0.5, 1, 0),
  ability = rnorm(10000),
  discrimination = female,  # Discrimination linked to being female
  occupation = 1 + 2 * ability - 2 * discrimination + rnorm(10000),
  wage = 1 - discrimination + occupation + 2 * ability + rnorm(10000)
)

# Regressions
lm_1 <- lm(wage ~ female, tb)                          # Unconditional, biased
lm_2 <- lm(wage ~ female + occupation, tb)             # Conditional, still biased
lm_3 <- lm(wage ~ female + occupation + ability, tb)   # Fully conditional, unbiased

# Display results with stargazer
stargazer(lm_1, lm_2, lm_3, type = "text",
          column.labels = c("Unconditional (Biased)", 
                            "Conditional (Biased)", 
                            "Fully Conditional (Unbiased)"),
          dep.var.caption = "Dependent variable: Wage",
          omit.stat = c("f", "ser"))