Hierarchical binary logit via bayesm (R) — rich-state / poor-state paradox¶
Companion cross-check to hier_logit_python.ipynb¶
References: Rossi, Allenby & McCulloch (2005), Bayesian Statistics and Marketing; Gelman, Shor, Bafumi & Park (2007).
| Section | Content |
|---|---|
| Model | Hierarchical binary logit; rhierBinLogit (random-walk Metropolis); references |
| Section 1 | Gelman 2004 exit-poll data — recover the paradox; compare to the Python from-scratch & PyMC fits |
Uses the same cleaned data as the Python notebook (gelman_2004.csv, gelman_2004_states.csv). bayesm fits the multilevel binary logit by random-walk Metropolis within a conjugate Gibbs hierarchy — a direct match to the from-scratch Pólya–Gamma model, with a different logit-likelihood sampler.
Model — hierarchical binary logit (rhierBinLogit)¶
The same multilevel model as the Python notebook:
$$\Pr(y_i=1)=\text{logit}^{-1}(\alpha_j + \beta_j\,\text{income}_i),\qquad (\alpha_j,\beta_j)' \sim N(\Delta' z_j,\ V_\beta),\quad z_j=(1,u_j).$$
bayesm::rhierBinLogit fits this with a per-state random-walk Metropolis step for each $(\alpha_j,\beta_j)$ and conjugate Gibbs updates for the hierarchy ($\Delta$ matrix-normal, $V_\beta$ inverse-Wishart) — the same structure as our from-scratch sampler, but using Metropolis (not Pólya–Gamma augmentation) for the logit likelihood.
Data format: lgtdata = list (one per state) of list(y = n_j×1 matrix of 0/1, X = [1, income]); Z = $[1, u]$ ($J\times2$). The output Deltadraw has columns in vec(Delta) order: [D00, D10, D01, D11] = [baseline intercept, state-income→intercept, mean income slope, state-income→slope].
References: Rossi, Allenby & McCulloch (2005); Gelman, Shor, Bafumi & Park (2007).
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
suppressMessages(library(bayesm))
d <- read.csv("C:/Users/user/project/gelman_2004.csv")
st <- read.csv("C:/Users/user/project/gelman_2004_states.csv")
J <- nrow(st)
lgt <- vector("list", J)
for (j in 1:J) {
ix <- which(d$state == (j - 1)) # state 0-indexed in CSV
lgt[[j]] <- list(y = matrix(d$y[ix], ncol = 1), # y MUST be a 1-col matrix here
X = cbind(1, d$income[ix])) # [intercept, income]
}
Z <- cbind(1, st$u) # J x 2 group-level design [1, u]
cat(sprintf("n=%d J=%d states Republican (Bush) share=%.3f\n", nrow(d), J, mean(d$y)))
n=67381 J=50 states Republican (Bush) share=0.497
set.seed(1)
out <- rhierBinLogit(Data = list(lgtdata = lgt, Z = Z), Mcmc = list(R = 10000, keep = 5))
Dd <- out$Deltadraw; Dd <- Dd[(nrow(Dd) %/% 4 + 1):nrow(Dd), ]
mn <- colMeans(Dd); lo <- apply(Dd, 2, quantile, .025); hi <- apply(Dd, 2, quantile, .975)
nm <- c("baseline alpha (D00)", "state-income->intercept (D10)",
"mean income slope (D01)", "state-income->slope (D11)")
py <- c(0.107, -0.310, 0.210, -0.042) # Python from-scratch (notebook §2)
cat(sprintf("\n%-32s %10s %10s %s\n", "parameter", "Python", "bayesm", "95% CI"))
cat(strrep("-", 72), "\n")
for (k in 1:4)
cat(sprintf("%-32s %10.3f %10.3f [%6.3f, %6.3f]\n", nm[k], py[k], mn[k], lo[k], hi[k]))
Warning message in rhierBinLogit(Data = list(lgtdata = lgt, Z = Z), Mcmc = list(R = 10000, : "'rhierBinLogit' is depricated Use 'rhierMnlRwMixture' instead"
Table of Y values pooled over all units
ypooled
0 1
33897 33484
Attempting MCMC Inference for Hierarchical Binary Logit:
2 variables in X
2 variables in Z
for 50 cross-sectional units
Prior Parms:
nu = 5
V
[,1] [,2]
[1,] 5 0
[2,] 0 5
Deltabar
[,1] [,2]
[1,] 0 0
[2,] 0 0
ADelta
[,1] [,2]
[1,] 0.01 0.00
[2,] 0.00 0.01
MCMC Parms:
sbeta= 0.2 R= 10000 keep= 5
MCMC Iteration (est time to end - min)
100 ( 2.9 )
200 ( 2.9 )
300 ( 2.9 )
400 ( 2.9 )
500 ( 2.9 )
600 ( 2.9 )
700 ( 2.8 )
800 ( 2.8 )
900 ( 2.8 )
1000 ( 2.7 )
1100 ( 2.7 )
1200 ( 2.7 )
1300 ( 2.6 )
1400 ( 2.6 )
1500 ( 2.6 )
1600 ( 2.6 )
1700 ( 2.5 )
1800 ( 2.5 )
1900 ( 2.5 )
2000 ( 2.4 )
2100 ( 2.4 )
2200 ( 2.4 )
2300 ( 2.4 )
2400 ( 2.3 )
2500 ( 2.3 )
2600 ( 2.3 )
2700 ( 2.2 )
2800 ( 2.2 )
2900 ( 2.2 )
3000 ( 2.2 )
3100 ( 2.1 )
3200 ( 2.1 )
3300 ( 2.1 )
3400 ( 2 )
3500 ( 2 )
3600 ( 2 )
3700 ( 1.9 )
3800 ( 1.9 )
3900 ( 1.9 )
4000 ( 1.8 )
4100 ( 1.8 )
4200 ( 1.8 )
4300 ( 1.8 )
4400 ( 1.7 )
4500 ( 1.7 )
4600 ( 1.7 )
4700 ( 1.6 )
4800 ( 1.6 )
4900 ( 1.6 )
5000 ( 1.6 )
5100 ( 1.6 )
5200 ( 1.5 )
5300 ( 1.5 )
5400 ( 1.5 )
5500 ( 1.5 )
5600 ( 1.4 )
5700 ( 1.4 )
5800 ( 1.4 )
5900 ( 1.4 )
6000 ( 1.3 )
6100 ( 1.3 )
6200 ( 1.3 )
6300 ( 1.3 )
6400 ( 1.2 )
6500 ( 1.2 )
6600 ( 1.2 )
6700 ( 1.1 )
6800 ( 1.1 )
6900 ( 1.1 )
7000 ( 1 )
7100 ( 1 )
7200 ( 1 )
7300 ( 1 )
7400 ( 0.9 )
7500 ( 0.9 )
7600 ( 0.9 )
7700 ( 0.8 )
7800 ( 0.8 )
7900 ( 0.8 )
8000 ( 0.7 )
8100 ( 0.7 )
8200 ( 0.7 )
8300 ( 0.6 )
8400 ( 0.6 )
8500 ( 0.5 )
8600 ( 0.5 )
8700 ( 0.5 )
8800 ( 0.4 )
8900 ( 0.4 )
9000 ( 0.4 )
9100 ( 0.3 )
9200 ( 0.3 )
9300 ( 0.3 )
9400 ( 0.2 )
9500 ( 0.2 )
9600 ( 0.2 )
9700 ( 0.1 )
9800 ( 0.1 )
9900 ( 0 )
10000 ( 0 )
Total Time Elapsed: 3.86
parameter Python bayesm 95% CI
------------------------------------------------------------------------
baseline alpha (D00) 0.107 0.106 [-0.024, 0.248]
state-income->intercept (D10) -0.310 -0.310 [-0.448, -0.176]
mean income slope (D01) 0.210 0.211 [ 0.117, 0.302]
state-income->slope (D11) -0.042 -0.043 [-0.134, 0.055]
# Per-state random effects (posterior means) + the paradox, state-labeled
bd <- out$betadraw; nd <- dim(bd)[3]; ki <- (nd %/% 4 + 1):nd
B <- apply(bd[, , ki, drop = FALSE], c(1, 2), mean) # J x 2 = (alpha_j, beta_j)
u <- st$u; ab <- st$abbr; repsh <- st$rep_share
options(repr.plot.width = 14, repr.plot.height = 6)
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
# Between states: Republican share vs state income, state-labeled (the famous scatter)
plot(u, repsh, pch = 19, col = ifelse(repsh > 0.5, "firebrick", "navy"), cex = 1.4,
xlab = "state mean income (standardized)", ylab = "Republican (Bush 2004) vote share",
main = "Richer STATES vote more Democratic")
abline(lm(repsh ~ u), lty = 2, lwd = 2); abline(h = 0.5, col = "gray", lty = 3)
key <- c("CT","NJ","MD","MA","NY","NH","CA","MS","ID","WY","WV","AR","OK","UT","TX","OH","FL")
sel <- ab %in% key; text(u[sel], repsh[sel], ab[sel], pos = 3, cex = 0.75, font = 2)
# Within states: P(Republican) vs income, CT/MA vs MS/WV highlighted
xi <- seq(-2, 2, length = 50)
plot(NA, xlim = c(-2, 2), ylim = c(0, 1), xlab = "individual income (-2 = poor .. +2 = rich)",
ylab = "P(vote Republican)", main = "Within EVERY state, richer -> more Republican")
for (j in 1:nrow(B)) lines(xi, plogis(B[j, 1] + B[j, 2] * xi), col = "gray80", lwd = 0.7)
for (pr in list(c("MS","firebrick"), c("WV","salmon"), c("CT","navy"), c("MA","royalblue"))) {
j <- which(ab == pr[1]); lines(xi, plogis(B[j, 1] + B[j, 2] * xi), col = pr[2], lwd = 2.6)
}
abline(h = 0.5, col = "gray", lty = 3)
legend("topleft", c("MS","WV","CT","MA"), col = c("firebrick","salmon","navy","royalblue"),
lwd = 2.6, bty = "n")
Results — bayesm rhierBinLogit vs. the Python fits¶
| parameter | Python from-scratch | bayesm (R) |
|---|---|---|
| baseline α (Δ₀₀) | 0.107 | 0.106 |
| mean income slope (Δ₀₁) | 0.210 | 0.211 [0.117, 0.302] |
| state-income → intercept (Δ₁₀) | −0.310 | −0.310 [−0.448, −0.176] |
| state-income → slope (Δ₁₁) | −0.042 | −0.043 |
Near-identical means. bayesm samples the logit likelihood by random-walk Metropolis (not Pólya–Gamma) within the same conjugate hierarchy, yet recovers the same group-level regression to 0.001 — and the same paradox: $\Delta_{01} > 0$ (within states, richer → Republican) and $\Delta_{10} < 0$ (between states, richer state → Democratic).
Reading the figure. Left — the famous between-state scatter: each point is a state, plotted at its mean income against its Bush vote share, and the fitted line slopes down. The rich states (CT, NJ, MD, MA) sit in the lower right — Democratic; the poor ones (MS, WV, OK, AR, ID) in the upper left — Republican. Right — one logit curve per state from the fitted random effects, with MS/WV (red) and CT/MA (blue) picked out. Every one of the 50 curves slopes up: inside any given state, richer voters are more Republican. But the rich states' curves sit lower — a Connecticut voter is less Republican than a Mississippi voter at the same income. Both panels are true at once, and that is the whole paradox: the sign of the income effect depends on whether you look across states or within them.
Why are bayesm's intervals wider? The credible intervals here are roughly 1.2–1.8× wider than the Python ones (Δ₀₁: [0.117, 0.302] vs [0.160, 0.265]) even though the means coincide. This is a prior difference, not a disagreement. The two fits place different inverse-Wishart priors on the random-effect covariance $V_\beta$: the from-scratch sampler uses $\mathrm{IW}(\nu = 4,\ S = I)$, whose prior mean is $I$, while bayesm's default is $\mathrm{IW}(\nu = 5,\ V = 5I)$, whose prior mean is $2.5\,I$. bayesm therefore expects larger state-to-state dispersion; the posterior $V_\beta$ comes out bigger, and since the conditional precision of $\Delta$ scales with $V_\beta^{-1}$, its intervals widen to match. With only $J = 50$ states the $V_\beta$ prior still has real influence — the same lesson as the $J$-vs-n_per sweep in the Python notebook. The location of the posterior is driven by the data (hence the agreement to 0.001); its width still remembers the prior.
On the deprecation warning. bayesm prints 'rhierBinLogit' is depricated — Use 'rhierMnlRwMixture' instead [sic]. rhierBinLogit is the purpose-built binary-logit routine and is kept here because it maps one-to-one onto the model written above; rhierMnlRwMixture is the more general multinomial routine, which fits a binary logit as the 2-choice special case with a mixture-of-normals (rather than a single normal) prior on the random effects. For this cross-check the older, simpler routine is the cleaner comparison. It still runs correctly — the warning is about maintenance, not validity.
Three-way cross-validation complete: the from-scratch hierarchical Pólya–Gamma Gibbs (Python), PyMC NUTS, and bayesm RW-Metropolis (R) all agree, on synthetic data and on the real Gelman–Shor–Bafumi–Park 2004 exit poll. Three different samplers, two languages, one conclusion.