lmerb()Chapter 01 fit a model with no group-level structure: every subject shrank toward the same grand mean. This chapter uses a dataset where the groups themselves differ in known ways, which is where the level-2 design earns its keep.
The data is big_word_club from the
bayesrules package (Johnson et
al. 2022) — a vocabulary intervention study in 47 schools.
Because bayesrules is only suggested and a 1000-draw
fit takes a few seconds, the model is fitted ahead of time by
data-raw/make_vignette_outputs.R and the results are loaded
from inst/extdata/. Every code chunk that fits is shown but
not evaluated.
data(big_word_club, package = "bayesrules")
dat <- big_word_club
dat$school_id <- factor(dat$school_id)
keep <- c("score_ppvt", "distracted_a1", "distracted_ppvt",
"private_school", "title1", "free_reduced_lunch", "school_id")
dat <- subset(
dat,
!is.na(score_ppvt) & !is.na(invalid_ppvt) & invalid_ppvt == 0L &
complete.cases(dat[, keep])
)
dat$school_id <- droplevels(dat$school_id)That leaves 551 students in 47 schools. The outcome
score_ppvt is a vocabulary score. Two student-level
covariates record whether the child was distracted — during the
vocabulary assessment (distracted_ppvt) and during an
earlier assessment (distracted_a1). Three school-level
covariates describe the school: private_school,
title1, and free_reduced_lunch (the percentage
of students on subsidized meals).
The distinction between the two kinds of covariate is the whole chapter. Student-level variables can have random slopes; school-level variables cannot, because they do not vary within a school.
form_simple <- score_ppvt ~ private_school + title1 +
distracted_ppvt + distracted_a1 +
(1 + distracted_ppvt + distracted_a1 || school_id)
model_setup(form_simple, data = dat)#> Call:
#> model_setup(formula = form_simple, data = dat)
#> --- Level 1 (Within-Group) ---
#> score_ppvt ~ (Intercept) + distracted_ppvt + distracted_a1
#> Observations : 551
#> Group-effect coefficients: 3
#> Group : school_id [47 levels]
#> Full-rank D_j: 34 of 47 groups
#> rank-deficient: 1, 4, 7, 8, 16, 21, 23, 27, 28, 32, ... (3 more)
#> Full-rank & estimable: 34 of 34 full-rank (34 of 47 total school_id)
#> --- Level 2 (Across-Group) ---
#> (Intercept) ~ 1 + private_school + title1
#> distracted_ppvt ~ 1
#> distracted_a1 ~ 1
#> --- Level 2 Rank ---
#> (Restricted to 34 estimable school_id)
#> (Intercept) groups=34 predictors=3 full-rank
#> distracted_ppvt groups=34 predictors=1 full-rank
#> distracted_a1 groups=34 predictors=1 full-rank
#> Rank-deficient flags:
#> (Intercept) FALSE
#> distracted_ppvt FALSE
#> distracted_a1 FALSE
#> popef.rank_ok: TRUE -- model rank looks OKLook at the Level 2 block. model_setup() has taken the
two school-level fixed effects and attached them to the
intercept’s level-2 row:
(Intercept) ~ 1 + private_school + title1
distracted_ppvt ~ 1
distracted_a1 ~ 1That is the right thing to do, and it is worth being explicit about
why. private_school does not vary within a school, so a
“fixed effect of private_school” is a statement about how
school intercepts differ. In the centered parameterization of
Chapter 02, it belongs in \(\mathcal{W}_j\), shifting each school’s
prior mean intercept. The two distraction slopes keep bare
~ 1 rows: their prior means are common across schools.
Also note the rank report. Only 34 of the 47 schools have a full-rank
within-group design — the rest lack the variation in
distracted_ppvt or distracted_a1 needed to
identify their own slopes.
They are not removed. All 47 schools remain in the fit and all 47 get coefficients; a rank-deficient school simply contributes no likelihood information in the deficient direction, so there it is drawn from its prior — shrunk all the way to \(\mathcal{W}_j\gamma\). That is the statistically sensible outcome, and it is why the model tolerates schools that a no-pooling analysis could not fit at all.
What the 34 are used for is the level-2 rank check, which runs on the estimable subset. That is the correct order: a school-level covariate that happens to be constant across the estimable schools would leave \(\gamma\) unidentified even though it varied in the raw data, and no amount of information from the other 13 would fix it.
Now the interesting specification. Does the effect of being distracted during the earlier assessment depend on the school’s poverty level?
form_full <- score_ppvt ~ private_school + title1 + free_reduced_lunch +
distracted_ppvt + distracted_a1 +
free_reduced_lunch:distracted_a1 +
(1 + distracted_ppvt + distracted_a1 || school_id)
design <- model_setup(form_full, data = dat)
design#> Call:
#> model_setup(formula = form_full, data = dat)
#> --- Level 1 (Within-Group) ---
#> score_ppvt ~ (Intercept) + distracted_ppvt + distracted_a1
#> Observations : 551
#> Group-effect coefficients: 3
#> Group : school_id [47 levels]
#> Full-rank D_j: 34 of 47 groups
#> rank-deficient: 1, 4, 7, 8, 16, 21, 23, 27, 28, 32, ... (3 more)
#> Full-rank & estimable: 34 of 34 full-rank (34 of 47 total school_id)
#> --- Level 2 (Across-Group) ---
#> (Intercept) ~ 1 + private_school + title1 + free_reduced_lunch
#> distracted_ppvt ~ 1
#> distracted_a1 ~ 1 + free_reduced_lunch
#> --- Level 2 Rank ---
#> (Restricted to 34 estimable school_id)
#> (Intercept) groups=34 predictors=4 full-rank
#> distracted_ppvt groups=34 predictors=1 full-rank
#> distracted_a1 groups=34 predictors=2 full-rank
#> Rank-deficient flags:
#> (Intercept) FALSE
#> distracted_ppvt FALSE
#> distracted_a1 FALSE
#> popef.rank_ok: TRUE -- model rank looks OKThe interaction has landed in the level-2 row of
distracted_a1:
(Intercept) ~ 1 + private_school + title1 + free_reduced_lunch
distracted_ppvt ~ 1
distracted_a1 ~ 1 + free_reduced_lunchThis is cross-level moderation. A school-level
variable is modifying an individual-level random slope. In
lme4 you would write the same interaction and get a single
fixed-effect coefficient; here you can see structurally what it does —
each school’s prior mean distracted_a1 slope now
sits on a line in free_reduced_lunch instead of at a common
value. Schools still shrink toward that line, just not toward the same
point.
Reading the Level 2 block before fitting is the fastest way to confirm you specified the model you intended. If an interaction shows up on the intercept row when you meant it on a slope row, the formula is wrong.
#> Call: Prior_Setup_GLMM()
#> --- Setup ---
#> Model family, prior-mean sources, and reference fit used for calibration.
#> family : gaussian (identity link)
#> pop.intercept_source : null_model
#> pop.effects_source : null_effects
#> dispformula : ~1
#> calibration_source : lme4 (fixef / tau^2_k / sd_tau)
#> --- Group dispersion (sigma^2) ---
#> Pooled Gamma prior on sigma^2 is stored in group.ing_prior.
#> (dGamma_list() requires a per-group dispformula.)
#> group.dispersion : 374.7 [pooled S_marg (A12 3.3.4)]
#> classical lmer sigma^2 : 357.2 (diagnostic)
#> --- Design check (from model_setup) ---
#> Identifiability flags copied from the design object (not re-checked here).
#> Full-rank groups (algebraic D_j): 34 of 47 school_id
#> Full-rank & estimable : 34 of 34 full-rank (34 of 47 total school_id)
#> --- Group-effect covariance (group.Sigma) ---
#> Diagonal Psi with entries tau^2_k for each group-effect coefficient.
#> (Intercept) distracted_ppvt distracted_a1
#> (Intercept) 199.2891 0.0000 0.0000
#> distracted_ppvt 0.0000 22.5391 0.0000
#> distracted_a1 0.0000 0.0000 42.0689
#> --- Population priors (pop.prior_list) ---
#> One prior list per group-effect coefficient (input to pfamily_list()).
#> [(Intercept)]
#> pop.pwt : 0.01
#> pop.dispersion.pwt : 0.01 [derived from pwt]
#> pop.dispersion.nprior : 0.4747 (= J * pwt_disp / (1 - pwt_disp))
#> mu:
#> (Intercept) private_school title1 free_reduced_lunch
#> 89.9099 0.0000 0.0000 0.0000
#> Sigma:
#> (Intercept) private_school title1 free_reduced_lunch
#> (Intercept) 4771.0901 -4724.7432 -411.9306 -56.4158
#> private_school -4724.7432 7337.0729 399.4243 56.0455
#> title1 -411.9306 399.4243 7177.9095 -77.5481
#> free_reduced_lunch -56.4158 56.0455 -77.5481 1.7389
#> dispersion : 168.7142 (tau^2_k; A12 3.3.4 marginal)
#> pop.dispersion.ref : 199.2891 (classical VarCorr / Psi)
#> ING tau^2 window : [105.3, 275.4] (0.01/0.99 limiting-posterior quantiles; upper/tau2 = 1.63)
#> ING shape, rate : 2.737, 293.1 (Gamma on 1/tau^2_k; for ptypes = "dIndependent_Normal_Gamma")
#> [distracted_ppvt]
#> pop.pwt : 0.01
#> pop.dispersion.pwt : 0.01 [derived from pwt]
#> pop.dispersion.nprior : 0.4747 (= J * pwt_disp / (1 - pwt_disp))
#> mu:
#> (Intercept)
#> 0
#> Sigma:
#> (Intercept)
#> (Intercept) 256.5316
#> dispersion : 4.3020 (tau^2_k; A12 3.3.4 marginal)
#> pop.dispersion.ref : 22.5391 (classical VarCorr / Psi)
#> ING tau^2 window : [2.686, 7.023] (0.01/0.99 limiting-posterior quantiles; upper/tau2 = 1.63)
#> ING shape, rate : 1.237, 1.021 (Gamma on 1/tau^2_k; for ptypes = "dIndependent_Normal_Gamma")
#> [distracted_a1]
#> pop.pwt : 0.01
#> pop.dispersion.pwt : 0.01 [derived from pwt]
#> pop.dispersion.nprior : 0.4747 (= J * pwt_disp / (1 - pwt_disp))
#> mu:
#> (Intercept) free_reduced_lunch
#> 0 0
#> Sigma:
#> (Intercept) free_reduced_lunch
#> (Intercept) 1646.1098 -20.0626
#> free_reduced_lunch -20.0626 0.3227
#> dispersion : 9.6366 (tau^2_k; A12 3.3.4 marginal)
#> pop.dispersion.ref : 42.0689 (classical VarCorr / Psi)
#> ING tau^2 window : [6.016, 15.73] (0.01/0.99 limiting-posterior quantiles; upper/tau2 = 1.63)
#> ING shape, rate : 1.737, 7.106 (Gamma on 1/tau^2_k; for ptypes = "dIndependent_Normal_Gamma")Three prior blocks, one per random coefficient, and their dimensions
follow the level-2 rows exactly: the intercept block has a 4-by-4
Sigma (intercept plus three school covariates),
distracted_ppvt has a scalar, and
distracted_a1 has 2-by-2. Every non-intercept
mu is zero, per the default
pop.effects_source = "null_effects" described in Chapter
03.
Each block also reports dispersion, which is \(\tau^2_k\) — how much individual school
coefficients scatter around their level-2 fitted value — alongside
pop.dispersion.ref, the classical VarCorr
estimate. The gap between them (168.7 versus 199.3 for the intercept) is
the marginal correction described in
vignette("Chapter-A12", package = "glmbayes").
set.seed(2024)
fit <- lmerb(
form_full,
data = dat,
pfamily_list = pfamily_list(ps),
dispersion_ranef = ps$group.dispersion,
n = 1000L,
progbar = FALSE
)
fit#> Call:
#> lmerb(formula = form_full, data = dat, pfamily_list = pfamily_list(ps),
#> dispersion_ranef = ps$group.dispersion, n = 1000L, progbar = FALSE)
#> Bayesian linear mixed model [1000 draws, exact iid]
#> Formula: score_ppvt ~ private_school + title1 + free_reduced_lunch + distracted_ppvt + distracted_a1 + free_reduced_lunch:distracted_a1 + (1 + distracted_ppvt + distracted_a1 || school_id)
#> Random effects (variance components fixed at lmer estimates):
#> Groups Name Std.Dev.
#> school_id (Intercept) 14.117
#> school_id.1 distracted_ppvt 4.748
#> school_id.2 distracted_a1 6.486
#> Residual 18.901
#> Number of obs: 551, groups: school_id, 47
#> --- Posterior means (ICM exact, under fixed variance components) ---
#> RE component parameter fixef.mode fixef.means draws SD
#> --------------- ------------------ ------------ ------------ ----------
#> (Intercept) (Intercept) 101.3426 101.3791 6.4806
#> (Intercept) private_school -3.5598 -3.4532 7.7906
#> (Intercept) title1 17.9754 17.8345 7.8408
#> (Intercept) free_reduced_lunch -0.3642 -0.3627 0.1189
#> distracted_ppvt (Intercept) -1.3948 -1.3151 1.3746
#> distracted_a1 (Intercept) 0.3282 0.3427 3.6902
#> distracted_a1 free_reduced_lunch -0.0865 -0.0880 0.0496exact iid again: the observation variance is fixed and
all three \(\tau^2_k\) are
dNormal, so the joint posterior is multivariate normal and
the 1000 draws are independent. m_convergence is 1 — no
Gibbs sweeps at all.
The population table has seven rows, one per level-2 parameter, and they group by which random coefficient they modify.
Intercept block. title1 is large and
positive (about 17.8 points) with a posterior SD near 7.8, so Title I
schools score higher after adjustment. private_school is
around -3.5 with an SD of 7.8 — the data does not distinguish it from
zero. free_reduced_lunch is small and negative, but its
scale is per percentage point, so a 40-point difference in
subsidized-meal rate moves the intercept by about 14 points.
Slope blocks. distracted_ppvt has a
population mean near -1.3: children distracted during the vocabulary
assessment score slightly lower, though the posterior SD of 1.4 makes
that weak evidence.
The moderation. The
distracted_a1 / free_reduced_lunch row is about -0.09 with
an SD of 0.05. The main distracted_a1 effect is near zero
(0.34, SD 3.7), so the story is not “distraction hurts” but “distraction
hurts more in higher-poverty schools” — and even that is modest
evidence. This is exactly the kind of claim that a single point estimate
plus a p-value tends to overstate and a posterior distribution does
not.
#> Call:
#> lmerb(formula = form_full, data = dat, pfamily_list = pfamily_list(ps),
#> dispersion_ranef = ps$group.dispersion, n = 1000L, progbar = FALSE)
#> Bayesian linear mixed model fit [1000 draws, exact iid]
#> Formula: score_ppvt ~ private_school + title1 + free_reduced_lunch + distracted_ppvt + distracted_a1 + free_reduced_lunch:distracted_a1 + (1 + distracted_ppvt + distracted_a1 || school_id)
#> Random effects (variance components fixed at lmer estimates):
#> Groups Name Std.Dev.
#> school_id (Intercept) 14.117
#> school_id.1 distracted_ppvt 4.748
#> school_id.2 distracted_a1 6.486
#> Residual 18.901
#> Number of obs: 551, groups: school_id, 47
#> === Block 2 dispersion (RE variance tau^2_k) ===
#> Prior and lmer reference:
#> Prior 1/E[1/tau2] E[tau2] sqrt(E[tau2]) disp_lower
#> (Intercept) dNormal 168.7142 168.7142 12.9890 NA
#> distracted_ppvt dNormal 4.3020 4.3020 2.0741 NA
#> distracted_a1 dNormal 9.6366 9.6366 3.1043 NA
#> Residual fixed 374.6978 374.6978 19.3571 NA
#> disp_upper lmer lmer SD
#> (Intercept) NA 199.2894 14.1170
#> distracted_ppvt NA 22.5390 4.7475
#> distracted_a1 NA 42.0690 6.4861
#> Residual NA 357.2463 18.9010
#> Overview:
#> Post.Mode Post.Mean Post.Sd MC Error Mean SD Pr(Prior_tail)
#> (Intercept) 168.714 168.714 0.000 0.000 12.989 NA
#> distracted_ppvt 4.302 4.302 0.000 0.000 2.074 NA
#> distracted_a1 9.637 9.637 0.000 0.000 3.104 NA
#> Residual 374.698 374.698 NA NA 19.357 NA
#> Cand/draw
#> (Intercept) 1
#> distracted_ppvt 1
#> distracted_a1 1
#> Residual NA
#> Distribution percentiles (tau^2):
#> 1.0% 2.5% 5.0% Median 95.0% 97.5% 99.0%
#> (Intercept) 168.714 168.714 168.714 168.714 168.714 168.714 168.714
#> distracted_ppvt 4.302 4.302 4.302 4.302 4.302 4.302 4.302
#> distracted_a1 9.637 9.637 9.637 9.637 9.637 9.637 9.637
#> Residual 374.698 374.698 374.698 374.698 374.698 374.698 374.698
#> SD credible interval (sqrt(tau^2) draws; lmer SD for reference):
#> 2.5% Median 97.5% lmer SD
#> (Intercept) 12.989 12.989 12.989 14.117
#> distracted_ppvt 2.074 2.074 2.074 4.747
#> distracted_a1 3.104 3.104 3.104 6.486
#> Residual 19.357 19.357 19.357 18.901
#> === Block 2: Level-2 fixed effects (hyperparameters) ===
#> Prior and lmer reference:
#> Prior Mean Prior.sd lmer lmer.se
#> (Intercept)::(Intercept) 89.90989 69.07308 101.75790 6.942
#> (Intercept)::private_school 0.00000 85.65672 -4.46840 8.609
#> (Intercept)::title1 0.00000 84.72254 18.01398 8.515
#> (Intercept)::free_reduced_lunch 0.00000 1.31866 -0.36767 0.133
#> distracted_ppvt::(Intercept) 0.00000 16.01660 -1.29086 1.610
#> distracted_a1::(Intercept) 0.00000 40.57228 1.05474 4.078
#> distracted_a1::free_reduced_lunch 0.00000 0.56806 -0.09111 0.057
#> Overview:
#> fixef.mode fixef.means Post.Sd MC.Error
#> (Intercept)::(Intercept) 101.34264 101.37912 6.48057 0.205
#> (Intercept)::private_school -3.55979 -3.45316 7.79061 0.246
#> (Intercept)::title1 17.97539 17.83450 7.84078 0.248
#> (Intercept)::free_reduced_lunch -0.36417 -0.36268 0.11892 0.004
#> distracted_ppvt::(Intercept) -1.39477 -1.31509 1.37460 0.043
#> distracted_a1::(Intercept) 0.32824 0.34273 3.69024 0.117
#> distracted_a1::free_reduced_lunch -0.08646 -0.08804 0.04960 0.002
#> Pr(Prior_tail)
#> (Intercept)::(Intercept) 0.037 *
#> (Intercept)::private_school 0.311
#> (Intercept)::title1 0.010 **
#> (Intercept)::free_reduced_lunch <2e-16 ***
#> distracted_ppvt::(Intercept) 0.173
#> distracted_a1::(Intercept) 0.466
#> distracted_a1::free_reduced_lunch 0.037 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Distribution percentiles:
#> 1.0% 2.5% 5.0% Median
#> (Intercept)::(Intercept) 87.214263 89.101783 90.631954 101.504704
#> (Intercept)::private_school -20.912827 -18.504666 -16.478518 -3.353277
#> (Intercept)::title1 0.060004 2.586376 5.425008 17.646881
#> (Intercept)::free_reduced_lunch -0.642895 -0.610479 -0.560587 -0.361539
#> distracted_ppvt::(Intercept) -4.408009 -4.017220 -3.638561 -1.254531
#> distracted_a1::(Intercept) -8.494159 -6.880089 -5.632714 0.284308
#> distracted_a1::free_reduced_lunch -0.204275 -0.184988 -0.172450 -0.085781
#> 95.0% 97.5% 99.0%
#> (Intercept)::(Intercept) 112.165188 114.087525 116.892
#> (Intercept)::private_school 9.313438 12.042670 15.392
#> (Intercept)::title1 31.094320 34.005236 35.731
#> (Intercept)::free_reduced_lunch -0.168963 -0.134516 -0.103
#> distracted_ppvt::(Intercept) 0.951829 1.180326 1.664
#> distracted_a1::(Intercept) 6.594520 7.874248 8.838
#> distracted_a1::free_reduced_lunch -0.005626 0.008956 0.023
#> === Block 1: Random effects (group-level) ===
#> Summary of posterior mode (ranef.mode) across groups:
#> Mean SD Min Q1 Median Q3 Max
#> (Intercept) 92.1210 14.0151 60.6991 82.5311 95.7448 103.3069 114.1528
#> distracted_ppvt -1.3953 0.5416 -2.8066 -1.5155 -1.3948 -1.1528 0.2402
#> distracted_a1 -4.4294 3.3158 -10.8201 -6.9601 -5.1185 -2.3465 1.3546
#> MCMC.mean
#> (Intercept) 92.143
#> distracted_ppvt -1.319
#> distracted_a1 -4.485
#> Per-group random effects: inspect fit$groupef.mode or fit$groupef,
#> or call summary(fit, groups = <level ids>) for selected groups.Three sections, and each answers a different question.
Block 2 dispersion compares each \(\tau^2_k\) against the lmer
reference. All three read dNormal in the Prior column and
show identical percentiles, because they are fixed rather than sampled.
Chapter 07 is what makes those rows interesting.
Block 2 level-2 fixed effects is the population
table again, with prior means and lmer reference values
alongside, plus Pr(Prior_tail). That last column is a
directional tail probability: how far the posterior has moved into the
prior’s tail. Values near 0.5 mean the data barely moved the
coefficient; small values mean it moved decisively.
Block 1 random effects summarizes the per-school
coefficients across schools rather than listing all 34. Use
summary(fit, groups = ...) for specific schools, or go to
the draws directly.
#> grpvar term grp condval condsd
#> 1 school_id (Intercept) 1 96.90886 NA
#> 2 school_id (Intercept) 2 82.51961 NA
#> 3 school_id (Intercept) 3 84.30990 NA
#> 4 school_id (Intercept) 4 103.22831 NA
#> 5 school_id (Intercept) 5 86.85338 NA
#> 6 school_id (Intercept) 6 70.46332 NA
#> 7 school_id (Intercept) 7 96.93794 NA
#> 8 school_id (Intercept) 8 87.66989 NA
#> 9 school_id (Intercept) 9 64.07830 NA
#> 10 school_id (Intercept) 10 60.76138 NA
#> 11 school_id (Intercept) 11 95.37196 NA
#> 12 school_id (Intercept) 12 89.51163 NARemember from Chapter 02 that these are full coefficients, not deviations: school 1’s intercept of about 97 is its expected score, not an offset. The long format has one row per school and term, which makes it convenient for plotting or joining back to school-level metadata.
Chapter 01 (sleepstudy) |
Chapter 04 (big_word_club) |
|
|---|---|---|
| Level-2 rows | all ~ 1 |
intercept and one slope carry covariates |
| Shrinkage target | one grand mean per coefficient | a fitted line per coefficient |
| Groups retained | 18 of 18 | 34 of 47 |
| Cross-level terms | none | free_reduced_lunch:distracted_a1 |
The sampler did not change at all. Both fits are exact iid draws from
a multivariate normal posterior. What changed is the model, and
model_setup() is where you see it.
lme4 and glmmTMB.?lmerb, ?model_setup