Binary outcomes with group structure. The data is
book_banning from bayesrules (Johnson et al. 2022): challenges to books in US
libraries, with an indicator for whether the book was removed, grouped
by state.
The mechanics are the same as Chapter 09 — glmerb() with
family = binomial() — so this chapter concentrates on what
is different about binary data: separation, estimability, and why the
sweep count is larger.
data(book_banning, package = "bayesrules")
bb <- book_banning[, c("state", "removed", "violent")]
bb <- bb[complete.cases(bb), ]
bb$removed_i <- as.integer(bb$removed == 1L | bb$removed == "1")
bb$violent_i <- as.integer(
bb$violent == TRUE | bb$violent == 1L | bb$violent == "TRUE"
)
bb$state <- factor(bb$state)
## Twelve largest states, to keep the stored fit small.
keep <- names(sort(table(bb$state), decreasing = TRUE))[seq_len(12L)]
bb <- droplevels(subset(bb, state %in% keep))
form <- removed_i ~ violent_i + (1 + violent_i || state)
model_setup(form, data = bb, family = binomial())#> Call:
#> model_setup(formula = form, data = bb, family = binomial())
#> --- Level 1 (Within-Group) ---
#> removed_i ~ (Intercept) + violent_i
#> Observations : 629
#> Group-effect coefficients: 2
#> Group : state [12 levels]
#> Full-rank D_j: 12 of 12 groups
#> Full-rank & estimable: 11 of 12 full-rank (11 of 12 total state)
#> full-rank but not estimable: MI
#> --- Level 2 (Across-Group) ---
#> (Intercept) ~ 1
#> violent_i ~ 1
#> --- Level 2 Rank ---
#> (Restricted to 11 estimable state)
#> (Intercept) groups=11 predictors=1 full-rank
#> violent_i groups=11 predictors=1 full-rank
#> Rank-deficient flags:
#> (Intercept) FALSE
#> violent_i FALSE
#> popef.rank_ok: TRUE -- model rank looks OKNote the explicit conversion to integers. Logical and factor
predictors get expanded by model.matrix() in ways that can
produce correlated random-effect terms, which this package rejects;
converting to 0/1 integers keeps the || structure intact.
If you hit Correlated random effects are not supported,
check for logical or factor variables inside the random-effects term
first.
The rank report shows something that did not appear in earlier chapters:
Full-rank D_j: 12 of 12 groups
Full-rank & estimable: 11 of 12 full-rank (11 of 12 total state)
full-rank but not estimable: MIMichigan’s design matrix has full column rank — it contains both
violent and non-violent challenges — but the model is still not
estimable there. That is separation: within Michigan,
the outcome is perfectly predicted by violent_i, so the
maximum-likelihood coefficient runs off to infinity. The likelihood is
log-concave but has no interior maximum in that direction.
Separation is a binary-data problem specifically; it cannot happen
with Gaussian or (in this form) Poisson responses.
model_setup() runs a separate estimability check for
exactly this reason, and it is the reason the check is reported on its
own line rather than folded into the rank test.
Michigan is not thrown away. It keeps its coefficients, drawn from
the prior in the direction the likelihood cannot pin down — which is
exactly what a proper prior is for. What it loses is standing in the
level-2 checks: an infinite-in-the-limit estimate carries no
usable information about the population parameters, so
model_setup() restricts the across-group rank condition to
the 11 estimable states, and Prior_Setup_GLMM() calibrates
from those. Details in
vignette("Chapter-B08", package = "lmebayesCore").
ps <- Prior_Setup_GLMM(form, data = bb, family = binomial(), pop.pwt = 0.01)
set.seed(2024)
fit <- glmerb(
form,
data = bb,
family = binomial(),
pfamily_list = pfamily_list(ps),
n = 200L,
progbar = FALSE
)
fit#> Call:
#> glmerb(formula = form, data = bb, family = binomial(), pfamily_list = pfamily_list(ps),
#> n = 200L, progbar = FALSE)
#> Bayesian generalized linear mixed model [binomial; 200 draws, two-block Gibbs]
#> Formula: removed_i ~ violent_i + (1 + violent_i || state)
#> Random effects (variance components fixed at glmer estimates):
#> Groups Name Std.Dev.
#> state (Intercept) 0.8630
#> state.1 violent_i 0.8198
#> Number of obs: 629, groups: state, 12
#> --- Block 2 hyperparameters (ICM posterior mode; MCMC means when simulated) ---
#> RE component parameter ICM mode popef.means draws SD
#> ------------ ----------- ------------ ------------ ----------
#> (Intercept) (Intercept) -1.1770 -1.2250 0.3586
#> violent_i (Intercept) 0.3342 0.2685 0.3985Coefficients are log-odds. Two useful transformations:
pm <- unlist(ch10$popef_means)
round(rbind(log_odds = pm, odds_ratio = exp(pm)), 3)
#> (Intercept).(Intercept) violent_i.(Intercept)
#> log_odds -1.225 0.268
#> odds_ratio 0.294 1.308
plogis(pm[[1]]) # P(removed), non-violent
#> [1] 0.2270659
plogis(pm[[1]] + pm[[2]]) # P(removed), violent
#> [1] 0.2775859In an average state, a challenged non-violent book has roughly a 23
percent chance of removal, rising modestly for books challenged as
violent. The posterior SD on the violent_i effect is larger
than the effect itself, so this is weak evidence of any population-level
difference.
Note that plogis() of the mean log-odds is not
the mean probability — the transformation is nonlinear. For a proper
posterior on the probability scale, transform the draws and then
summarize:
#> grpvar term grp condval condsd
#> 1 state (Intercept) CA -0.8438486 NA
#> 2 state (Intercept) CO -1.9240385 NA
#> 3 state (Intercept) FL -0.9772028 NA
#> 4 state (Intercept) IL -0.9855758 NA
#> 5 state (Intercept) IN -1.0040763 NA
#> 6 state (Intercept) MI -1.3526249 NA
#> 7 state (Intercept) NY -0.4239438 NA
#> 8 state (Intercept) OH -0.8093770 NA
#> 9 state (Intercept) OK -0.7953365 NA
#> 10 state (Intercept) OR -2.8705434 NA
#> 11 state (Intercept) PA -2.5431215 NA
#> 12 state (Intercept) VA -0.0915660 NAThese are log-odds intercepts per state, and they are full coefficients rather than deviations (Chapter 02). Oregon at about -2.87 and Virginia at about -0.09 correspond to removal probabilities of roughly 5 percent and 48 percent — a large practical difference that a state-level fixed effect would have estimated far more noisily, and that complete pooling would have erased entirely.
8 sweeps, against 3 for the Poisson model of Chapter 09 and 1 for the Gaussian fits.
The reason is information per observation. A binary response carries at most about 0.25 units of Fisher information per observation, and much less when the probability is near 0 or 1 — which it is here, since most challenges do not result in removal. With only 12 states and weak per-observation information, the group coefficients and the population coefficients are tightly coupled, the rate \(\lambda^*\) is close to one, and more sweeps are needed.
This is the same relationship as in Chapter 02: models that shrink hard mix slowly. Binary outcomes with few groups and rare events sit at the difficult end of that spectrum, and the sweep count reports it honestly rather than hiding it.
#> Call:
#> glmerb(formula = form, data = bb, family = binomial(), pfamily_list = pfamily_list(ps),
#> n = 200L, progbar = FALSE)
#> Bayesian linear mixed model fit [200 draws, two-block Gibbs]
#> Formula: removed_i ~ violent_i + (1 + violent_i || state)
#> Random effects (variance components fixed at glmer estimates):
#> Groups Name Std.Dev.
#> state (Intercept) 0.8630
#> state.1 violent_i 0.8198
#> Number of obs: 629, groups: state, 12
#> === Block 2 dispersion (RE variance tau^2_k) ===
#> Prior and glmer reference:
#> Prior 1/E[1/tau2] E[tau2] sqrt(E[tau2]) disp_lower disp_upper
#> (Intercept) dNormal 0.7448 0.7448 0.8630 NA NA
#> violent_i dNormal 0.6721 0.6721 0.8198 NA NA
#> glmer glmer SD
#> (Intercept) 0.7448 0.8630
#> violent_i 0.6721 0.8198
#> Overview:
#> Post.Mode Post.Mean Post.Sd MC Error Mean SD Pr(Prior_tail)
#> (Intercept) 0.7448 0.7448 0.0000 0.0000 0.8630 NA
#> violent_i 0.6721 0.6721 0.0000 0.0000 0.8198 NA
#> Cand/draw
#> (Intercept) 1
#> violent_i 1
#> Distribution percentiles (tau^2):
#> 1.0% 2.5% 5.0% Median 95.0% 97.5% 99.0%
#> (Intercept) 0.7448 0.7448 0.7448 0.7448 0.7448 0.7448 0.745
#> violent_i 0.6721 0.6721 0.6721 0.6721 0.6721 0.6721 0.672
#> SD credible interval (sqrt(tau^2) draws; glmer SD for reference):
#> 2.5% Median 97.5% glmer SD
#> (Intercept) 0.8630 0.8630 0.8630 0.863
#> violent_i 0.8198 0.8198 0.8198 0.820
#> === Block 2: Level-2 fixed effects (hyperparameters) ===
#> Prior and glmer reference:
#> Prior Mean Prior.sd glmer glmer.se
#> (Intercept)::(Intercept) -1.1702 2.7541 -1.2075 0.277
#> violent_i::(Intercept) 0.0000 4.1137 0.2859 0.413
#> Overview:
#> fixef.mode fixef.means Post.Sd MC.Error Pr(Prior_tail)
#> (Intercept)::(Intercept) -1.1770 -1.2250 0.3586 0.025 0.43
#> violent_i::(Intercept) 0.3342 0.2685 0.3985 0.028 0.22
#> Distribution percentiles:
#> 1.0% 2.5% 5.0% Median 95.0% 97.5% 99.0%
#> (Intercept)::(Intercept) -1.9542 -1.8585 -1.7845 -1.2609 -0.6676 -0.5536 -0.44
#> violent_i::(Intercept) -0.6092 -0.5794 -0.4291 0.2895 0.8636 0.9277 1.28
#> === Block 1: Random effects (group-level) ===
#> Summary of posterior mode (ranef.mode) across groups:
#> Mean SD Min Q1 Median Q3 Max
#> (Intercept) -1.17703 0.80772 -2.76751 -1.47585 -0.92871 -0.77131 -0.07223
#> violent_i 0.33533 0.50096 -0.23663 0.11177 0.20765 0.44171 1.73505
#> MCMC.mean
#> (Intercept) -1.218
#> violent_i 0.288
#> Per-group random effects: inspect fit$groupef.mode or fit$groupef,
#> or call summary(fit, groups = <level ids>) for selected groups.Both random-effect standard deviations are near 0.85 on the log-odds scale, which is substantial: a one-standard-deviation state differs from average by a factor of about 2.4 in odds. As in Chapter 09, the between-group variation exceeds the population-level effect, which is the signature of a model that needed random slopes.
m_convergence
climbs into the tens, the cause is usually too few groups or too rare an
outcome, not a sampler problem.lambda* > 0.95 warning.
The chain is still valid, but the model is close to the boundary where
the group and population levels cannot be separated.vignette("Chapter-B08", package = "lmebayesCore") — the
rank and estimability conditions.?glmerb, ?model_setup