1. Modeling HW with Data from the Current Population Survey (CPS)

This note is a Bayesian hourly-wage regression with weakly informative priors elicited via quantile-parameterized distributions (GLD / J-QPD helpers in Stan). The directed structure is the simple regression sketched below: log real wages depend on gender, race, marital status, and age (in decades), plus an intercept and Gaussian noise.

Hand-drawn model diagram from the original write-up.

Hand-drawn model diagram from the original write-up.

Data: CEPR CPS ORG Uniform Extracts (March 2019; see data/README.md). Variable names follow the CEPR conventions used in the original project (rw, female, wbho, married, age). Exact vintage from the 2020 HTML is unavailable, so summaries will not match that page digit-for-digit; the elicitation workflow and model class do.

1.1 Prior predictive distribution: a simple model

rstan::expose_stan_functions("stan/quantile_functions.stan")
source("R/GLD_helpers.R")

cps_raw <- readr::read_csv("data/cps_cepr_org_2019_march.csv", show_col_types = FALSE)
summary(cps_raw$rw)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
   1.00   13.25   19.23   25.68   30.77  206.50    6743

NB I am recoding the wbho variable (Race) to be a Boolean (Non-White) 1/0. That recode is already applied in the extract.

cps <- cps_raw |>
  filter(
    !is.na(rw), rw > 0,
    !is.na(female), !is.na(wbho), !is.na(married), !is.na(age_dec)
  )

tibble(
  n_raw = nrow(cps_raw),
  n_complete = nrow(cps),
  mean_rw = mean(cps$rw),
  median_rw = median(cps$rw)
)
# A tibble: 1 × 4
  n_raw n_complete mean_rw median_rw
  <int>      <int>   <dbl>     <dbl>
1 18678      11935    25.7      19.2

Now, using the theoretical model above — where I think that one’s gender, race, marital status, and age (in decades) have a bearing on one’s hourly wages, plus some intercept and an error term — let’s construct some prior predictive distribution…

set.seed(42)
# y = alpha + beta * x
# b1/Female, b2/Race (Non-White), b3/Married, b4/Age^2 (in decades)
# NB Median should be about equidistant between the lower and upper quartiles.

(a_s1 <- GLD_solver(
  lower_quartile = -0.3, median = -0.2, upper_quartile = -0.1,
  other_quantile = 0.2, alpha = 0.9
))
asymmetry steepness 
0.0000000 0.9653298 
(a_s2 <- GLD_solver(
  lower_quartile = -0.01, median = 0, upper_quartile = 0.01,
  other_quantile = 0.03, alpha = 0.9
))
asymmetry steepness 
0.0000000 0.9283763 
(a_s3 <- GLD_solver(
  lower_quartile = 0.3, median = 0.5, upper_quartile = 0.65,
  other_quantile = 0.85, alpha = 0.9
))
 asymmetry  steepness 
-0.4226581  0.9033965 
(a_s4 <- GLD_solver(
  lower_quartile = 0, median = 0.2, upper_quartile = 0.3,
  other_quantile = 0.5, alpha = 0.9
))
 asymmetry  steepness 
-0.6771304  0.9722460

1.2 Drawing from the prior predictive distribution

m_alpha <- log(16)
s_alpha <- 0.1 # Normal prior for alpha near typical median wage.
q_sigma <- c(lower = 0.25, median = 1, upper = 2)

d <- model.matrix(
  log(rw) ~ female + wbho + married + age_dec,
  data = cps
)[, -1]
d <- sweep(d, MARGIN = 2, STATS = colMeans(d), FUN = `-`)

n_draws <- 1000L
ppd <- matrix(NA_real_, nrow = nrow(d), ncol = n_draws)
r_squared <- numeric(n_draws)

for (j in seq_len(n_draws)) {
  alpha_ <- rnorm(1, mean = m_alpha, sd = s_alpha)
  sigma_ <- JQPDS_rng(lower_bound = 0, alpha = 0.25, quantiles = q_sigma)
  beta1_ <- GLD_rng(
    median = -0.2, IQR = -0.1 - (-0.3),
    asymmetry = a_s1[1], steepness = a_s1[2]
  )
  beta2_ <- GLD_rng(
    median = 0, IQR = 0.01 - (-0.01),
    asymmetry = a_s2[1], steepness = a_s2[2]
  )
  beta3_ <- GLD_rng(
    median = 0.5, IQR = 0.65 - 0.3,
    asymmetry = a_s3[1], steepness = a_s3[2]
  )
  beta4_ <- GLD_rng(
    median = 0.2, IQR = 0.3 - 0,
    asymmetry = a_s4[1], steepness = a_s4[2]
  )
  mu_ <- alpha_ +
    beta1_ * d[, "female"] +
    beta2_ * d[, "wbho"] +
    beta3_ * d[, "married"] +
    beta4_ * (d[, "age_dec"]^2)
  epsilon_ <- rnorm(n = length(mu_), mean = 0, sd = sigma_)
  y_ <- mu_ + epsilon_
  ppd[, j] <- y_
  r_squared[j] <- stats::var(mu_) / stats::var(y_)
}

summary(r_squared)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.002021 0.099360 0.457775 0.511951 0.970441 1.000289 
round(quantile(c(exp(ppd)), probs = 1:9 / 10), digits = 3)
    10%     20%     30%     40%     50%     60%     70%     80%     90% 
  0.720   4.607   9.753  14.667  20.039  27.033  39.557  71.673 243.398

1.3 Drawing from the posterior distribution

The 2020 write-up used Gaussian stan_lm on \(\log(\mathrm{rw})\) with an \(R^2\) prior at the prior-predictive median, high adapt_delta, and age entering only as \(\mathrm{age\_dec}^2\). That specification is fit first as the historical baseline. LOO comparison then prefers a richer age mean (\(\mathrm{age\_dec} + \mathrm{age\_dec}^2\)); that richer fit is the preferred model for interpretation below.

r2_loc <- median(r_squared, na.rm = TRUE)

post_quad <- stan_lm(
  log(rw) ~ female + wbho + married + I(age_dec^2),
  data = cps,
  prior_intercept = normal(location = m_alpha, scale = s_alpha, autoscale = FALSE),
  prior = R2(location = r2_loc, what = "median"),
  chains = 4,
  iter = 2000,
  warmup = 1000,
  adapt_delta = 0.999,
  seed = 42,
  refresh = 0
)

post <- stan_lm(
  log(rw) ~ female + wbho + married + age_dec + I(age_dec^2),
  data = cps,
  prior_intercept = normal(location = m_alpha, scale = s_alpha, autoscale = FALSE),
  prior = R2(location = r2_loc, what = "median"),
  chains = 4,
  iter = 2000,
  warmup = 1000,
  adapt_delta = 0.999,
  seed = 43,
  refresh = 0
)

cat("Historical baseline (age_dec^2 only):\n")
Historical baseline (age_dec^2 only):
print(post_quad, digits = 4)
stan_lm
 family:       gaussian [identity]
 formula:      log(rw) ~ female + wbho + married + I(age_dec^2)
 observations: 11935
 predictors:   5
------
             Median  MAD_SD 
(Intercept)   2.8587  0.0128
female       -0.1945  0.0108
wbho         -0.1228  0.0117
married       0.2331  0.0114
I(age_dec^2)  0.0102  0.0005

Auxiliary parameter(s):
              Median MAD_SD
R2            0.1228 0.0053
log-fit_ratio 0.0002 0.0065
sigma         0.5917 0.0039

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg
cat("\nPreferred model (age_dec + age_dec^2):\n")

Preferred model (age_dec + age_dec^2):
print(post, digits = 4)
stan_lm
 family:       gaussian [identity]
 formula:      log(rw) ~ female + wbho + married + age_dec + I(age_dec^2)
 observations: 11935
 predictors:   6
------
             Median  MAD_SD 
(Intercept)   1.7489  0.0508
female       -0.2024  0.0108
wbho         -0.1379  0.0112
married       0.1581  0.0122
age_dec       0.6127  0.0268
I(age_dec^2) -0.0627  0.0032

Auxiliary parameter(s):
              Median  MAD_SD 
R2             0.1545  0.0061
log-fit_ratio -0.0089  0.0067
sigma          0.5756  0.0039

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg
rh <- rhat(post)
ne <- neff_ratio(post)
n_div <- sum(nuts_params(post) |>
  filter(Parameter == "divergent__") |>
  pull(Value))

tibble(
  model = "preferred (age_dec + age_dec^2)",
  max_rhat = max(rh, na.rm = TRUE),
  min_neff_ratio = min(ne, na.rm = TRUE),
  n_divergent = n_div,
  r2_prior_location = r2_loc
)
# A tibble: 1 × 5
  model                    max_rhat min_neff_ratio n_divergent r2_prior_location
  <chr>                       <dbl>          <dbl>       <dbl>             <dbl>
1 preferred (age_dec + ag…     1.00          0.426           0             0.458
mcmc_rhat(rh) + yaxis_text(hjust = 1) + theme_minimal(base_size = 10)

1.4 Interpreting the posterior coefficients

Coefficients are on the log-wage scale (preferred model). A negative female coefficient means lower expected log wages for women, holding the other covariates fixed; wbho (Non-White) is similarly signed in the original write-up; married and the age terms pick up the remaining associations.

posterior_interval(post, prob = 0.9) |>
  as.data.frame() |>
  tibble::rownames_to_column("parameter") |>
  as_tibble()
# A tibble: 9 × 3
  parameter        `5%`    `95%`
  <chr>           <dbl>    <dbl>
1 (Intercept)    1.67    1.83   
2 female        -0.220  -0.185  
3 wbho          -0.157  -0.119  
4 married        0.139   0.178  
5 age_dec        0.569   0.655  
6 I(age_dec^2)  -0.0679 -0.0575 
7 sigma          0.569   0.582  
8 log-fit_ratio -0.0198  0.00203
9 R2             0.145   0.164

Somewhat predictably different, perhaps, in the direction but not so in the magnitude, one might say — relative to the GLD prior medians used above.

# Same ridge style as the original note. Age terms are plotted separately so a
# shared axis does not flatten them into a spike at zero.
mcmc_areas_ridges(
  as.array(post),
  pars = c("female", "wbho", "married", "sigma")
)

mcmc_areas_ridges(
  as.array(post),
  pars = c("age_dec", "I(age_dec^2)")
)

# Original pairs view: coefficients + mean_PPD; exclude intercept / R2 / sigma labels soup.
pairs(post, regex_pars = "^[^(lRs]")

1.5 Model diagnostics and calibration

loo_quad <- loo(post_quad, save_psis = TRUE)
loo_pref <- loo(post, save_psis = TRUE)
print(loo_compare(loo_quad, loo_pref))
          elpd_diff se_diff
post         0.0       0.0 
post_quad -329.3      20.2 
tibble(
  model = c("age_dec^2 only (2020 baseline)", "age_dec + age_dec^2 (preferred)"),
  elpd_loo = c(
    loo_quad$estimates["elpd_loo", "Estimate"],
    loo_pref$estimates["elpd_loo", "Estimate"]
  ),
  p_loo = c(
    loo_quad$estimates["p_loo", "Estimate"],
    loo_pref$estimates["p_loo", "Estimate"]
  )
) |>
  arrange(desc(elpd_loo))
# A tibble: 2 × 3
  model                           elpd_loo p_loo
  <chr>                              <dbl> <dbl>
1 age_dec + age_dec^2 (preferred)  -10347.  7.40
2 age_dec^2 only (2020 baseline)   -10676.  6.39
pp_check(post, plotfun = "dens_overlay", nreps = 50)

pp_check(post, plotfun = "stat_2d", stat = c("mean", "sd"))

# LOO-PIT: floating-point / PSIS edge cases can yield PIT slightly above 1
# (bayesplot rounds those to 1). A wavy curve is the substantive signal.
pit_quad <- suppressWarnings(pp_check(post_quad, plotfun = "loo_pit_overlay"))
Some PIT values larger than 1! Largest:  1 
Rounding PIT > 1 to 1.
pit_pref <- suppressWarnings(pp_check(post, plotfun = "loo_pit_overlay"))
Some PIT values larger than 1! Largest:  1 
Rounding PIT > 1 to 1.
pit_quad

pit_pref

We see that the quadratic-only model — as in the original note — is often a bit lacking in the lower tail. Adding a linear age term improves ELPD substantially and is the preferred fit here; residual lower-tail tension on the PIT plot is the honest remaining limitation (measurement error and heteroskedasticity in CPS wages are the usual suspects).

Divergent transitions are mitigated with adapt_delta = 0.999. If any remain in the pairs plot, further simplifying the geometry would be the next modeling step.