Quantile Regression

# Quantile Regression: Income returns to education across distribution
# Demonstrates heterogeneous effects that OLS misses

library(quantreg)    # For quantile regression
library(ggplot2)     # For visualization
library(dplyr)

# Simulate realistic income data (right-skewed, heterogeneous education effects)
set.seed(2025)
n <- 500
data <- data.frame(
  education = rnorm(n, 14, 3)  # Years of education, M=14, SD=3
)
data$education <- pmax(8, pmin(22, data$education))  # Bound 8-22 years

# Experience (correlated with education)
data$experience <- pmax(0, rnorm(n, 15 - 0.5*data$education + 20, 5))

# Income: heterogeneous returns (larger effect for high earners)
data$income <- 15 + 
  2.5 * data$education +              # Base education effect
  0.3 * data$experience +             # Experience effect
  0.15 * data$education * data$experience/10 +  # Interaction
  exp(rnorm(n, 0.8, 0.6))            # Right-skewed errors

# === STEP 1: Compare OLS vs Quantile Regression ===

# OLS regression (standard approach)
ols_model <- lm(income ~ education + experience, data = data)
summary(ols_model)

# Median regression (τ = 0.5)
median_model <- rq(income ~ education + experience, tau = 0.5, data = data)
summary(median_model, se = "boot")  # Bootstrap SE for inference

# === STEP 2: Quantile Regression at Multiple Quantiles ===

taus <- c(0.1, 0.25, 0.5, 0.75, 0.9)
qr_models <- rq(income ~ education + experience, tau = taus, data = data)
summary(qr_models, se = "boot", R = 1000)  # 1000 bootstrap replicates

# Extract coefficients
coefs <- coef(qr_models)
print(round(coefs, 3))

# === STEP 3: Visualize Heterogeneous Effects ===

# Plot education coefficient across quantiles
qr_full <- rq(income ~ education + experience, tau = seq(0.05, 0.95, 0.05), data = data)
qr_summary <- summary(qr_full, se = "boot", R = 500)

# Extract education coefficients and CIs
edu_coefs <- sapply(qr_summary, function(x) x$coefficients["education", "Value"])
edu_lower <- sapply(qr_summary, function(x) x$coefficients["education", "lower bd"])
edu_upper <- sapply(qr_summary, function(x) x$coefficients["education", "upper bd"])
taus_full <- seq(0.05, 0.95, 0.05)

plot_data <- data.frame(tau = taus_full, coef = edu_coefs, 
                        lower = edu_lower, upper = edu_upper)

ggplot(plot_data, aes(x = tau, y = coef)) +
  geom_line(color = "blue", size = 1) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = "blue") +
  geom_hline(yintercept = coef(ols_model)["education"], 
             linetype = "dashed", color = "red", size = 1) +
  labs(title = "Education Returns Across Income Distribution",
       subtitle = "Blue: Quantile regression coefficients | Red: OLS coefficient",
       x = "Income Quantile(τ)",
       y = "Education Coefficient($/year)") +
  theme_classic() +
  annotate("text", x = 0.5, y = coef(ols_model)["education"] + 0.5, 
           label = "OLS(constant)", color = "red")

# === STEP 4: Check Quantile Crossing ===

# Predict quantiles at different education levels
edu_range <- seq(10, 20, by = 2)
pred_data <- expand.grid(education = edu_range, tau = taus)
pred_data$experience <- mean(data$experience)  # Hold experience at mean

# Get predictions
predictions <- NULL
for (tau in taus) {
  temp <- data.frame(
    education = edu_range,
    experience = mean(data$experience)
  )
  temp$income_pred <- predict(rq(income ~ education + experience, tau = tau, data = data),
                              newdata = temp)
  temp$tau <- tau
  predictions <- rbind(predictions, temp)
}

ggplot(predictions, aes(x = education, y = income_pred, color = factor(tau), group = tau)) +
  geom_line(size = 1) +
  labs(title = "Predicted Income Quantiles by Education",
       subtitle = "Lines should not cross(monotonicity check)",
       x = "Years of Education",
       y = "Predicted Income($1000s)",
       color = "Quantile(τ)") +
  scale_color_brewer(palette = "RdYlBu") +
  theme_classic()

# === STEP 5: Test for Coefficient Equality Across Quantiles ===

# Wald test: Are education effects equal across quantiles?
anovatest <- anova(qr_models, test = "Wald", joint = FALSE)
print(anovatest)

# === APA-Style Reporting ===
cat("\n=== Results Summary ===\n")
cat("Quantile regression revealed heterogeneous education returns across the\n")
cat("income distribution. At the 10th percentile, each additional year of\n")
cat("education was associated with $", round(coefs["education", "tau= 0.10"], 2), 
    "k higher income.\n")
cat("At the 90th percentile, the return increased to $", 
    round(coefs["education", "tau= 0.90"], 2), "k per year(95% CI [X, X]).\n")
cat("Wald tests confirmed coefficients differed significantly across quantiles\n")
cat(", highlighting distributional heterogeneity missed by OLS(constant\n")
cat("effect = $", round(coef(ols_model)["education"], 2), "k).\n")