Load the data

The data and some of the examples are derivatives of the sources noted in data/README.md (plus a couple more).

What we’ve got as the dependent variable is the percentage change in GOP support from 2012 to 2016.

d <- readr::read_csv("data/election_counties.csv", show_col_types = FALSE)
glimpse(d)
Rows: 3,111
Columns: 11
$ gop_support_change <dbl> 1.1, 0.0, 8.1, 5.3, 3.9, 3.1, 5.1, 5.8, 8.6, 9.2, 3…
$ pop_change         <dbl> 1.5, 9.8, -2.1, -1.8, 0.7, -1.4, -3.1, -2.3, -0.3, …
$ age_65_plus        <dbl> 13.8, 18.7, 16.5, 14.8, 17.0, 14.9, 18.0, 16.0, 18.…
$ black              <dbl> 18.7, 9.6, 47.6, 22.1, 1.8, 70.1, 44.0, 21.1, 39.5,…
$ hispanic           <dbl> 2.7, 4.6, 4.5, 2.1, 8.7, 7.5, 1.2, 3.5, 2.0, 1.5, 7…
$ hs_grad            <dbl> 85.6, 89.1, 73.7, 77.5, 77.0, 67.8, 76.3, 78.6, 75.…
$ undergrad          <dbl> 20.9, 27.7, 13.4, 12.1, 12.1, 12.5, 14.0, 16.1, 11.…
$ homeownership_rate <dbl> 76.8, 72.6, 67.7, 79.0, 81.0, 74.3, 70.3, 68.7, 67.…
$ median_home_value  <dbl> 136200, 168600, 89200, 90500, 117100, 70600, 74700,…
$ median_income      <dbl> 53682, 50221, 32911, 36447, 44145, 32033, 29918, 39…
$ poverty_rate       <dbl> 12.1, 13.9, 26.7, 18.1, 15.8, 21.6, 28.4, 21.9, 24.…

The independent variables (all continuous): pop_change, age_65_plus, black, hispanic, hs_grad, undergrad, homeownership_rate, median_home_value, median_income, poverty_rate.

Frequentist lm()

You could think about what model to use, but as an exemplar we are more interested in comparing notes than in the actual analysis. We’ll therefore proceed to lm().

lm_freq <- lm(gop_support_change ~ ., data = d)
summary(lm_freq)

Call:
lm(formula = gop_support_change ~ ., data = d)

Residuals:
    Min      1Q  Median      3Q     Max 
-40.268  -4.292  -0.703   3.870  52.981 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -1.088e+01  4.161e+00  -2.614  0.00898 ** 
pop_change         -2.693e-01  3.920e-02  -6.869 7.78e-12 ***
age_65_plus         2.066e-01  4.508e-02   4.583 4.76e-06 ***
black              -1.095e-01  1.161e-02  -9.434  < 2e-16 ***
hispanic           -1.526e-01  1.283e-02 -11.890  < 2e-16 ***
hs_grad             2.606e-01  3.697e-02   7.051 2.19e-12 ***
undergrad          -7.169e-01  3.000e-02 -23.895  < 2e-16 ***
homeownership_rate -1.572e-04  2.509e-02  -0.006  0.99500    
median_home_value  -1.748e-05  3.007e-06  -5.815 6.70e-09 ***
median_income       1.541e-04  3.166e-05   4.866 1.19e-06 ***
poverty_rate        2.275e-01  4.416e-02   5.151 2.75e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.482 on 3100 degrees of freedom
Multiple R-squared:  0.4929,    Adjusted R-squared:  0.4912 
F-statistic: 301.3 on 10 and 3100 DF,  p-value: < 2.2e-16

On this fit, a one-percentage-point higher Hispanic population share is associated with about \(-0.15\) points of GOP-support change (holding other covariates fixed), with a conventional \(p\)-value far below \(0.001\). The Bayesian sections ask what a full posterior says about the same slope.

Bayesian methods

rstanarm::stan_glm()

Useful references: stan_glm, continuous outcomes, bayesplot MCMC, Stan warnings.

rstanarm emulates familiar R model formulas but samples with Stan. Current default priors are autoscaled Student-\(t\) / normal families (not the older fixed normal(0, 10) / normal(0, 5) values from early docs). Leaving priors unspecified is still a valid weakly informative choice.

post1 <- stan_glm(
  gop_support_change ~ .,
  data = d,
  family = gaussian(),
  seed = 42,
  refresh = 0
)
print(post1, digits = 3)
stan_glm
 family:       gaussian [identity]
 formula:      gop_support_change ~ .
 observations: 3111
 predictors:   11
------
                   Median  MAD_SD 
(Intercept)        -10.965   4.088
pop_change          -0.269   0.039
age_65_plus          0.206   0.044
black               -0.109   0.012
hispanic            -0.152   0.012
hs_grad              0.261   0.037
undergrad           -0.716   0.030
homeownership_rate   0.001   0.025
median_home_value    0.000   0.000
median_income        0.000   0.000
poverty_rate         0.227   0.043

Auxiliary parameter(s):
      Median MAD_SD
sigma 7.480  0.097 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg
loo1 <- loo(post1)
loo1

Computed from 4000 by 3111 log-likelihood matrix.

         Estimate    SE
elpd_loo -10685.3  72.7
p_loo        19.0   2.5
looic     21370.6 145.4
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume independent draws (r_eff=1).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.
diag1 <- tibble(
  max_rhat = max(rhat(post1), na.rm = TRUE),
  min_neff_ratio = min(neff_ratio(post1), na.rm = TRUE)
)
diag1
# A tibble: 1 × 2
  max_rhat min_neff_ratio
     <dbl>          <dbl>
1     1.00          0.618
draws <- as.data.frame(as.matrix(post1))
colnames(draws)[1] <- "intercept"

ggplot(d, aes(x = hispanic, y = gop_support_change)) +
  geom_point(size = 0.3, color = "#003366", alpha = 0.35) +
  geom_abline(
    data = draws, aes(intercept = intercept, slope = hispanic),
    color = "skyblue", linewidth = 0.2, alpha = 0.05
  ) +
  geom_abline(
    intercept = coef(post1)["(Intercept)"],
    slope = coef(post1)["hispanic"],
    color = "skyblue4", linewidth = 1
  ) +
  theme_minimal()

pp_check(post1, plotfun = "hist", nreps = 5)

mcmc_intervals() defaults to 50% (thick) and 90% (thin) central intervals with the point at the posterior median.

slope_pars <- setdiff(names(coef(post1)), c("(Intercept)"))
mcmc_intervals(
  as.matrix(post1),
  pars = slope_pars,
  prob = 0.5,
  prob_outer = 0.9
)

m1 <- as.matrix(post1)
median(m1[, "hispanic"])
[1] -0.1520994
round(quantile(m1[, "hispanic"], probs = c(0.1, 0.5, 0.9)), 3)
   10%    50%    90% 
-0.168 -0.152 -0.136 
mcmc_intervals(
  as.matrix(post1),
  pars = slope_pars,
  prob = 0.8,
  prob_outer = 0.99
)

Conditional on the data and this model, the posterior for the Hispanic slope is tightly concentrated near \(-0.15\): a 99% central interval is roughly \(-0.17\) to \(-0.13\).

mcmc_areas(as.matrix(post1), pars = c("age_65_plus", "undergrad", "poverty_rate"))

Age \(65+\) and poverty rate pull toward positive GOP-support change in this period; bachelor’s share pulls negative, with a narrower central mass.

pp_check(post1)

The density overlay shows mild mismatch in the bulk (medians slightly off), which could be remedied e.g. by changing our priors.

Custom Stan (rstan)

Writing the likelihood and priors in Stan makes the model explicit. Below is one centered/scaled design with Student-\(t\) priors (stan/mlr.stan). Flat-prior versions of this regression mix poorly on modern Stan; weakly informative priors are the pedagogical default here.

X_scaled <- scale(as.matrix(d |> select(-gop_support_change)))
stan_dat <- list(
  n = nrow(d),
  k = ncol(X_scaled),
  X = X_scaled,
  Y = d$gop_support_change
)
post2 <- stan(
  file = "stan/mlr.stan",
  data = stan_dat,
  chains = 4,
  iter = 2000,
  warmup = 1000,
  seed = 42,
  refresh = 0,
  control = list(adapt_delta = 0.95)
)
sum2 <- summary(post2)$summary
diag2 <- tibble(
  max_rhat = max(sum2[, "Rhat"], na.rm = TRUE),
  min_n_eff = min(sum2[, "n_eff"], na.rm = TRUE)
)
diag2
# A tibble: 1 × 2
  max_rhat min_n_eff
     <dbl>     <dbl>
1     1.00     1756.
beta_pars <- paste0("beta[", seq_len(stan_dat$k), "]")
color_scheme_set("yellow")
mcmc_intervals(post2, pars = beta_pars) +
  scale_y_discrete(labels = colnames(X_scaled))

color_scheme_set("yellow")
mcmc_areas(post2, pars = c("beta[4]", "beta[7]")) +
  labs(title = "Posterior Distributions") +
  scale_y_discrete(labels = c("hispanic", "homeownership_rate"))

color_scheme_set("blue")

Slopes here are in scaled-predictor units, so magnitudes are not directly comparable to the stan_glm coefficients on the raw scale. The qualitative signs for Hispanic and homeownership still line up.

Prior elicitation with GLD / J-QPD

Johnson quantile-parameterized distributions (J-QPD) and the generalized lambda distribution (GLD) let you state priors through quantiles rather than a named family. Helpers live in stan/quantile_functions.stan and R/GLD_helpers.R.

Draw \(S\) prior predictive \(R^2\) values:

\[ R_s^2=\frac{\frac{1}{N-1}\sum_{n=1}^{N}(\tilde\mu_n-\bar\mu)^2}{\tilde\sigma^2+\frac{1}{N-1}\sum_{n=1}^{N}(\tilde\mu_n-\bar\mu)^2}. \]

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

X_c <- model.matrix(gop_support_change ~ ., data = d)[, -1, drop = FALSE]
X_c <- sweep(X_c, 2, colMeans(X_c), `-`)
m_alpha <- -10.9
s_alpha <- 4.1

lowers <- c(
  PC = -0.0005, A = 0.0001, B = -0.0004, H = -0.0004, HS = 0.0023,
  U = -0.0009, HR = -0.0003, MHV = -0.0003, MI = -0.0003, PR = -0.0004
)
medians <- c(
  PC = -0.0004, A = 0.00015, B = -0.0003, H = -0.0003, HS = 0.0025,
  U = -0.0008, HR = -0.0002, MHV = -0.0002, MI = -0.0002, PR = -0.0003
)
uppers <- c(
  PC = -0.0003, A = 0.0002, B = -0.0002, H = -0.0002, HS = 0.0027,
  U = -0.0007, HR = -0.0001, MHV = -0.0001, MI = -0.0001, PR = -0.0002
)
ninety <- c(
  PC = -0.0002, A = 0.00025, B = -0.0001, H = -0.0001, HS = 0.0029,
  U = -0.0006, HR = 0.00, MHV = 0.00, MI = 0.00, PR = -0.0001
)

a_s <- lapply(seq_along(medians), function(i) {
  GLD_solver(lowers[i], medians[i], uppers[i], ninety[i], alpha = 0.9)
})
names(a_s) <- names(medians)

q_sigma <- c(lower = 0.25, median = 0.5, upper = 1)

These quantile targets are deliberately tight and inspired by the stan_glm posterior from the preceding chunks.

feat <- c("pop_change", "age_65_plus", "black", "hispanic", "hs_grad",
          "undergrad", "homeownership_rate", "median_home_value",
          "median_income", "poverty_rate")
keys <- c("PC", "A", "B", "H", "HS", "U", "HR", "MHV", "MI", "PR")

R_squared <- replicate(1000, {
  alpha_ <- rnorm(1, mean = m_alpha, sd = s_alpha)
  betas <- vapply(seq_along(keys), function(j) {
    k <- keys[j]
    GLD_rng(
      medians[k],
      IQR = uppers[k] - lowers[k],
      asymmetry = a_s[[k]][1],
      steepness = a_s[[k]][2]
    )
  }, numeric(1))
  mu_ <- as.numeric(alpha_ + X_c[, feat, drop = FALSE] %*% betas)
  sigma_ <- JQPDS_rng(lower_bound = 0, alpha = 0.25, quantiles = q_sigma)
  numer <- var(mu_)
  numer / (sigma_^2 + numer)
})
summary(R_squared)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.1082  0.9938  0.9989  0.9796  0.9998  1.0000

Posterior with an \(R^2\) prior (stan_lm)

post4 <- stan_lm(
  gop_support_change ~ .,
  data = d,
  prior = R2(location = median(R_squared), what = "median"),
  prior_intercept = normal(location = m_alpha, scale = s_alpha, autoscale = FALSE),
  seed = 42,
  refresh = 0
)
print(post4, digits = 2)
stan_lm
 family:       gaussian [identity]
 formula:      gop_support_change ~ .
 observations: 3111
 predictors:   11
------
                   Median MAD_SD
(Intercept)        -10.95   2.89
pop_change          -0.27   0.04
age_65_plus          0.20   0.04
black               -0.11   0.01
hispanic            -0.15   0.01
hs_grad              0.26   0.03
undergrad           -0.72   0.03
homeownership_rate   0.00   0.02
median_home_value    0.00   0.00
median_income        0.00   0.00
poverty_rate         0.23   0.04

Auxiliary parameter(s):
              Median MAD_SD
R2            0.49   0.01  
log-fit_ratio 0.00   0.01  
sigma         7.47   0.10  

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg
diag4 <- tibble(
  max_rhat = max(rhat(post4), na.rm = TRUE),
  min_neff_ratio = min(neff_ratio(post4), na.rm = TRUE)
)
diag4
# A tibble: 1 × 2
  max_rhat min_neff_ratio
     <dbl>          <dbl>
1     1.00          0.591
loo4 <- loo(post4)
loo_compare(loo1, loo4)
      elpd_diff se_diff
post4  0.0       0.0   
post1 -0.3       0.2   
ppc_intervals(
  y = d$gop_support_change,
  yrep = posterior_predict(post4),
  x = d$hispanic
) +
  xlab("hispanic")

m4 <- as.matrix(post4)
slope_names <- setdiff(colnames(m4), c("(Intercept)", "sigma", "R2", "log-fit_ratio"))
slope_names <- slope_names[slope_names %in% colnames(d)]

posterior_quantiles <- purrr::map_dfr(slope_names, function(nm) {
  q <- quantile(m4[, nm], probs = c(0.025, 0.25, 0.5, 0.75, 0.975))
  tibble(
    term = nm,
    q025 = unname(q[1]),
    q25 = unname(q[2]),
    q50 = unname(q[3]),
    q75 = unname(q[4]),
    q975 = unname(q[5])
  )
})
posterior_quantiles
# A tibble: 10 × 6
   term                     q025        q25        q50        q75       q975
   <chr>                   <dbl>      <dbl>      <dbl>      <dbl>      <dbl>
 1 pop_change         -0.344     -0.294     -0.268     -0.241     -0.190    
 2 age_65_plus         0.120      0.175      0.205      0.234      0.293    
 3 black              -0.132     -0.117     -0.109     -0.102     -0.0865   
 4 hispanic           -0.176     -0.160     -0.152     -0.144     -0.129    
 5 hs_grad             0.200      0.240      0.261      0.282      0.322    
 6 undergrad          -0.774     -0.736     -0.716     -0.696     -0.660    
 7 homeownership_rate -0.0462    -0.0152     0.000473   0.0159     0.0472   
 8 median_home_value  -0.0000234 -0.0000194 -0.0000174 -0.0000154 -0.0000116
 9 median_income       0.0000947  0.000133   0.000153   0.000175   0.000213 
10 poverty_rate        0.153      0.202      0.228      0.254      0.302

Conclusion

Three fits on the same county panel:

  1. OLS — point estimates and \(p\)-values.
  2. stan_glm — default weakly informative priors, full posterior, LOO, PPC.
  3. Custom Stan + stan_lm($R^2$) — explicit Student-\(t\) priors on a scaled design, then a quantile-elicited \(R^2\) prior that folds GLD / J-QPD into rstanarm.

It may not be quite trivial to “guess” the prior stuff using quantile functions. But the tradeoff is that we can dispense with pre-defined distributions — a cleverer approach, but more difficult to implement.