Week 7: Multilevel Models

Previously

• Markov Chain Monte Carlo
  • Construct a Markov chain whose stationary distribution is the posterior
  • Metropolis-Hastings, Gibbs sampling
• Stan
  • Hamiltonian Monte Carlo - use the gradient of the log posterior to get good M-H proposals
    • Avoids proposals that are either too auto-correlated or too likely to be rejected
  • Common structure to statistical modeling: data, parameters, model
  • Diagnosing convergence: trace plots, \(\hat{R}\), effective sample size

This week

• Multilevel regression
  • Extending the Bayesian GLM by adding additional grouping and structure
  • “Random effects” models/“Partial pooling”
• Diagnosing model fit
  • Posterior predictive checks
  • Information criteria vs. Cross-validation
• Post-stratification
  • Adjusting non-representative surveys to match population targets
  • Combining with multilevel regression: “MrP”

Hierarchical Models

• In our regression models so far, we assumed a common intercept and slope for all of the data.

  • The underlying model is the same regardless of which group (e.g. state) an observation is located in.

• We also assumed that observations were completely exchangeable

  • \(Y_i \sim \text{Normal}(X_i^{\prime}\beta, \sigma^2)\) implies no residual correlation in the outcomes across all observations.

• This may be a poor modeling choice if states are very heterogeneous and information about the state in which an observation is located would provide additional predictive power beyond knowing \(X_i\).

  • Additionally, exchangeability of observations may be violated - unobserved factors may lead observations within a particular state to be correlated (e.g. if all respondents in a state are exposed to some common shock).

Hierarchical Models

• There are lots of choices about how we want to incorporate hierarchy into the model the coefficients across groupings
  • What levels of aggregation and how many?
  • Which parameters should be allowed to vary (all of the coefficients or just the intercepts)?
  • Where does the structure enter into the model (the means? the variances? the covariances?)

Application: Education polarization

• A growing divide in American politics centers on education
  • Voters with a college degree have moved towards the Democrats over the last decade
  • Voters without one have moved towards the Republicans
• Is the education gap the same everywhere?
  • Gelman et al. (2007) found that the relationship between income and vote choice differs across states: steep in poor states like Mississippi, nearly flat in rich states like Connecticut
  • Does the relationship between education and vote choice differ across states too?

Data: CES 2016

• The 2016 Cooperative Congressional Election Study (Ansolabehere and Schaffner 2017; via Kuriwaki 2026)
  • A large (64,600 respondents) opt-in online survey run by YouGov
  • Respondents matched to state voter files.
• Our sample: respondents with a validated vote record who say they voted for Clinton or Trump

ces <- read_csv("data/ces2016.csv") %>%
  mutate(college = as.integer(educ %in% c("BA", "Postgrad"))) # Has a 4-year degree
voters <- ces %>% filter(!is.na(clinton))

• 28,462 voters; Clinton won 65% of college graduates and 47% of everyone else

Data: CES 2016

• The “raw” education gap in each state: Clinton’s share among college graduates minus her share among non-graduates

raw_gaps <- voters %>% group_by(state) %>%
  summarize(n = n(), gap = mean(clinton[college == 1]) - mean(clinton[college == 0]))

• Is this variation real? Or is it noise from the small samples in small states?

Complete pooling

• Let’s start with a single logistic regression for everyone.

  • \(Y_{ij} = 1\) if respondent \(i\) in state \(j\) voted for Clinton, \(0\) if Trump.

  \[Y_{ij} \sim \text{Bernoulli}(\pi_{ij})\] \[\log\bigg(\frac{\pi_{ij}}{1 - \pi_{ij}}\bigg) = \beta_0 + \beta_1 \text{college}_{ij}\] \[\beta_0, \beta_1 \sim \text{Normal}(0, 5^2)\]

• Our only predictors are state and education, so every respondent in the same state \(\times\) education cell has the same \(\pi_{ij}\)

  • We’ll fit the model to all 28,462 respondents
  • But anything computed after estimation (predictions, model checks) only needs to be computed once for each of the 102 cells

Complete pooling

• Set up the respondent-level data, and a table of the cells

voters <- voters %>% mutate(state_num = as.numeric(as.factor(state)))
X_mat <- model.matrix(~ college, data = voters)
K <- ncol(X_mat) # Number of beta parameters
cells <- voters %>% group_by(state, state_num, college) %>%
  summarize(n = n(), y = sum(clinton), .groups = "drop")
X_cell <- model.matrix(~ college, data = cells)

Complete pooling

pooled_model <- "
data {
  int N; // number of respondents
  int K; // number of covariates
  matrix[N, K] X; // covariate matrix
  array[N] int<lower=0, upper=1> y; // 1 = Clinton, 0 = Trump
  int M; // number of cells
  matrix[M, K] X_cell; // covariates for each cell
  array[M] int n_cell; // number of respondents in each cell
}
parameters {
  vector[K] beta; // regression coefficients
}
model {
  beta ~ normal(0, 5); // normal prior
  y ~ bernoulli_logit(X * beta); // logistic regression
}
generated quantities {
  vector[M] eta_cell = X_cell * beta; // linear predictor for each cell
  array[M] int y_rep = binomial_rng(n_cell, inv_logit(eta_cell)); // replicated Clinton voters in each cell
  vector[M] ll_clinton = log_inv_logit(eta_cell); // log-likelihood of one Clinton voter in each cell
  vector[M] ll_trump = log1m_inv_logit(eta_cell); // log-likelihood of one Trump voter in each cell
}
"

Complete pooling

• Pass the data as a list and run Stan!

stan_data <- list(N = nrow(voters), K = K, X = X_mat, y = voters$clinton,
                  J = max(voters$state_num), state = voters$state_num,
                  M = nrow(cells), X_cell = X_cell, n_cell = cells$n, cell_state = cells$state_num)

pooled_fit <- stan(
  model_code = pooled_model,  # Stan code
  data = stan_data,       # named list of data (Stan ignores the elements it doesn't use)
  chains = 4,             # number of Markov chains
  warmup = 500,           # number of warmup iterations per chain
  iter = 2500,            # total number of iterations per chain
  cores = 4,              # number of cores
  refresh = 0,            # no progress shown
  seed = 53706
  )

print(pooled_fit, pars = "beta")

Inference for Stan model: anon_model.
4 chains, each with iter=2500; warmup=500; thin=1; 
post-warmup draws per chain=2000, total post-warmup draws=8000.

         mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
beta[1] -0.10       0 0.02 -0.13 -0.11 -0.10 -0.09 -0.07  3090    1
beta[2]  0.71       0 0.03  0.67  0.70  0.71  0.73  0.76  3112    1

Samples were drawn using NUTS(diag_e) at Wed Oct 21 11:37:11 2026.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Diagnosing model fit

Posterior predictive distribution

• A very useful quantity to obtain is the predictive distribution of \(y\) given the model.

  \[f(\tilde{y} | \mathbf{Y}) =\int f(\tilde{y}, \theta | \mathbf{Y}) d\theta = \int f(\tilde{y} | \mathbf{Y}, \theta) f(\theta | \mathbf{Y}) d\theta\]

• Suppose our prediction target is the in-sample distribution of the covariates

  • We can get the posterior predictive distribution directly from our MCMC algorithm by drawing a value of \(y\) for each value of \(\theta\) that we sample.

• (e.g.) for the logit

generated quantities {
  array[N] int y_rep = bernoulli_logit_rng(X * beta); //generate coin flips for each sample
}

• It is common to use the posterior predictive distribution as a model diagnostic (“posterior predictive checking”)
  1. Check in-sample fit: Does the model actually reproduce the distribution of the data?
  2. Check against simulated data: Does the model recover the simulated distributions?

Posterior predictive check

• In our model, y_rep is a replicated count of Clinton voters for every cell
  • Drawing \(n_m\) Bernoulli outcomes with the same \(\pi_m\) and adding them up is the same as one Binomial(\(n_m\), \(\pi_m\)) draw
• Our model should be able to reproduce Clinton’s share of the vote in each state

state_n <- tapply(cells$n, cells$state, sum) # Respondents per state
state_share_ppc <- function(fit){
  y_rep <- rstan::extract(fit, "y_rep")$y_rep # draws x cells
  share_rep <- sweep(t(rowsum(t(y_rep), cells$state)), 2, state_n, "/") # draws x states
  tibble(state = names(state_n), obs = tapply(cells$y, cells$state, sum)/state_n, mean = colMeans(share_rep),
         lower = apply(share_rep, 2, quantile, .025), upper = apply(share_rep, 2, quantile, .975))
}
ppc_pooled <- state_share_ppc(pooled_fit)

Posterior predictive check

• 24 of 51 states fall outside their 95% predictive intervals
  • The model thinks every state looks like the nation (up to its share of college graduates)

Varying intercepts model

• Our regression from before is an example of a fully pooled regression.

  • Within individual states though, the regression line may be a poor predictor
  • The simplest fix is to allow the “intercept” to shift across units.

• Assume:

  \[\log\bigg(\frac{\pi_{ij}}{1 - \pi_{ij}}\bigg) = X_{ij}^{\prime}\beta + \alpha_{j}\] \[\alpha_j \sim \text{Normal}(0, \sigma^2_\alpha)\]

  And keep the same prior as before on \(\beta\), which still contains the intercept. For \(\sigma_\alpha\) we’ll use a Half-Cauchy prior (a Cauchy with the constraint \(\sigma_\alpha > 0\)).

• \(\alpha_j\) can be interpreted as the group-specific “shift” in the intercept from the “grand” intercept \(\beta_0\)

  • An equivalent parameterization would be to remove the intercept from \(\beta\) and write \(\alpha_j \sim \text{Normal}(\mu_\alpha, \sigma^2_\alpha)\)

Varying intercepts model

• Let’s implement this in Stan - first, our data block

data{
  int N; // number of respondents
  int K; // number of covariates
  matrix[N, K] X; //matrix of covariates
  array[N] int<lower=0, upper=1> y; // 1 = Clinton, 0 = Trump
  int J; // Number of groups
  array[N] int state; // group membership indicator
  int M; // number of cells
  matrix[M, K] X_cell; // covariates for each cell
  array[M] int n_cell; // number of respondents in each cell
  array[M] int cell_state; // group of each cell
}

Varying intercepts model

• Next, our parameters.

parameters{
  vector[K] beta; // beta coefficients
  vector[J] alpha; // random intercepts
  real<lower=0> sigma_a; // standard deviation of intercepts
}

• Finally, our model

model{
  beta ~ normal(0, 5); // normal prior on coefficients
  alpha ~ normal(0, sigma_a); // normal distribution on random intercepts
  sigma_a ~ cauchy(0, 2.5); // half-Cauchy, since sigma_a is constrained to be positive
  y ~ bernoulli_logit(alpha[state] + X*beta);
}

Varying intercepts model

• Put it into Stan! The generated quantities block computes the same cell-level quantities as before

generated quantities {
  vector[M] eta_cell = alpha[cell_state] + X_cell*beta;
  array[M] int y_rep = binomial_rng(n_cell, inv_logit(eta_cell));
  vector[M] ll_clinton = log_inv_logit(eta_cell);
  vector[M] ll_trump = log1m_inv_logit(eta_cell);
}

model_var_intercept <- stan(
  model_code = varying_intercepts_model,  # Stan code
  data = stan_data,       # named list of data
  chains = 4,             # number of Markov chains
  warmup = 500,           # number of warmup iterations per chain
  iter = 2500,            # total number of iterations per chain
  cores = 4,              # number of cores (this is much slower if = 1)
  refresh = 0,            # no progress shown
  seed = 53706            # random seed
  )

Varying intercepts model

• Summarize the fit

print(model_var_intercept, pars = c("beta", "sigma_a", "alpha[1]"))

Inference for Stan model: anon_model.
4 chains, each with iter=2500; warmup=500; thin=1; 
post-warmup draws per chain=2000, total post-warmup draws=8000.

          mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
beta[1]  -0.13       0 0.05 -0.22 -0.16 -0.13 -0.10 -0.04   658    1
beta[2]   0.69       0 0.03  0.64  0.68  0.69  0.71  0.75 12103    1
sigma_a   0.31       0 0.04  0.24  0.28  0.31  0.34  0.40  6578    1
alpha[1] -0.06       0 0.21 -0.47 -0.20 -0.06  0.08  0.34 10426    1

Samples were drawn using NUTS(diag_e) at Wed Oct 21 11:39:31 2026.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Varying intercepts model

• Plot the deviations from the “grand mean” intercept by state

Posterior predictive check

• Does the varying intercepts model reproduce the state shares?

ppc_varint <- state_share_ppc(model_var_intercept)

• Now 1 of 51 states fall outside their 95% predictive intervals

Model Diagnostics

• There are two general approaches to evaluating model quality and comparing the predictive power of two or more statistical models given the same data
  • Information criteria
  • Cross-validation
• Both approaches actually share common goals
  • Reward quality of predictions
  • Penalize model complexity
• We want to avoid models that overfit to the observed data
  • Information criteria address this with an explicit penalty term for the number of parameters
  • Cross-validation methods do this by evaluating fit on a “held-out” part of the data that’s not used to estimate the model.

Information Criteria

• Frequentist/MLE Inference - the Akaike Information Criterion

  \[\text{AIC} = 2k - 2\ell(\hat{\theta})\]

  • \(k\) is the number of parameters and \(\ell(\hat{\theta})\) is the value of the maximized log-likelihood.
  • Pick the model with the lowest AIC

• Bayesian Inference - the Deviance Information Criterion

  \[\text{DIC} = 2p_{DIC} - 2\log p(y | \hat{\theta}_{\text{PM}})\]

  • \(\hat{\theta}_{\text{PM}}\) is the posterior mean and
  • the “effective number of parameters” is \(p_{DIC} = 2(\log p(y | \hat{\theta}_{\text{PM}}) - E_{\text{post}}[\log p(y|\theta)])\)

Log-predictive density

• In general, a very common objective function for evaluating the quality of model “fit” in a Bayesian setting is the log-predictive density

  • Think of it as a more general version of criteria like “mean-squared error”

• For a new “out-of-sample” \(\tilde{y}_i\), the log posterior predictive density is

  \[\log p_{\text{post}} (\tilde{y}_i) = \log \int p(\tilde{y}_i | \theta) \times p_{\text{post}}(\theta)d\theta\]

• We want to pick a model that has a high expected log pointwise-predictive density – for a “new” size \(n\) sample of observation drawn from a common distribution \(f(\tilde{y}_i)\),

  \[\text{elppd} = \sum_{i=1}^n E_f[\log p_{\text{post}} (\tilde{y}_i)] = \sum_{i=1}^n\int \log p_{\text{post}} (\tilde{y}_i)f(\tilde{y}_i)d\tilde{y}_i\]

• Our goal is to try to estimate this from the data that we have.

Estimating elppd

• One approach would be to use the “in-sample” fit criterion using in-sample observations \(y_i\)

  \[\text{lppd} = \sum_{i=1}^n \log p_{\text{post}} (y_i)\]

• We can calculate it by averaging over our \(S\) posterior samples

  \[\widehat{\text{lppd}} = \sum_{i=1}^n \log \bigg( \frac{1}{S} \sum_{s=1}^S p(y_i | \theta^s)\bigg)\]

• But here we know that it’s going to be biased (too optimistic) for the “out-of-sample” predictive accuracy due to over-fitting.

Estimating elppd

• An improvement would be to correct the “in-sample” fit using a measure of the model complexity
  • This is the intuition behind AIC, DIC, (and other information criteria like WAIC)
• An alternative would be to use cross-validation:
  • Split the sample into \(K\) “folds” (typical choice is 10)
  • For each fold \(k\), estimate the model on the \(K-1\) other folds
  • Using that model, calculate the log-predictive density for the observations in the \(k\)th fold
  • The \(k\)th fold is the “held-out” (or “test” set), the \(K-1\) remaining folds are the “training” set
• When \(K = N\), this is often known as “leave-one out” cross-validation

Leave-one-out CV

• From Vehtari, Gelman, and Gabry (2017), our “leave-one-out” estimate of the expected log predictive density is:

  \[\text{elpd}_{\text{loo}} = \sum_{i=1}^n \log p(y_i | y_{-i})\]

• In principle, this requires re-fitting the model to \(N\) datasets (with a single observation removed). However, it is possible to obtain an estimate of \(\log p(y_i | y_{-i})\) using the full model fit. Note:

  \[p(y_i | y_{-i}) = \int p(y_i | \theta) p(\theta | y_{-i}) d\theta\]

• We can calculate \(p(y_i | \theta)\) for a given draw of theta.

  • However, we only have draws from \(p(\theta | y)\) (our posterior).
  • The solution is to re-weight our draws such that, averaging over them yields an integral over \(p(\theta | y_{-i})\)

Leave-one-out CV

• In order to use loo, we need the log-likelihood of each observation \(\log p(y_i | \theta^s)\), evaluated at each posterior draw \(s\). This comes by default in canned rstanarm or brms models but for custom models we need to specify to return it in the generated quantities block
• For our logistic regression, each respondent’s log-likelihood:

generated quantities {
  vector[N] log_lik;
  for (i in 1:N) log_lik[i] = bernoulli_logit_lpmf(y[i] | X[i] * beta);
}

• loo() then computes PSIS-LOO from this \(N \times S\) matrix, and loo_compare() compares models by summing the pointwise differences in \(\log p(y_i | y_{-i})\)

Leave-one-out CV

• For our example: \(N \times S\) is 28,462 respondents \(\times\) 8,000 draws \(\approx\) 228 million numbers (1.8 GB) for each model

• So we’ll instead do the same calculation by state x education cell.

  • Every Clinton voter in cell \(m\) has the same log-likelihood, \(\log \pi_m\) (ll_clinton); every Trump voter has \(\log(1 - \pi_m)\) (ll_trump)
  • So they also have the same leave-one-out density
  • Run loo once for each cell \(\times\) vote, then replicate each result once per respondent who falls in that cell.

Leave-one-out CV

loo_respondents <- function(fit){
  # draws x chains x (cell, vote) combinations
  ll <- abind::abind(extract_log_lik(fit, "ll_clinton", merge_chains = F),
                     extract_log_lik(fit, "ll_trump", merge_chains = F), along = 3)
  loo_types <- loo(ll, r_eff = relative_eff(exp(ll)))
  rep(loo_types$pointwise[, "elpd_loo"], c(cells$y, cells$n - cells$y)) # one entry per respondent
}
elpd_summary <- function(pointwise) c(elpd_loo = sum(pointwise), se = sqrt(length(pointwise)*var(pointwise)))

loo_pooled <- loo_respondents(pooled_fit)
elpd_summary(loo_pooled)

elpd_loo       se 
-19215.0     31.6 

loo_varint <- loo_respondents(model_var_intercept)
elpd_summary(loo_varint)

elpd_loo       se 
-19041.6     36.2 

Leave-one-out CV

• Compare the two models:

elpd_summary(loo_varint - loo_pooled)

elpd_loo       se 
   173.3     18.4 

• The varying intercepts model predicts held-out respondents much better

Posterior predictive check

• The varying intercepts model still assumes the same education gap (on the logit scale) in every state. Does it reproduce the gap in each state?

state_gap_ppc <- function(fit){
  y_rep <- rstan::extract(fit, "y_rep")$y_rep
  col <- cells$college == 1
  gap_rep <- sweep(y_rep[, col], 2, cells$n[col], "/") - sweep(y_rep[, !col], 2, cells$n[!col], "/")
  tibble(state = cells$state[col], obs = cells$y[col]/cells$n[col] - cells$y[!col]/cells$n[!col], mean = colMeans(gap_rep),
         lower = apply(gap_rep, 2, quantile, .025), upper = apply(gap_rep, 2, quantile, .975))
}
gap_varint <- state_gap_ppc(model_var_intercept)

• 6 states fall outside: FL, KY, MI, MN, OR, WV

Varying slopes model

• What if we allow all of the parameters to vary by group?

  \[\beta_j \sim \text{Normal}(\mu_{\beta}, \Sigma_{\beta})\] \[\log\bigg(\frac{\pi_{ij}}{1 - \pi_{ij}}\bigg) = X_{ij}^{\prime}\beta_{j}\]

• This is a more general and flexible hierarchical model

  • But introduces a slight complication – we now need a prior distribution on an entire matrix \(\Sigma_{\beta}\) (the variances and covariances of the parameters).

• The conjugate prior is the inverse-Wishart, the multivariate extension of the inverse-gamma.

  • An alternative non-conjugate prior (that appears to have better performance) is a scaled correlation matrix with an LKJ prior on the correlation matrix (see Stan documentation for more).
  • The latter is now preferred (especially by the Stan programmers)

Varying slopes model

• Our prior on \(\Sigma_{\beta}\) starts by decomposing \(\Sigma_{\beta}\) into a diagonal scaling matrix \(\text{diag}(\tau)\) and a correlation matrix \(\Omega\)

  \[\Sigma_{\beta} = \text{diag}(\tau)\Omega\text{diag}(\tau)\]

• For each \(\tau_k\) we’ll use a Half-Cauchy distribution (Cauchy with the constraint \(\tau_k > 0\))

• For \(\Omega\), we’ll use the Lewandowski-Kurowicka-Joe (LKJ) distribution which defines a distribution over correlation matrices

  • Intuition: Similar to a beta distribution in multiple dimensions.

Varying slopes model

• Write each state’s coefficients as the national average plus a state-specific deviation

  \[\beta_j = \mu_{\beta} + e_j, \qquad e_j \sim \text{Normal}(0, \Sigma_{\beta})\] \[\log\bigg(\frac{\pi_{ij}}{1 - \pi_{ij}}\bigg) = \underbrace{X_{ij}^{\prime}\mu_{\beta}}_{\text{same in every state}} + \underbrace{X_{ij}^{\prime}e_j}_{\text{state deviation}}\]

• Same setup as above, but keeping the part that’s the same in every state separate will make it easy to add state-level predictors later

Varying slopes model

• Add a second design matrix X_fix for the terms that are the same in every state (here, \(X_{ij}^{\prime}\mu_{\beta}\))
  • For now it’s just X again - but it will get more columns later

stan_data$P <- K
stan_data$X_fix <- X_mat
stan_data$X_fix_cell <- X_cell

• Add them to the data block

  int P; // number of population-level coefficients
  matrix[N, P] X_fix; // population-level design matrix
  matrix[M, P] X_fix_cell; // population-level design matrix for each cell

Varying slopes model

• The parameters block

parameters{
  vector[P] b; // population-level coefficients (here, mu_beta)
  corr_matrix[K] Omega; //prior correlation
  vector<lower=0>[K] tau; //prior scale
  array[J] vector[K] e; // state deviations
}

• The model block

model{
  vector[N] eta; // linear predictor for each respondent
  tau ~ cauchy(0, 2.5);
  Omega ~ lkj_corr(2);
  b ~ normal(0, 5);
  e ~ multi_normal(rep_vector(0, K), quad_form_diag(Omega, tau));
  for (i in 1:N) {
    eta[i] = X_fix[i] * b + X[i] * e[state[i]];
  }
  y ~ bernoulli_logit(eta);
}

Varying slopes model

• The generated quantities compute the same cell-level quantities as before

generated quantities {
  vector[M] eta_cell;
  for (m in 1:M) {
    eta_cell[m] = X_fix_cell[m] * b + X_cell[m] * e[cell_state[m]];
  }
  array[M] int y_rep = binomial_rng(n_cell, inv_logit(eta_cell));
  vector[M] ll_clinton = log_inv_logit(eta_cell);
  vector[M] ll_trump = log1m_inv_logit(eta_cell);
}

model_var_slope <- stan(
  model_code = varying_slopes_model,  # Stan code
  data = stan_data,       # named list of data
  chains = 4,             # number of Markov chains
  warmup = 500,           # number of warmup iterations per chain
  iter = 2500,            # total number of iterations per chain
  cores = 4,
  refresh = 0,            # no progress shown
  seed = 53706            # random seed
  )

Varying slopes model

• Summarize the fit

print(model_var_slope, pars=c("b", "tau", "Omega[1,2]"))

Inference for Stan model: anon_model.
4 chains, each with iter=2500; warmup=500; thin=1; 
post-warmup draws per chain=2000, total post-warmup draws=8000.

            mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
b[1]       -0.13    0.00 0.04 -0.22 -0.16 -0.13 -0.10 -0.05  1283    1
b[2]        0.69    0.00 0.04  0.62  0.67  0.69  0.72  0.76  4428    1
tau[1]      0.28    0.00 0.04  0.21  0.25  0.28  0.31  0.38  3088    1
tau[2]      0.15    0.00 0.05  0.06  0.12  0.15  0.18  0.24   342    1
Omega[1,2]  0.31    0.01 0.27 -0.24  0.13  0.33  0.51  0.80  1025    1

Samples were drawn using NUTS(diag_e) at Wed Oct 21 20:25:51 2026.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

• b is \(\mu_{\beta}\), the national average intercept and education slope
• tau[2], the standard deviation of the education slopes across states, is clearly above zero

Varying slopes model

• Each state’s coefficients are \(\mu_{\beta} + e_j\): b plus that state’s e

b_vs <- rstan::extract(model_var_slope, "b")$b # draws x K
e_vs <- rstan::extract(model_var_slope, "e")$e # draws x states x K
int_draws <- b_vs[, 1] + e_vs[, , 1] # state intercepts
slope_draws <- b_vs[, 2] + e_vs[, , 2] # state education slopes

Posterior predictive check

• Does the varying slopes model reproduce the education gap in each state?

gap_varslope <- state_gap_ppc(model_var_slope)

• Now 1 of 51 states fall outside

Leave-one-out CV

loo_varslope <- loo_respondents(model_var_slope)
elpd_summary(loo_varslope)

elpd_loo       se 
-19036.7     36.3 

• Compare against the varying intercepts model

elpd_summary(loo_varslope - loo_varint)

elpd_loo       se 
    4.91     3.75 

• Some improvement - compared to how much the state intercepts helped, the slopes add relatively little predictive power
  • Most of what we learn about any one respondent comes from their state’s overall lean
  • But we do a lot better on the gap in education polarization by state.

Partial pooling

• Raw (no pooling) gaps (black) vs. the varying slopes model’s estimates (red)

model_gaps <- tibble(state = levels(as.factor(cells$state)), # model-implied gap in each state
                     model = colMeans(plogis(int_draws + slope_draws) - plogis(int_draws)))

• States with few respondents get pulled the most towards the national gap

Group-level predictors

• So far, every state’s coefficients are pooled towards the same national mean \(\mu_{\beta}\)

  • But we know a lot about the states themselves! (e.g. how they voted in 2012)

• Let the mean of each state’s coefficients depend on a vector of \(L\) state-level predictors \(Z_j\) (including a 1 for the intercept)

  \[\beta_j \sim \text{Normal}(\Gamma^{\prime}Z_j, \Sigma_{\beta})\]

  • \(\Gamma\) is an \(L \times K\) matrix: one column of coefficients for the state intercepts, one for the state slopes
  • Each state is pooled towards a prediction for that state rather than towards the national mean

• Gelman et al. (2007) model the state intercepts and income slopes using each state’s average income and region

  • We’ll use Obama’s 2012 share of the state’s two-party vote (on the logit scale) and the state’s share of adults with a college degree (both standardized), so \(L = 3\)

Group-level predictors

• As before, write each state’s coefficients as a mean plus a deviation: \(\beta_j = \Gamma^{\prime}Z_j + e_j\), with \(e_j \sim \text{Normal}(0, \Sigma_{\beta})\)

• Plug this into the linear predictor

  \[\begin{aligned}\log\bigg(\frac{\pi_{ij}}{1 - \pi_{ij}}\bigg) = &\underbrace{(\Gamma_{11} + \Gamma_{21}\text{obama}_j + \Gamma_{31}\text{coll}_j + e_{0j})}_{\beta_{0j}\text{: intercept}} \\ &+ \underbrace{(\Gamma_{12} + \Gamma_{22}\text{obama}_j + \Gamma_{32}\text{coll}_j + e_{1j})}_{\beta_{1j}\text{: education slope}} \times \text{college}_{ij}\end{aligned}\]

• The state-level predictors are just more regression terms

  • Main effects of \(\text{obama}_j\) and \(\text{coll}_j\) (shifting the intercepts), plus their interactions with \(\text{college}_{ij}\) (shifting the slopes)
  • Plus the same zero-mean state deviations \(e_j\) as in the varying slopes model

Group-level predictors

• Add the state-level predictors to the respondents and the cells, and rebuild X_fix with their main effects and interactions with college

state_preds <- read_csv("data/state_predictors.csv") %>% arrange(state)
state_preds_z <- state_preds %>% transmute(state, state_obama = as.numeric(scale(qlogis(obama_2012))),
                                           state_college = as.numeric(scale(college_share)))
X_fix <- model.matrix(~ college * (state_obama + state_college), data = voters %>% left_join(state_preds_z, by = "state"))
X_fix_cell <- model.matrix(~ college * (state_obama + state_college), data = cells %>% left_join(state_preds_z, by = "state"))
stan_data$P <- ncol(X_fix)
stan_data$X_fix <- X_fix
stan_data$X_fix_cell <- X_fix_cell
colnames(X_fix)

[1] "(Intercept)"           "college"               "state_obama"          
[4] "state_college"         "college:state_obama"   "college:state_college"

• And run the varying slopes code again

model_group <- stan(
  model_code = varying_slopes_model,  # same Stan code
  data = stan_data,       # with the new X_fix
  chains = 4,             # number of Markov chains
  warmup = 500,           # number of warmup iterations per chain
  iter = 2500,            # total number of iterations per chain
  cores = 4,
  refresh = 0,            # no progress shown
  seed = 53706,           # random seed
  control = list(adapt_delta = 0.95) # smaller steps (default is 0.8) - see next slide
  )

Group-level predictors

• Label the coefficients with the columns of X_fix

b_summary <- summary(model_group, pars = "b")$summary[, c("mean", "2.5%", "97.5%", "n_eff", "Rhat")]
rownames(b_summary) <- colnames(X_fix)
round(b_summary, 2)

                       mean  2.5% 97.5% n_eff Rhat
(Intercept)           -0.13 -0.17 -0.09  6031    1
college                0.69  0.62  0.76  5344    1
state_obama            0.32  0.24  0.40  4499    1
state_college          0.04 -0.03  0.10  4491    1
college:state_obama   -0.02 -0.16  0.12  4126    1
college:state_college  0.08 -0.04  0.19  4325    1

print(model_group, pars = c("tau", "Omega[1,2]"))

Inference for Stan model: anon_model.
4 chains, each with iter=2500; warmup=500; thin=1; 
post-warmup draws per chain=2000, total post-warmup draws=8000.

            mean se_mean   sd  2.5%   25%   50%  75% 97.5% n_eff Rhat
tau[1]      0.07    0.00 0.02  0.03  0.05  0.07 0.08  0.11   569 1.00
tau[2]      0.13    0.00 0.04  0.05  0.11  0.13 0.16  0.23   432 1.02
Omega[1,2] -0.07    0.01 0.35 -0.66 -0.33 -0.10 0.17  0.65   678 1.00

Samples were drawn using NUTS(diag_e) at Wed Oct 21 20:44:40 2026.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Group-level predictors

• state_obama: states that voted more for Obama in 2012 have higher intercepts
• college:state_college: states with more college graduates may have somewhat larger education gaps - but it’s uncertain

loo_group <- loo_respondents(model_group)
elpd_summary(loo_group - loo_varslope)

elpd_loo       se 
    18.7      5.3 

• A clear improvement over the varying slopes model without state-level predictors

Group-level predictors

• The state intercepts \(\beta_{0j}\) against 2012 vote: without (black) and with (red) the state-level predictors

• Small states like DC and Wyoming are no longer pulled towards the national mean, but towards states that voted like them

Multilevel Regression and Post-stratification

• Often we’re interested in estimating political attitudes at very small areas
  • Within a nationally-representative sample, we may only have a minor fraction of respondents who are even in a particular region.
• Beyond that, our sample may not be demographically representative for a given region
  • Can we do better than just using the sub-group means for our outcome of interest?

MrP

• Multilevel Regression and Post-stratification combines two well-known tools to try to accomplish this task (Park, Gelman, and Bafumi 2004)
  1. Multilevel regression: Using the entire dataset, specify a flexible model for the outcome of interest
  2. Post-stratification: Using known population characteristics, generate predictions for each covariate “cell” in each region
• The national-level regression allows us to improve precision in small-area estimation
  • e.g. if some of the variability between regions is driven by individual or group-level characteristics, we can model the relationship between those characteristics and the outcome using all the data.
• The population-level data allows us to address non-response/sampling bias
  • Ideally all of the variables that would go in the survey weighting model also go in the MrP model (so in principle we don’t need to include survey weights in the analysis)

MrP Workflow

1. Get survey data that contains individual and group-level predictors (e.g. demographic characteristics + “small area” of residence)
2. Get information on the “small areas” themselves
3. Estimate a multilevel regression model using these two sets of data
4. Construct a post-stratification frame that has the joint distribution of our demographics for each “small area”
5. Predict using the model for each cell in the post-stratification frame
6. Aggregate to get predictions for each small area

Post-stratification

• The easiest approach to adjusting a non-representative survey is to weight to match the known joint distribution of \(x\) in the population \(f_{\mathcal{U}}(\mathbf{x})\)
  • This requires a lot of auxiliary information about the target - typically obtained from high-quality census data
  • (e.g.) U.S. Census Public-Use Microdata: What is the share of Black, college educated, 30-45 year olds in Massachusetts?
• In post-stratification, we divide our sample up into \(C\) cells \(c\) that are mutually exclusive and exhaustive.
  • \(\tilde{T}_{\mathbf{x}} = \{\tilde{P}_{1}, \tilde{P}_2, \dotsc, \tilde{P}_{C}\}\) is our population distribution of these cells
• Our sampling/response model assumes that within each of these cells we have a simple-random sample from that particular stratum of the population
  • There may be variation in non-response or over/under-sampling, but only across cells
  • “post”-stratification because the intuition is akin to a design where we actually did stratify ex-ante
• Our measurement model assumes we observe the joint distribution of auxiliary variables \(\mathbf{x}\) in the target population
  • Difficult in many cases!

Post-stratification

• Two ways to think of post-stratification:

• First - Let \(\hat{\theta}_{y}^{c}\) denote our estimator for the population mean within cell \(c\)

  • Then, our post-stratification estimator \(\hat{\theta}_{y}^{PS}\) is:

  \[\hat{\theta}_{y}^{PS} = \sum_{c=1}^C \tilde{P}_{c} \hat{\theta}_{y}^{c}\]

• Alternatively, we’ll sometimes write \(\tilde{P}_c = \frac{\tilde{N}_c}{\tilde{N}}\) where \(\tilde{N}\) is the size of the population and \(\tilde{N}_c\) is the number of units in that cell.

  \[\hat{\theta}_{y}^{PS} = \sum_{c=1}^C \frac{\tilde{N}_c}{\tilde{N}} \hat{\theta}_{y}^{c}\]

• Fixed “constant” weights on the within-cell estimators.

  • This is how we’ll approach it with MRP.

Post-stratification

• Second - We can think of it as weighting individual observations using our adjustment weights \(\tilde{w}^{PS}_i\)

  • This is more common in survey weighting approaches to balancing sample and population.

• Let \(c(i)\) denote the class to which unit \(i\) belongs. Then the post-stratification weights are

  \[\tilde{w}^{PS}_i = (\tilde{P}_{c(i)}/\hat{P}^{S}_{c(i)}) \times d_i\]

  where \(\hat{P}^{S}_{c(i)}\) is the estimated in-sample proportion of observations in class \(c\) and \(d_i\) is unit \(i\)’s design weight

• When design weights are constant, \(\hat{P}^{S}_{c(i)}\) is just the sample mean of the indicator of class membership \(\mathbf{1}(i \in c)\)

  • More generally, you can write it as \(\hat{P}^{S}_{c(i)} = (\sum_{i=1}^N d_i \mathbf{1}(i \in c))/(\sum_{i=1}^N d_i)\)

• Intuition

  • The weights up-weight observations that are under-represented in the sample relative to the target population
  • The weights down-weight observations that are over-represented in the sample relative to the target population.

Example: CES 2016

• How well does the unweighted CES do at estimating Clinton’s share of the two-party vote in each state? (Actual county returns from MIT Election Data and Science Lab (2018))

pres <- read_csv("data/county_pres_2012_2016.csv")
state_results <- pres %>% group_by(state) %>% summarize(actual = sum(dem_2016)/sum(dem_2016 + rep_2016))
state_raw <- voters %>% group_by(state) %>% summarize(n = n(), raw = mean(clinton)) %>%
  inner_join(state_results, by = "state") # Alaska doesn't report results by county

Example: CES 2016

• Let’s fit a regression model
  • We’ll use the brms package which has pre-existing implementations of popular models like the random slopes/random intercepts generalized linear mixed model
• We’ll fit a logistic regression with coefficients on gender, age and college degree as well as our varying intercepts and education slopes by state
  • The brms formula syntax: (1 + college | state) is our varying slopes model from before

library(brms)

vote_cells <- voters %>% group_by(state, female, age, college) %>%
  summarize(n = n(), clinton = sum(clinton), .groups = "drop")

state_mlm <- brm(clinton | trials(n) ~ female + age + college + (1 + college | state),
                 family = binomial, data = vote_cells, cores = 4, seed = 53706, refresh = 0)

Model diagnostics

• Summarize the results

summary(state_mlm)$random$state[, 1:5] # state-level standard deviations and correlation

                       Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)             0.295    0.0445   0.2190    0.391    1
sd(college)               0.143    0.0494   0.0443    0.244    1
cor(Intercept,college)    0.257    0.3088  -0.3445    0.845    1

summary(state_mlm)$fixed[, 1:5]

          Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept    0.507    0.0882    0.331    0.676    1
female       0.411    0.0253    0.361    0.461    1
age25M34    -0.380    0.0859   -0.546   -0.212    1
age35M44    -0.647    0.0832   -0.806   -0.484    1
age45M64    -0.939    0.0763   -1.088   -0.790    1
age65P      -1.096    0.0781   -1.245   -0.947    1
college      0.679    0.0365    0.608    0.752    1

Who is the target population?

• Our post-stratification frame comes from the Census - it counts adults
  • But our target is Clinton’s share among voters
  • Turnout differs a lot across cells (e.g. by age and education)
• So we need a second model: the probability that someone in each cell turns out
  • We can fit it on the whole CES sample, voters and non-voters (validated turnout)
  • Then weight each cell by \(\tilde{N}_c \times \Pr(\text{vote} | c)\)

turnout_cells <- ces %>% group_by(state, female, age, college) %>%
  summarize(n = n(), voted = sum(voted), .groups = "drop")

turnout_mlm <- brm(voted | trials(n) ~ female + age + college + (1 | state),
                   family = binomial, data = turnout_cells, cores = 4, seed = 53706, refresh = 0)

Generating the post-stratification frame

• We’ll use the ACS 2012-2016 5-year estimates of sex \(\times\) age \(\times\) education among adults in each county (table B15001)

acs <- read_csv("data/acs2016_county_frame.csv") %>%
  mutate(college = as.integer(educ %in% c("BA", "Postgrad")))

• Each row contains a row for each county x gender x age x education combination

Predicting with the model

• Our models only depend on state and demographics: predict once for each state x gender x age x college degree cell

state_frame <- acs %>% group_by(state, female, age, college) %>%
  summarize(N = sum(n), .groups = "drop") %>% mutate(n = 1) # n = 1: predict a probability
state_frame$p_turnout <- colMeans(posterior_epred(turnout_mlm, newdata = state_frame))
state_frame$p_clinton <- colMeans(posterior_epred(state_mlm, newdata = state_frame))

• Aggregate the predictions by the stratum shares (adults \(\times\) turnout) to get state-level estimates!

mrp_states <- state_frame %>% group_by(state) %>%
  summarize(mrp = sum(N * p_turnout * p_clinton)/sum(N * p_turnout)) %>%
  inner_join(state_raw, by = "state")

Compare

mrp_states %>% summarize(rmse_raw = sqrt(mean((raw - actual)^2)), rmse_mrp = sqrt(mean((mrp - actual)^2)),
                         bias_raw = mean(raw - actual), bias_mrp = mean(mrp - actual))

# A tibble: 1 × 4
  rmse_raw rmse_mrp bias_raw bias_mrp
     <dbl>    <dbl>    <dbl>    <dbl>
1   0.0774   0.0865   0.0518   0.0488

Pitfall: Pooling towards the national mean

• MrP barely reduces the average bias, and its RMSE is worse than the raw means

  • The state intercepts are partially pooled towards a common grand mean
  • So the most extreme states get pulled towards the middle

• DC: raw 97%, MrP 78%, actual 96%

• Wyoming: raw 48%, MrP 50%, actual 24%

• But we know there are some extreme states out there!

  • The solution is to include state-level predictors (e.g. past vote share), so that each state is pooled towards a prediction for that state, not towards the nation as a whole (Park, Gelman, and Bafumi 2004)

Group-level predictors in brms

• In brms (and lme4), a group-level predictor is just a population-level term: a column that’s constant within each group

  y ~ college * (state_obama + state_college) + (1 + college | state)

• This is exactly our model from before, (though brms has slightly different default priors)
  • Mechanically, brms also avoids divergent transitions by not sampling \(\tau\) directly (a “non-centered” parameterization)

Group-level predictors in brms

• Only the interactions move the slopes: main effects of the state-level predictors shift the intercepts

• Every varying term needs a population-level counterpart

  • (1 + college | state) centers the state slopes at zero, so college has to be in the population-level part too
  • Otherwise we’re assuming that the average education gap is zero

• Merge the group-level predictors onto every row of the data (they must be constant within each group - brms won’t check!)

• This works with varying intercepts, but not with state fixed effects

  • With a dummy for each state, state_obama would be perfectly collinear with the dummies
  • With partially pooled intercepts, its coefficient is learned from the differences between states

Group-level predictors in brms

• brms can also write the model literally as \(\beta_j = \Gamma^{\prime}Z_j + e_j\) with its bf() formula syntax

bf(clinton | trials(n) ~ b0 + b1 * college,
   b0 ~ state_obama + state_college + (1 | s | state), # intercepts: Gamma[, 1] + e_0j
   b1 ~ state_obama + state_college + (1 | s | state), # slopes: Gamma[, 2] + e_1j
   nl = TRUE) # the shared "s" lets e_0j and e_1j be correlated

• Here we’ll stick with the formula version. Add the state-level predictors to our cells

vote_cells_st <- vote_cells %>% left_join(state_preds_z, by = "state")

state_mlm2 <- brm(clinton | trials(n) ~ female + age + college * (state_obama + state_college) +
                    (1 + college | state),
                  family = binomial, data = vote_cells_st, cores = 4, seed = 53706, refresh = 0)

State-level predictors

• Predict and post-stratify as before

state_frame <- state_frame %>% left_join(state_preds_z, by = "state")
state_frame$p_clinton2 <- colMeans(posterior_epred(state_mlm2, newdata = state_frame))
mrp_states <- state_frame %>% group_by(state) %>%
  summarize(mrp_statepred = sum(N * p_turnout * p_clinton2)/sum(N * p_turnout)) %>%
  inner_join(mrp_states, by = "state")
mrp_states %>% pivot_longer(c(raw, mrp, mrp_statepred), names_to = "estimate") %>% group_by(estimate) %>%
  summarize(rmse = sqrt(mean((value - actual)^2)), bias = mean(value - actual))

# A tibble: 3 × 3
  estimate        rmse   bias
  <chr>          <dbl>  <dbl>
1 mrp           0.0865 0.0488
2 mrp_statepred 0.0681 0.0485
3 raw           0.0774 0.0518

• DC is now at 88% and Wyoming at 39%: pooled towards states that voted like them

Smaller areas

• MrP’s value is typically in making preductions for small areas where we have few or no respondents

voters %>% filter(county_fips %in% pres$county_fips) %>% count(county_fips) %>%
  summarize(counties_with_respondents = n(), median_respondents = median(n))

# A tibble: 1 × 2
  counties_with_respondents median_respondents
                      <int>              <dbl>
1                      2290                  4

nrow(pres) # counties (outside Alaska)

[1] 3112

• We’ll analyze at the county level now (our frame was already county-level) - let’s combine the predictions (from the model with state-level predictors).

mrp_counties <- acs %>% left_join(state_frame %>% select(state, female, age, college, p_turnout, p_clinton2),
                                  by = c("state", "female", "age", "college")) %>%
  group_by(county_fips) %>% summarize(mrp_state = sum(n * p_turnout * p_clinton2)/sum(n * p_turnout)) %>%
  inner_join(pres %>% select(county_fips, state, dem_share_2016), by = "county_fips")

Visualize

Visualize

• Zoom in on the Midwest

Pitfall: State lines

• In the model, the only thing that distinguishes two counties in the same state is their demographics (gender, age, education)
  • Everything else is the state effect
  • Rural Illinois gets the same state effect as Chicago, so it looks far more Democratic than neighboring Indiana or Missouri
• Estimates jump at state boundaries
  • Even though there is likely a lot more smoothness across the border
• Common feature of small-area MRP estimates!

Pitfall: State lines

• Office of the Assistant Secretary for Planning and Evaluation (2021): county estimates of COVID-19 vaccine hesitancy, from a logistic regression on demographics + state applied to Census microdata

County-level predictors

• The solution is to give the model information about the small areas themselves

  • A county-level predictor: Obama’s 2012 share of the two-party vote in the county (on the logit scale)
  • A varying intercept for each county, partially pooled towards the prediction from the county-level predictor

• In brms, a group-level predictor is just another column that’s constant within each group: obama_2012 + (1 | county_fips) pools each county’s intercept towards \(\gamma_0 + \gamma_1 \times \text{(2012 vote)}_c\)

• The county’s own 2012 vote takes the place of the state-level predictors

  • Counties with no respondents still get a prediction from their demographics, their state, and their 2012 vote

county_info <- pres %>% transmute(county_fips, obama_2012 = qlogis(dem_share_2012))
county_cells <- voters %>% inner_join(county_info, by = "county_fips") %>%
  group_by(state, county_fips, obama_2012, female, age, college) %>%
  summarize(n = n(), clinton = sum(clinton), .groups = "drop")

county_mlm <- brm(clinton | trials(n) ~ female + age + college + obama_2012 +
                    (1 + college | state) + (1 | county_fips),
                  family = binomial, data = county_cells, cores = 4, seed = 53706, refresh = 0)

County-level predictors

summary(county_mlm)$random$county_fips[, 1:5] # county-level standard deviation

              Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)    0.153    0.0288   0.0965    0.208    1

summary(county_mlm)$random$state[, 1:5]

                       Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)            0.0611    0.0291  0.00666    0.118    1
sd(college)              0.1042    0.0501  0.01114    0.206    1
cor(Intercept,college)  -0.1938    0.4854 -0.90368    0.865    1

summary(county_mlm)$fixed[, 1:5]

           Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept     0.470    0.0816    0.308    0.630    1
female        0.434    0.0265    0.382    0.485    1
age25M34     -0.374    0.0875   -0.545   -0.204    1
age35M44     -0.628    0.0868   -0.793   -0.452    1
age45M64     -0.901    0.0799   -1.058   -0.744    1
age65P       -1.047    0.0814   -1.212   -0.889    1
college       0.608    0.0340    0.542    0.678    1
obama_2012    0.836    0.0247    0.789    0.884    1

Predicting with the model

• Now predictions depend on the county: predict for every row of the county frame
  • allow_new_levels: counties with no respondents draw their county effect from the fitted distribution

county_frame <- acs %>% inner_join(county_info, by = "county_fips") %>% mutate(n_adults = n, n = 1)
county_frame$p_clinton <- colMeans(posterior_epred(county_mlm, newdata = county_frame,
                                                   allow_new_levels = TRUE, ndraws = 500, seed = 53706))

county_frame <- county_frame %>% left_join(state_frame %>% select(state, female, age, college, p_turnout),
                                           by = c("state", "female", "age", "college"))
mrp_counties <- mrp_counties %>% left_join(county_frame %>% group_by(county_fips) %>%
  summarize(mrp_county = sum(n_adults * p_turnout * p_clinton)/sum(n_adults * p_turnout)), by = "county_fips")

Visualize

Visualize

Compare

• How close are the county estimates to the actual results?

mrp_counties %>% pivot_longer(c(mrp_state, mrp_county), names_to = "model", values_to = "estimate") %>%
  group_by(model) %>% summarize(rmse = sqrt(mean((estimate - dem_share_2016)^2)),
                                bias = mean(estimate - dem_share_2016),
                                correlation = cor(estimate, dem_share_2016))

# A tibble: 2 × 4
  model        rmse   bias correlation
  <chr>       <dbl>  <dbl>       <dbl>
1 mrp_county 0.0979 0.0848       0.964
2 mrp_state  0.205  0.153        0.555

• Aggregate the county predictions up to the states too

mrp_states <- county_frame %>% group_by(state) %>%
  summarize(mrp_county = sum(n_adults * p_turnout * p_clinton)/sum(n_adults * p_turnout)) %>%
  inner_join(mrp_states, by = "state")
mrp_states %>% pivot_longer(c(raw, mrp, mrp_statepred, mrp_county), names_to = "estimate") %>% group_by(estimate) %>%
  summarize(rmse = sqrt(mean((value - actual)^2)), bias = mean(value - actual))

# A tibble: 4 × 3
  estimate        rmse   bias
  <chr>          <dbl>  <dbl>
1 mrp           0.0865 0.0488
2 mrp_county    0.0595 0.0461
3 mrp_statepred 0.0681 0.0485
4 raw           0.0774 0.0518

MrP can only fix what’s in the model

• Even with the county predictor, our estimates still lean towards Clinton
  • Clinton’s estimated share is 4.6 points too high in the average state, and 8.5 points too high in the average county
• The sampling/response model assumption: within each cell, respondents look like non-respondents
  • CES respondents are more Democratic than people with the same gender, age, education and county who didn’t take the survey

How to improve this model

• We’ve provided a very simple example of an MrP model but there are a lot of ways it can be improved
  1. More group-level predictors (e.g. county demographics, urbanicity)
  2. Individual-level interactions in the model + partial-pooling on the coefficients
  3. More covariates + knowing the true joint distribution in the population.
  4. Many more possibilities…
• MrP is most powerful when there are very few respondents from a particular small area in our data
  • But it can only correct for differences between the sample and the population that are in the model

Next week

• Mixture models
  • Exploratory data analysis and clustering models
  • MLE and MAP estimation via the “Expectation-Maximization” algorithm
• Readings: Chapters 13 and 22 of Gelman et al. (2013)