library(haven) # read .dta files
library(arm) # sim()
library(lme4) # glmer()
library(rstanarm) # stan_glmer()
library(bayesplot)
library(ggplot2)
library(dplyr)
library(tidyr)
library(purrr)
library(knitr)
library(skimr)
options(mc.cores = min(4, parallel::detectCores()))
theme_set(bayesplot::theme_default(base_size = 11))Part 4 led us into Bayesian regression with MCMC. We used rstanarm to fit logistic, probit, and ordered models, and the main lesson was computational: once the posterior is no longer available in closed form, we sample from it.
This post moves from single-level regression to multilevel modeling. The main idea is: when observations are nested within groups, we should recognize group differences without treating each group as entirely separate.
Multilevel models give us a principled compromise called partial pooling. Multilevle or hierarchical models provide a principled approach to estimating group-specific effects through empirical Bayes, or shrinkage, estimation: each group’s estimate is informed by its own data while also being adjusted toward the value predicted by the higher-level model.
We will use the 1992 American National Election Study (NES) to predict whether a voter supported George H. W. Bush. Voters are nested in states, so the question is:
How do we estimate state differences in Bush support without overreacting to noisy state-level samples?
Along the way, we will compare:
- Complete pooling: ignore state entirely.
- No pooling: estimate a separate fixed effect for every state.
- Partial pooling: estimate state effects as draws from a shared distribution.
- Bayesian partial pooling: fit the same multilevel model with
stan_glmer()and summarize posterior uncertainty directly.
Setup
1. The Pooling Problem
Suppose voters are nested in states. Let \(y_i = 1\) if voter \(i\) voted for Bush and \(0\) otherwise.
A simple logistic regression ignores state:
\[ \Pr(y_i = 1) = \text{logit}^{-1}(X_i \beta). \]
This is complete pooling. Every voter contributes to one national regression line. The model is stable, but it cannot show residual state-to-state differences.
The opposite approach estimates a separate intercept for every state:
\[ \Pr(y_i = 1) = \text{logit}^{-1}(\alpha_{s[i]} + X_i \beta). \]
This is no pooling. Each state gets its own intercept. The model is flexible, but small states can get noisy, extreme estimates.
The multilevel model lands between those extremes:
\[ \Pr(y_i = 1) = \text{logit}^{-1}(\alpha + a_{s[i]} + X_i \beta), \qquad a_s \sim \mathcal{N}(0, \sigma_a^2). \]
This is what we called partial pooling. State effects are estimated as coming from a shared distribution. States with little data are pulled toward the national average; states with more data are allowed to move farther away.
NoteUnderstaning partial pooling in education context
Example: Suppose you want to estimate each school’s average achievement. A school with only 10 students should not be allowed to look extremely high or low based only on a tiny sample. Partial pooling says: “Use what we know from this school, but also borrow strength from all schools.” So the small school’s estimate is pulled more toward the overall pattern, while a large school’s estimate stays closer to its own observed average.
In Raudenbush and Bryk’s (2002) terminology, this is closely related to empirical Bayes shrinkage: group-specific estimates are shrunk more when they are less reliable, especially when the group has a small sample size or more noise.
2. Why This Is Bayesian in Spirit
The random-effect assumption
\[ a_s \sim \mathcal{N}(0, \sigma_a^2) \]
resembles a prior distribution on state effects. In a fully Bayesian model, it would be interpreted in that way. In glmer(), this structure is estimated using maximum likelihood rather than full Bayesian inference, but it still carries a similar regularizing logic.
In this sense, even before moving to Stan, multilevel modeling already has a Bayesian flavor:
- Group effects are treated as uncertain quantities.
- Extreme estimates are moderated toward a shared center.
- Predictions combine individual-level information with group-level information.
- Uncertainty can be explored by simulating from the fitted model.
This is why multilevel modeling offers a natural bridge from classical regression to Bayesian data analysis.
3. NES 1992 Data
The NES is a long-running survey of U.S. voters. We use the 1992 presidential election, keep only two-party voters, and code:
\[ y_i = \begin{cases} 1, & \text{voted for Bush} \\ 0, & \text{voted for Clinton.} \end{cases} \]
nes_raw <- read_dta(
"https://faculty.stat.ucla.edu/handcock/216/examples/nes/nes5200_processed_voters_realideo.dta"
)
scale_2sd <- function(x) {
as.numeric((x - mean(x, na.rm = TRUE)) / (2 * sd(x, na.rm = TRUE)))
}
nes92 <- nes_raw %>%
filter(year == 1992, presvote < 3) %>%
mutate(vote_bush = ifelse(presvote == 2, 1, 0)) %>%
select(vote_bush, female, black, age, educ1, income,
partyid7, real_ideo, ideo_feel, state) %>%
na.omit() %>%
mutate(
age_c = scale_2sd(age),
educ_c = scale_2sd(educ1),
income_c = scale_2sd(income),
party_c = scale_2sd(partyid7),
ideo_c = scale_2sd(real_ideo),
state = factor(state)
)
dim(nes92)[1] 921 15skim(nes92)| Name | nes92 |
| Number of rows | 921 |
| Number of columns | 15 |
| _______________________ | |
| Column type frequency: | |
| factor | 1 |
| numeric | 14 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| state | 0 | 1 | FALSE | 33 | 71: 100, 13: 57, 23: 55, 40: 53 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| vote_bush | 0 | 1 | 0.44 | 0.50 | 0.00 | 0.00 | 0.00 | 1.00 | 1.00 | ▇▁▁▁▆ |
| female | 0 | 1 | 0.52 | 0.50 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00 | ▇▁▁▁▇ |
| black | 0 | 1 | 0.10 | 0.30 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| age | 0 | 1 | 46.39 | 16.81 | 18.00 | 33.00 | 43.00 | 59.00 | 91.00 | ▅▇▅▃▁ |
| educ1 | 0 | 1 | 3.00 | 0.92 | 1.00 | 2.00 | 3.00 | 4.00 | 4.00 | ▁▆▁▅▇ |
| income | 0 | 1 | 3.21 | 1.05 | 1.00 | 3.00 | 3.00 | 4.00 | 5.00 | ▂▃▇▇▂ |
| partyid7 | 0 | 1 | 3.76 | 2.23 | 1.00 | 2.00 | 3.00 | 6.00 | 7.00 | ▇▃▁▂▆ |
| real_ideo | 0 | 1 | 4.16 | 1.48 | 1.00 | 3.00 | 4.00 | 5.00 | 7.00 | ▅▃▇▆▆ |
| ideo_feel | 0 | 1 | 51.44 | 18.14 | 0.00 | 39.00 | 49.00 | 61.00 | 97.00 | ▁▂▇▃▂ |
| age_c | 0 | 1 | 0.00 | 0.50 | -0.84 | -0.40 | -0.10 | 0.38 | 1.33 | ▅▇▅▃▁ |
| educ_c | 0 | 1 | 0.00 | 0.50 | -1.08 | -0.54 | 0.00 | 0.54 | 0.54 | ▁▆▁▅▇ |
| income_c | 0 | 1 | 0.00 | 0.50 | -1.05 | -0.10 | -0.10 | 0.38 | 0.85 | ▂▃▇▇▂ |
| party_c | 0 | 1 | 0.00 | 0.50 | -0.62 | -0.39 | -0.17 | 0.50 | 0.73 | ▇▃▁▂▆ |
| ideo_c | 0 | 1 | 0.00 | 0.50 | -1.07 | -0.39 | -0.05 | 0.28 | 0.96 | ▅▃▇▆▆ |
NoteWhy divide continuous predictors by 2 SD?
Binary predictors, such as female or black, usually change from 0 to 1. This makes their coefficients easy to interpret: they describe the difference between the two groups.
Continuous predictors, such as age or income, are different. A one-unit change may be very small, such as one year of age or one dollar of income. This can make their coefficients difficult to compare with binary predictors.
Gelman (2008) recommends dividing continuous predictors by two standard deviations. This puts continuous predictors on a scale that is more comparable to a 0-to-1 change in a binary variable.
For example, if we divide age by 2 * sd(age), the age coefficient no longer describes the effect of a one-year increase. Instead, it describes a larger and more meaningful difference, roughly comparing a relatively young voter with a relatively old voter.
This scaling makes it easier to compare which predictors have larger or smaller associations with the outcome.
4. Complete Pooling: Ignore State
Start with an ordinary logistic regression. This model uses individual-level predictors but ignores state.
fit_complete <- glm(
vote_bush ~ female + black + age_c + educ_c + income_c +
party_c + ideo_c + income_c:ideo_c + party_c:ideo_c,
family = binomial,
data = nes92
)
summary(fit_complete)
Call:
glm(formula = vote_bush ~ female + black + age_c + educ_c + income_c +
party_c + ideo_c + income_c:ideo_c + party_c:ideo_c, family = binomial,
data = nes92)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.5443 0.1786 -3.047 0.00231 **
female 0.2212 0.2406 0.919 0.35802
black -2.5920 0.6572 -3.944 8.01e-05 ***
age_c 0.0863 0.2401 0.359 0.71930
educ_c 0.2485 0.2687 0.925 0.35513
income_c -0.2088 0.2559 -0.816 0.41453
party_c 4.7002 0.3538 13.286 < 2e-16 ***
ideo_c 2.3310 0.3142 7.418 1.19e-13 ***
income_c:ideo_c 0.3552 0.5860 0.606 0.54440
party_c:ideo_c -1.4953 0.8110 -1.844 0.06521 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1261.87 on 920 degrees of freedom
Residual deviance: 485.73 on 911 degrees of freedom
AIC: 505.73
Number of Fisher Scoring iterations: 6exp(coef(fit_complete)[c("party_c", "ideo_c", "black")]) party_c ideo_c black
109.96477198 10.28767247 0.07487281 On the logit scale, positive coefficients indicate higher log-odds of voting for Bush. On the odds scale, exp(beta) represents the multiplicative change in the odds of voting for Bush.
The main findings are straightforward:
- Voters who identify as more Republican are much more likely to vote for Bush.
- Voters who identify as more conservative are also more likely to vote for Bush.
- Black respondents are much less likely to vote for Bush, after accounting for the other predictors.
However, this model cannot tell us whether some states remain more or less pro-Bush after adjusting for individual-level characteristics.
5. No Pooling: Fixed State Effects
The no-pooling model includes state as a fixed effect. This gives every state its own intercept.
fit_no_pool <- glm(
vote_bush ~ female + black + age_c + educ_c + income_c +
party_c + ideo_c + income_c:ideo_c + party_c:ideo_c +
state,
family = binomial,
data = nes92
)
summary(fit_no_pool)
Call:
glm(formula = vote_bush ~ female + black + age_c + educ_c + income_c +
party_c + ideo_c + income_c:ideo_c + party_c:ideo_c + state,
family = binomial, data = nes92)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.93438 0.84076 -2.301 0.021405 *
female 0.06659 0.25712 0.259 0.795639
black -2.60811 0.69436 -3.756 0.000173 ***
age_c -0.01770 0.26463 -0.067 0.946670
educ_c 0.23975 0.29509 0.812 0.416528
income_c -0.19912 0.28266 -0.704 0.481149
party_c 4.90394 0.39172 12.519 < 2e-16 ***
ideo_c 2.45450 0.33930 7.234 4.69e-13 ***
state3 1.58701 1.01331 1.566 0.117309
state4 0.05214 1.35517 0.038 0.969311
state12 0.75188 1.06012 0.709 0.478178
state13 1.25901 0.95611 1.317 0.187904
state14 0.94179 1.06080 0.888 0.374641
state21 2.42057 0.98959 2.446 0.014444 *
state22 2.63188 1.10135 2.390 0.016863 *
state23 1.11641 0.95502 1.169 0.242410
state24 1.84597 1.05719 1.746 0.080792 .
state25 1.08388 1.23581 0.877 0.380451
state31 0.59518 1.62224 0.367 0.713705
state32 1.44053 1.15704 1.245 0.213126
state33 0.08732 0.98124 0.089 0.929090
state34 2.33381 1.14244 2.043 0.041069 *
state35 3.06362 2.37235 1.291 0.196569
state40 1.21255 1.00525 1.206 0.227733
state41 1.50154 1.30427 1.151 0.249627
state42 0.91778 1.54124 0.595 0.551520
state43 1.92289 1.06020 1.814 0.069722 .
state44 2.01590 1.05558 1.910 0.056165 .
state45 0.29964 1.83012 0.164 0.869945
state47 2.65722 1.45308 1.829 0.067447 .
state49 1.78729 1.01894 1.754 0.079419 .
state52 1.80362 1.09948 1.640 0.100915
state54 2.18485 1.00441 2.175 0.029611 *
state56 2.70501 1.21284 2.230 0.025726 *
state61 1.68195 1.08147 1.555 0.119889
state62 -0.97162 1.55662 -0.624 0.532507
state68 2.38610 1.60866 1.483 0.138000
state71 1.18694 0.93217 1.273 0.202907
state72 1.58790 1.22167 1.300 0.193676
state73 1.42856 1.23015 1.161 0.245524
income_c:ideo_c 0.53192 0.63830 0.833 0.404652
party_c:ideo_c -1.37914 0.86683 -1.591 0.111606
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1261.87 on 920 degrees of freedom
Residual deviance: 451.59 on 879 degrees of freedom
AIC: 535.59
Number of Fisher Scoring iterations: 7The problem is not that state fixed effects are wrong. They are often useful. The problem is that they treat every state’s estimate as independent. A state with a small sample can get an extreme intercept even if the evidence is thin.
Let’s extract the state fixed effects and compare them to state sample sizes.
state_counts <- nes92 %>%
count(state, name = "n")
no_pool_effects <- coef(fit_no_pool) %>%
as.data.frame() %>%
tibble::rownames_to_column("term") %>%
rename(estimate = 2) %>%
filter(grepl("^state", term)) %>%
mutate(state = sub("^state", "", term)) %>%
select(state, estimate)
ggplot(no_pool_effects %>% left_join(state_counts, by = "state"),
aes(x = n, y = estimate)) +
geom_point(alpha = 0.75) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
labs(
x = "State sample size",
y = "No-pooling state fixed effect",
title = "Small states can have noisy fixed-effect estimates"
) +
theme_bw()
WarningA subtle issue with state fixed effects
Fixed effects do not know that state estimates are related to one another. If a state has little data, the model does not automatically borrow strength from other states. That is exactly what multilevel modeling fixes.
6. Partial Pooling: Varying Intercepts by State
Now fit the multilevel model. The formula (1 | state) means each state gets its own intercept, and those intercepts are modeled as coming from a shared normal distribution.
fit_partial <- glmer(
vote_bush ~ female + black + age_c + educ_c + income_c +
party_c + ideo_c + income_c:ideo_c + party_c:ideo_c +
(1 | state),
family = binomial,
data = nes92
)
summary(fit_partial)Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: vote_bush ~ female + black + age_c + educ_c + income_c + party_c +
ideo_c + income_c:ideo_c + party_c:ideo_c + (1 | state)
Data: nes92
AIC BIC logLik -2*log(L) df.resid
507.7 560.7 -242.8 485.7 910
Scaled residuals:
Min 1Q Median 3Q Max
-8.6838 -0.2273 -0.0516 0.2516 6.1013
Random effects:
Groups Name Variance Std.Dev.
state (Intercept) 0.02963 0.1721
Number of obs: 921, groups: state, 33
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.5386 0.1837 -2.932 0.00337 **
female 0.2120 0.2436 0.870 0.38404
black -2.6025 0.6600 -3.943 8.04e-05 ***
age_c 0.0750 0.2441 0.307 0.75870
educ_c 0.2460 0.2704 0.910 0.36277
income_c -0.2093 0.2571 -0.814 0.41567
party_c 4.7131 0.3583 13.153 < 2e-16 ***
ideo_c 2.3429 0.3185 7.355 1.91e-13 ***
income_c:ideo_c 0.3547 0.5882 0.603 0.54650
party_c:ideo_c -1.4880 0.8150 -1.826 0.06787 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) female black age_c educ_c incm_c prty_c ideo_c inc_:_
female -0.698
black -0.085 -0.014
age_c 0.048 -0.042 0.095
educ_c 0.008 0.059 -0.011 0.164
income_c -0.014 0.026 0.017 0.022 -0.255
party_c -0.156 0.066 -0.032 0.056 -0.024 -0.124
ideo_c -0.169 0.006 -0.135 -0.100 0.164 -0.064 0.033
incom_c:d_c -0.070 0.070 0.054 0.003 -0.086 -0.163 -0.004 0.124
party_c:d_c 0.031 0.005 -0.027 -0.030 -0.118 0.014 -0.418 -0.098 -0.166The state-level standard deviation tells us how much residual state-to-state variation remains after accounting for individual predictors.
as.data.frame(VarCorr(fit_partial)) %>%
select(grp, var1, sdcor) grp var1 sdcor
1 state (Intercept) 0.172135If this standard deviation is small, then most of the state differences are already explained by individual-level covariates. If it is large, similar voters still differ meaningfully across states.
7. Seeing Shrinkage
Partial pooling shrinks state effects toward zero. The amount of shrinkage is not arbitrary: it depends on state sample size and estimated state-level variance.
partial_effects <- ranef(fit_partial)$state[, 1] %>%
as.data.frame() %>%
tibble::rownames_to_column("state") %>%
rename(partial_pooling = 2)
shrink_df <- no_pool_effects %>%
rename(no_pooling = estimate) %>%
inner_join(partial_effects, by = "state") %>%
left_join(state_counts, by = "state")
ggplot(shrink_df, aes(x = no_pooling, y = partial_pooling, size = n)) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray60") +
geom_point(alpha = 0.75, color = "#5B8DB8") +
scale_size_continuous(range = c(2, 7)) +
labs(
x = "No-pooling state effect",
y = "Partial-pooling state effect",
size = "State n",
title = "Partial pooling pulls state effects toward the national average"
) +
theme_bw()
Points below the diagonal in absolute value have been pulled toward zero. Smaller states usually shrink more because the model has less state-specific information for them.
The gray dash line is the “no change” reference line, where partial-pooling state effect = no-pooling state effect. So:
- Points on the dashed line: partial pooling gives the same estimate as no pooling.
- Points below the dashed line: partial pooling gives a smaller estimate than no pooling.
- Points above the dashed line: partial pooling gives a larger estimate than no pooling.
8. Bayesian Partial Pooling with stan_glmer
Now we can fit a Bayesian version of the same multilevel logistic regression model. The model structure is very similar to the glmer() model: we still predict whether a respondent voted for Bush, include individual-level predictors, and allow the intercept to vary by state.
The key difference is that stan_glmer() estimates a posterior distribution for each parameter rather than only producing one point estimate and a standard error. This means we can describe not only the most likely value of each coefficient, but also the uncertainty around it.
fit_stan <- stan_glmer(
vote_bush ~ female + black + age_c + educ_c + income_c +
party_c + ideo_c + income_c:ideo_c + party_c:ideo_c +
(1 | state),
family = binomial(link = "logit"),
data = nes92,
prior = normal(0, 2.5),
prior_intercept = normal(0, 5),
iter = 2000,
chains = 4,
seed = 216,
refresh = 0
)
summary(fit_stan, probs = c(0.025, 0.5, 0.975), digits = 2)
Model Info:
function: stan_glmer
family: binomial [logit]
formula: vote_bush ~ female + black + age_c + educ_c + income_c + party_c +
ideo_c + income_c:ideo_c + party_c:ideo_c + (1 | state)
algorithm: sampling
sample: 4000 (posterior sample size)
priors: see help('prior_summary')
observations: 921
groups: state (33)
Estimates:
mean sd 2.5% 50% 97.5%
(Intercept) -0.54 0.19 -0.91 -0.53 -0.18
female 0.20 0.24 -0.26 0.20 0.65
black -2.57 0.64 -3.88 -2.53 -1.41
age_c 0.07 0.24 -0.39 0.07 0.54
educ_c 0.23 0.27 -0.30 0.23 0.79
income_c -0.20 0.26 -0.69 -0.20 0.31
party_c 4.73 0.35 4.06 4.72 5.45
ideo_c 2.36 0.32 1.76 2.36 3.02
income_c:ideo_c 0.34 0.58 -0.79 0.34 1.47
party_c:ideo_c -1.28 0.77 -2.82 -1.26 0.21
b[(Intercept) state:1] -0.16 0.31 -0.94 -0.08 0.30
b[(Intercept) state:3] 0.01 0.24 -0.50 0.00 0.54
b[(Intercept) state:4] -0.12 0.30 -0.87 -0.05 0.37
b[(Intercept) state:12] -0.10 0.28 -0.77 -0.05 0.37
b[(Intercept) state:13] -0.04 0.23 -0.57 -0.02 0.39
b[(Intercept) state:14] -0.08 0.27 -0.73 -0.03 0.41
b[(Intercept) state:21] 0.19 0.29 -0.23 0.11 0.89
b[(Intercept) state:22] 0.17 0.30 -0.29 0.08 0.91
b[(Intercept) state:23] -0.08 0.24 -0.65 -0.03 0.37
b[(Intercept) state:24] 0.06 0.27 -0.43 0.02 0.70
b[(Intercept) state:25] -0.02 0.27 -0.64 0.00 0.52
b[(Intercept) state:31] -0.04 0.29 -0.71 -0.01 0.53
b[(Intercept) state:32] 0.01 0.27 -0.56 0.00 0.60
b[(Intercept) state:33] -0.27 0.33 -1.10 -0.17 0.13
b[(Intercept) state:34] 0.09 0.28 -0.41 0.04 0.76
b[(Intercept) state:35] 0.04 0.30 -0.56 0.01 0.74
b[(Intercept) state:40] -0.03 0.25 -0.62 -0.01 0.47
b[(Intercept) state:41] 0.01 0.27 -0.57 0.00 0.61
b[(Intercept) state:42] -0.03 0.29 -0.71 -0.01 0.56
b[(Intercept) state:43] 0.07 0.25 -0.40 0.02 0.67
b[(Intercept) state:44] 0.10 0.27 -0.36 0.04 0.78
b[(Intercept) state:45] -0.03 0.30 -0.71 -0.01 0.56
b[(Intercept) state:47] 0.07 0.30 -0.47 0.02 0.79
b[(Intercept) state:49] 0.06 0.26 -0.45 0.02 0.65
b[(Intercept) state:52] 0.04 0.26 -0.49 0.01 0.62
b[(Intercept) state:54] 0.14 0.28 -0.29 0.07 0.84
b[(Intercept) state:56] 0.10 0.30 -0.43 0.04 0.85
b[(Intercept) state:61] 0.02 0.26 -0.50 0.01 0.62
b[(Intercept) state:62] -0.15 0.32 -0.95 -0.07 0.34
b[(Intercept) state:68] 0.04 0.30 -0.55 0.01 0.70
b[(Intercept) state:71] -0.07 0.23 -0.62 -0.03 0.35
b[(Intercept) state:72] 0.01 0.27 -0.57 0.00 0.61
b[(Intercept) state:73] -0.01 0.27 -0.63 -0.01 0.57
Sigma[state:(Intercept),(Intercept)] 0.09 0.11 0.00 0.05 0.40
Fit Diagnostics:
mean sd 2.5% 50% 97.5%
mean_PPD 0.44 0.01 0.41 0.44 0.46
The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
MCMC diagnostics
mcse Rhat n_eff
(Intercept) 0.00 1.00 4436
female 0.00 1.00 5291
black 0.01 1.00 4317
age_c 0.00 1.00 3434
educ_c 0.00 1.00 4916
income_c 0.00 1.00 4682
party_c 0.01 1.00 3588
ideo_c 0.01 1.00 3417
income_c:ideo_c 0.01 1.00 4883
party_c:ideo_c 0.01 1.00 2671
b[(Intercept) state:1] 0.01 1.00 2956
b[(Intercept) state:3] 0.00 1.00 5241
b[(Intercept) state:4] 0.01 1.00 3107
b[(Intercept) state:12] 0.00 1.00 3346
b[(Intercept) state:13] 0.00 1.00 4982
b[(Intercept) state:14] 0.00 1.00 4378
b[(Intercept) state:21] 0.01 1.00 2846
b[(Intercept) state:22] 0.01 1.00 2731
b[(Intercept) state:23] 0.00 1.00 4189
b[(Intercept) state:24] 0.00 1.00 5314
b[(Intercept) state:25] 0.00 1.00 4449
b[(Intercept) state:31] 0.00 1.00 4541
b[(Intercept) state:32] 0.00 1.00 4981
b[(Intercept) state:33] 0.01 1.00 2331
b[(Intercept) state:34] 0.00 1.00 4148
b[(Intercept) state:35] 0.00 1.00 4371
b[(Intercept) state:40] 0.00 1.00 4758
b[(Intercept) state:41] 0.00 1.00 4328
b[(Intercept) state:42] 0.00 1.00 3460
b[(Intercept) state:43] 0.00 1.00 4079
b[(Intercept) state:44] 0.00 1.00 3920
b[(Intercept) state:45] 0.00 1.00 4482
b[(Intercept) state:47] 0.00 1.00 4410
b[(Intercept) state:49] 0.00 1.00 4466
b[(Intercept) state:52] 0.00 1.00 4427
b[(Intercept) state:54] 0.00 1.00 3275
b[(Intercept) state:56] 0.00 1.00 4008
b[(Intercept) state:61] 0.00 1.00 4755
b[(Intercept) state:62] 0.01 1.00 3305
b[(Intercept) state:68] 0.00 1.00 4041
b[(Intercept) state:71] 0.00 1.00 4553
b[(Intercept) state:72] 0.00 1.00 4783
b[(Intercept) state:73] 0.00 1.00 4717
Sigma[state:(Intercept),(Intercept)] 0.00 1.00 1808
mean_PPD 0.00 1.00 5447
log-posterior 0.15 1.00 1310
For each parameter, mcse is Monte Carlo standard error, 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).The term (1 | state) means that each state is allowed to have its own baseline tendency to vote for Bush. These state-specific intercepts are not estimated completely separately. Instead, the model assumes that they come from a common state-level distribution. This is what creates partial pooling: each state estimate uses information from that state, but it is also stabilized by information from the overall distribution of states.
In practical terms, states with less information are pulled more strongly toward the national average, while states with more information are allowed to remain closer to their own observed pattern. This helps avoid overinterpreting noisy state-level estimates, especially when some states have small sample sizes.
The weakly informative prior normal(0, 2.5) is fairly wide on the log-odds scale. A coefficient of 2.5 corresponds to multiplying the odds by about 12, since exp(2.5) is approximately 12. This prior does not force coefficients to be small, but it discourages extremely large and implausible effects. In this way, the prior provides gentle regularization while still allowing the data to drive the results.
The prior on the intercept, normal(0, 5), is even wider. This gives the model flexibility in estimating the overall baseline log-odds of voting for Bush.
NoteConfidence interval vs credible interval
A 95% frequentist confidence interval describes the long-run behavior of an estimation procedure. If we repeatedly sampled data and built intervals in the same way, about 95% of those intervals would contain the true parameter.
A 95% Bayesian credible interval has a more direct interpretation. Given the model, prior, and data, 95% of the posterior distribution lies inside that interval.
These intervals may look similar when priors are weak and samples are large, but they come from different statistical frameworks and should not be interpreted in the same way.
9. Comparing glmer and stan_glmer
With weakly informative priors and enough data, the classical and Bayesian multilevel estimates are often very similar. This is what we see here. The glmer estimates and the stan_glmer posterior medians are close for most coefficients.
The main difference is in how we interpret the results. The glmer model gives a point estimate and standard error based on maximum likelihood estimation. The stan_glmer model gives a posterior distribution for each parameter. In the table below, stan_median is the middle of the posterior distribution, while stan_low and stan_high give the lower and upper bounds of the 95% credible interval.
glmer_coef <- fixef(fit_partial)
stan_draws <- as.matrix(fit_stan)
stan_fixed <- stan_draws[, names(glmer_coef), drop = FALSE]
compare_tbl <- tibble(
term = names(glmer_coef),
glmer_estimate = as.numeric(glmer_coef),
stan_median = apply(stan_fixed, 2, median),
stan_low = apply(stan_fixed, 2, quantile, probs = 0.025),
stan_high = apply(stan_fixed, 2, quantile, probs = 0.975)
)
kable(compare_tbl, digits = 2)| term | glmer_estimate | stan_median | stan_low | stan_high |
|---|---|---|---|---|
| (Intercept) | -0.54 | -0.53 | -0.91 | -0.18 |
| female | 0.21 | 0.20 | -0.26 | 0.65 |
| black | -2.60 | -2.53 | -3.88 | -1.41 |
| age_c | 0.07 | 0.07 | -0.39 | 0.54 |
| educ_c | 0.25 | 0.23 | -0.30 | 0.79 |
| income_c | -0.21 | -0.20 | -0.69 | 0.31 |
| party_c | 4.71 | 4.72 | 4.06 | 5.45 |
| ideo_c | 2.34 | 2.36 | 1.76 | 3.02 |
| income_c:ideo_c | 0.35 | 0.34 | -0.79 | 1.47 |
| party_c:ideo_c | -1.49 | -1.26 | -2.82 | 0.21 |
All coefficients are on the logit scale. Positive values mean higher log-odds of voting for Bush, while negative values mean lower log-odds of voting for Bush.
The comparison shows that the two models tell a very similar story. Party identification has the largest positive coefficient in both models. The glmer estimate is 4.71, and the Bayesian posterior median is 4.72. This means that voters who identify as more Republican are much more likely to vote for Bush. Ideology is also strongly positive: more conservative voters are more likely to vote for Bush.
The coefficient for Black respondents is strongly negative in both models. The glmer estimate is -2.60, and the Bayesian posterior median is -2.53. The Bayesian credible interval is also entirely below zero, from -3.88 to -1.41. This indicates that, after adjusting for the other predictors in the model, Black respondents were much less likely to vote for Bush.
For several other predictors, the Bayesian credible intervals include zero. For example, the intervals for female, age_c, educ_c, income_c, and income_c:ideo_c all cross zero. This means the data do not provide strong evidence that these predictors have clearly positive or clearly negative associations with voting for Bush, after accounting for the other variables in the model.
The interaction between party identification and ideology, party_c:ideo_c, is negative in both models. However, its Bayesian credible interval ranges from -2.82 to 0.21, which includes zero. This suggests that the interaction may be negative, but there is still considerable uncertainty about its size and direction.
Overall, the comparison is reassuring. The Bayesian model does not produce very different estimates from the classical multilevel model. Instead, it gives a similar set of results while adding a posterior interpretation. This allows us to describe uncertainty more directly and to use posterior simulations for predicted probabilities, state effects, and other quantities of interest.
10. Visualizing State-Level Uncertainty
So far, we have compared coefficients and state-level intercepts numerically. Another useful way to understand the multilevel model is to visualize how the predicted probability of voting for Bush changes across states.
Gelman and Hill’s (2007) Figure 14.2 in Analysis Using Regression and Multilevel/Hierarchical Models shows fitted logistic curves for several states. The idea is simple: we first summarize each respondent’s individual-level predictors into the fixed part of the linear predictor, and then we allow each state’s random intercept to shift the curve up or down.
In this plot, the x-axis represents the fixed part of the linear predictor:
\[ \eta_i = X_i\beta \]
This combines the individual-level information, such as gender, race, age, education, income, party identification, ideology, and the interaction terms.
Then the state-specific model adds a state random intercept:
\[ \eta_{ij} = X_{ij}\beta + a_j \]
where (a_j) is the estimated state effect for state (j). After applying the inverse-logit transformation, we get the predicted probability of voting for Bush:
\[ P(Y_{ij} = 1) = \frac{1}{1 + \exp[-(X_{ij}\beta + a_j)]} \]
The key point is that the state random intercept shifts the logistic curve. A state with a positive random intercept has a higher predicted probability of voting for Bush, holding individual characteristics constant. A state with a negative random intercept has a lower predicted probability.
First, we choose eight states with the largest sample sizes. Using states with larger samples makes the plots easier to read because there are more observed data points in each panel.
states_plot <- nes92 %>%
count(state) %>%
arrange(desc(n)) %>%
slice_head(n = 8) %>%
pull(state) %>%
as.character()
states_plot[1] "71" "13" "23" "40" "21" "49" "14" "44"10.1 Approximate posterior simulation with arm::sim()
The first version uses arm::sim() to simulate uncertainty from the fitted glmer() model. This is not full Bayesian MCMC, but it is useful for teaching because it gives a posterior-like simulation. The basic idea is:
- Fit the multilevel model using
glmer(). - Simulate plausible values of the model parameters.
- Use those simulated values to draw many possible prediction curves.
- Compare the curves across states.
The thin colored lines represent simulated curves. Each line is one plausible version of the state-specific relationship between the fixed linear predictor and the probability of voting for Bush. The black line shows the fitted curve using the estimated state random intercept. The jittered points show the observed vote outcomes, where 1 means voting for Bush and 0 means not voting for Bush.
beta <- fixef(fit_partial)
X <- model.matrix(fit_partial)
eta_fixed <- as.numeric(X %*% beta)
invlogit <- function(x) 1 / (1 + exp(-x))
x_grid <- seq(min(eta_fixed), max(eta_fixed), length.out = 200)
set.seed(216)
sims <- arm::sim(fit_partial, n.sims = 20)
state_names <- rownames(ranef(fit_partial)$state)
state_draw_mat <- sims@ranef$state[, , 1]
colnames(state_draw_mat) <- state_names
draw_curves_arm <- map_dfr(states_plot, function(s) {
re_draws <- state_draw_mat[, s]
expand.grid(x = x_grid, draw = seq_along(re_draws)) %>%
mutate(
state = s,
y = invlogit(x + re_draws[draw])
)
})
re_state <- ranef(fit_partial)$state[, 1]
names(re_state) <- state_names
median_curves <- expand.grid(x = x_grid, state = states_plot) %>%
mutate(y = invlogit(x + re_state[as.character(state)]))
plot_df <- nes92 %>%
mutate(eta_fixed = eta_fixed) %>%
filter(state %in% states_plot)
ggplot() +
geom_jitter(
data = plot_df,
aes(x = eta_fixed, y = vote_bush),
height = 0.03, alpha = 0.3, size = 0.6
) +
geom_line(
data = draw_curves_arm,
aes(x = x, y = y, group = draw),
color = "#5B8DB8", alpha = 0.25, linewidth = 0.3
) +
geom_line(
data = median_curves,
aes(x = x, y = y),
color = "black", linewidth = 0.9
) +
facet_wrap(~ state, nrow = 2) +
labs(
x = expression(eta[i] ~ "(fixed part of linear predictor)"),
y = "P(vote Bush)",
title = "State-specific curves from a varying-intercept logit model"
) +
theme_bw()
The plot should be read as follows. Moving from left to right on the x-axis means that the respondent’s individual characteristics make them more likely to vote for Bush. Moving upward on the y-axis means a higher predicted probability of voting for Bush.
The panels differ by state. If one state’s curve is generally higher than another state’s curve, this means that voters in that state are estimated to be more likely to vote for Bush, after holding individual-level characteristics constant. If the curve is lower, that state is estimated to be less pro-Bush, conditional on the same individual-level predictors.
The spread of the thin colored lines shows uncertainty. When the lines are close together, the model is relatively certain about that state’s curve. When the lines are more spread out, the model is less certain. States with less information usually have more uncertainty, and their estimated state effects are also more strongly pulled toward the overall average.
10.2 Full posterior curves from stan_glmer
Now we repeat the same idea using the Bayesian model from stan_glmer(). This time, the curves come from actual posterior draws. Each draw represents one plausible set of parameter values after combining the prior, the model, and the observed data.
This version has a more direct Bayesian interpretation. The thin curves show plausible state-specific prediction curves from the posterior distribution. The variation among these curves represents posterior uncertainty about the state-level effect.
post_mat <- as.matrix(fit_stan)
set.seed(216)
draw_idx <- sample(seq_len(nrow(post_mat)), 20)
find_state_col <- function(s) {
hits <- grep(paste0("state:", s, "\\]"), colnames(post_mat), value = TRUE)
if (length(hits) == 0) {
stop("Could not find posterior column for state ", s)
}
hits[1]
}
draw_curves_bayes <- map_dfr(states_plot, function(s) {
col <- find_state_col(s)
re_draws <- post_mat[draw_idx, col]
expand.grid(x = x_grid, draw = seq_along(re_draws)) %>%
mutate(
state = s,
y = invlogit(x + re_draws[draw])
)
})
ggplot() +
geom_jitter(
data = plot_df,
aes(x = eta_fixed, y = vote_bush),
height = 0.03, alpha = 0.3, size = 0.6
) +
geom_line(
data = draw_curves_bayes,
aes(x = x, y = y, group = draw),
color = "#8B3A8B", alpha = 0.25, linewidth = 0.3
) +
geom_line(
data = median_curves,
aes(x = x, y = y),
color = "black", linewidth = 0.9
) +
facet_wrap(~ state, nrow = 2) +
labs(
x = expression(eta[i] ~ "(fixed part of linear predictor)"),
y = "P(vote Bush)",
title = "State-specific curves from posterior draws"
) +
theme_bw()
The two plots should tell almost the same story. Both show how state-level intercepts shift the logistic curve upward or downward. Both also show that the model does not produce a single certain curve; instead, there is uncertainty around each state-specific prediction.
The difference is interpretive. The arm::sim() plot is based on a normal approximation to the fitted glmer() model. It is useful and often works well, but it is still an approximation. The stan_glmer() plot is based on draws from the posterior distribution, so it has a clearer Bayesian interpretation.
This tutorial is based on UCLA STATS C216: Applied Bayesian Social Statistics taught by Prof Mark Handcock.