This chapter fits one model end to end. The data is
lme4::sleepstudy: reaction times for 18 subjects over ten
days of restricted sleep. Everyone slows down; the interesting question
is by how much, and how much that varies from person to person.
data(sleepstudy, package = "lme4")
dat <- sleepstudy
dat$Days_c <- dat$Days - mean(dat$Days)
str(dat)
#> 'data.frame': 180 obs. of 4 variables:
#> $ Reaction: num 250 259 251 321 357 ...
#> $ Days : num 0 1 2 3 4 5 6 7 8 9 ...
#> $ Subject : Factor w/ 18 levels "308","309","310",..: 1 1 1 1 1 1 1 1 1 1 ...
#> $ Days_c : num -4.5 -3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 4.5 ...Days is centered. Centering is not required, but it
makes the intercept mean reaction time at the average day
rather than at day zero, which gives the random intercept and random
slope a cleaner interpretation and reduces their correlation — useful
here, because the package assumes they are uncorrelated.
Read this as: reaction time depends on days of deprivation, and both the baseline level and the per-day effect vary by subject.
The double bar || is required. lme4 would
let you write (1 + Days_c | Subject), which additionally
estimates a correlation between each subject’s intercept and slope.
lmebayes does not support that; it needs a diagonal
random-effect covariance so the Block 2 update factorizes across
coefficients. Chapter 02 explains why, and Chapter 13 discusses what you
give up.
Before fitting anything, model_setup() shows how the
formula is turned into matrices.
design <- model_setup(form, data = dat)
design
#> Call:
#> model_setup(formula = form, data = dat)
#>
#> --- Level 1 (Within-Group) ---
#> Reaction ~ (Intercept) + Days_c
#>
#> Observations : 180
#> Group-effect coefficients: 2
#> Group : Subject [18 levels]
#> Full-rank D_j: 18 of 18 groups
#> Full-rank & estimable: 18 of 18 full-rank (18 of 18 total Subject)
#>
#> --- Level 2 (Across-Group) ---
#> (Intercept) ~ 1
#> Days_c ~ 1
#>
#> --- Level 2 Rank ---
#> (Restricted to 18 estimable Subject)
#>
#> (Intercept) groups=18 predictors=1 full-rank
#> Days_c groups=18 predictors=1 full-rank
#>
#> Rank-deficient flags:
#> (Intercept) FALSE
#> Days_c FALSE
#>
#> popef.rank_ok: TRUE -- model rank looks OKThree things to check here.
Level 1 is the within-group model. Each subject gets
their own intercept and their own Days_c slope, so there
are two group-effect coefficients per subject.
Level 2 is the across-group model. Both lines read
~ 1, meaning neither coefficient has group-level
predictors: subjects are exchangeable, and all we estimate at the
population level is the average intercept and the average slope. Chapter
04 shows what happens when you add school-level covariates.
Rank matters.
Full-rank D_j: 18 of 18 groups says every subject has
enough distinct Days_c values to identify their own slope.
Subjects who fail this check are flagged rather than removed: they stay
in the fit, but in the deficient direction they carry no likelihood
information and are effectively drawn from their prior. What they
are excluded from is the across-group rank check and the prior
calibration. Chapter 04 shows a dataset where this happens, and
vignette("Chapter-B08", package = "lmebayesCore") gives the
conditions in full.
Prior_Setup_GLMM() fits a reference lmer()
model and uses it to calibrate priors on a scale that matches the
data.
pop.pwt = 0.01 sets the prior weight: the
population-level prior should carry about as much information as one
percent of the observations. Small values mean weak priors and results
close to lmer(); larger values pull the population effects
toward the prior mean. Chapter 03 unpacks every field of this
object.
pf <- pfamily_list(ps)
pf
#>
#> [[(Intercept)]]
#>
#> Call: dNormal(mu = mu_k, Sigma = Sig_k, dispersion = d_k)
#>
#> Prior Family: dNormal
#>
#> Prior List:
#>
#> $mu
#> [,1]
#> (Intercept) 298.5079
#>
#> $Sigma
#> (Intercept)
#> (Intercept) 8108.233
#>
#> $dispersion
#> [1] 1346.146
#>
#> $ddef
#> [1] FALSE
#>
#> attr(,"Prior Type")
#> [1] "dNormal"
#>
#> [[Days_c]]
#>
#> Call: dNormal(mu = mu_k, Sigma = Sig_k, dispersion = d_k)
#>
#> Prior Family: dNormal
#>
#> Prior List:
#>
#> $mu
#> [,1]
#> (Intercept) 0
#>
#> $Sigma
#> (Intercept)
#> (Intercept) 236.5569
#>
#> $dispersion
#> [1] 29.37223
#>
#> $ddef
#> [1] FALSE
#>
#> attr(,"Prior Type")
#> [1] "dNormal"pfamily_list() converts the calibration into one prior
object per random-effect coefficient — here, one for the intercept and
one for Days_c. Both are dNormal, meaning the
corresponding variance component tau^2_k is held fixed at
the lmer estimate. Chapter 07 shows how to sample it
instead.
set.seed(1)
fit <- lmerb(
form,
data = dat,
pfamily_list = pf,
dispersion_ranef = ps$group.dispersion,
n = 500L,
progbar = FALSE
)n = 500 is the number of posterior draws.
dispersion_ranef supplies the observation variance
sigma^2; here it is the fixed lmer
estimate.
fit
#> Call:
#> lmerb(formula = form, data = dat, pfamily_list = pf, dispersion_ranef = ps$group.dispersion,
#> n = 500L, progbar = FALSE)
#>
#> Bayesian linear mixed model [500 draws, exact iid]
#> Formula: Reaction ~ Days_c + (1 + Days_c || Subject)
#>
#> Random effects (variance components fixed at lmer estimates):
#> Groups Name Std.Dev.
#> Subject (Intercept) 37.533
#> Subject.1 Days_c 5.922
#> Residual 25.592
#> Number of obs: 180, groups: Subject, 18
#>
#> --- Posterior means (ICM exact, under fixed variance components) ---
#>
#> RE component parameter fixef.mode fixef.means draws SD
#> ------------ ----------- ------------ ------------ ----------
#> (Intercept) (Intercept) 298.5079 297.9952 9.0581
#> Days_c (Intercept) 10.3748 10.3326 1.4691The header says 500 draws, exact iid. Because
sigma^2 is fixed and both tau^2_k are known,
the joint posterior over the subject coefficients and the population
coefficients is exactly multivariate normal, so the sampler draws from
it directly. There is no burn-in, no autocorrelation, and no convergence
diagnostic to run — the 500 draws are genuinely independent. Chapter 06
covers the case where Gibbs sweeps are needed.
The table compares three quantities per parameter:
| Column | Meaning |
|---|---|
fixef.mode |
Exact posterior mean, computed by iterated conditional modes before any sampling |
fixef.means |
Average of the 500 draws |
draws SD |
Posterior standard deviation from the draws |
The first two should agree to within Monte Carlo error — that is the package’s own correctness check, and a useful one to glance at. Here the average subject reacts about 298 ms at mid-study and slows by roughly 10 ms per additional day of deprivation.
summary() gives the full report: prior versus reference
values, posterior percentiles, tail probabilities, and a per-group
section.
summary(fit)
#> Call:
#> lmerb(formula = form, data = dat, pfamily_list = pf, dispersion_ranef = ps$group.dispersion,
#> n = 500L, progbar = FALSE)
#>
#> Bayesian linear mixed model fit [500 draws, exact iid]
#> Formula: Reaction ~ Days_c + (1 + Days_c || Subject)
#>
#> Random effects (variance components fixed at lmer estimates):
#> Groups Name Std.Dev.
#> Subject (Intercept) 37.533
#> Subject.1 Days_c 5.922
#> Residual 25.592
#> Number of obs: 180, groups: Subject, 18
#>
#> === Block 2 dispersion (RE variance tau^2_k) ===
#>
#> Prior and lmer reference:
#>
#> Prior 1/E[1/tau2] E[tau2] sqrt(E[tau2]) disp_lower disp_upper
#> (Intercept) dNormal 1346.1456 1346.1456 36.6899 NA NA
#> Days_c dNormal 29.3722 29.3722 5.4196 NA NA
#> Residual fixed 709.8260 709.8260 26.6426 NA NA
#> lmer lmer SD
#> (Intercept) 1408.7541 37.5334
#> Days_c 35.0716 5.9221
#> Residual 654.9399 25.5918
#>
#> Overview:
#> Post.Mode Post.Mean Post.Sd MC Error Mean SD Pr(Prior_tail)
#> (Intercept) 1346.15 1346.15 0.00 0.00 36.69 NA
#> Days_c 29.37 29.37 0.00 0.00 5.42 NA
#> Residual 709.83 709.83 NA NA 26.64 NA
#> Cand/draw
#> (Intercept) 1
#> Days_c 1
#> Residual NA
#>
#> Distribution percentiles (tau^2):
#>
#> 1.0% 2.5% 5.0% Median 95.0% 97.5% 99.0%
#> (Intercept) 1346.15 1346.15 1346.15 1346.15 1346.15 1346.15 1346.15
#> Days_c 29.37 29.37 29.37 29.37 29.37 29.37 29.37
#> Residual 709.83 709.83 709.83 709.83 709.83 709.83 709.83
#>
#> SD credible interval (sqrt(tau^2) draws; lmer SD for reference):
#>
#> 2.5% Median 97.5% lmer SD
#> (Intercept) 36.69 36.69 36.69 37.533
#> Days_c 5.42 5.42 5.42 5.922
#> Residual 26.64 26.64 26.64 25.592
#>
#> === Block 2: Level-2 fixed effects (hyperparameters) ===
#>
#> Prior and lmer reference:
#>
#> Prior Mean Prior.sd lmer lmer.se
#> (Intercept)::(Intercept) 298.51 90.05 298.51 9.050
#> Days_c::(Intercept) 0.00 15.38 10.47 1.546
#>
#> Overview:
#> fixef.mode fixef.means Post.Sd MC.Error Pr(Prior_tail)
#> (Intercept)::(Intercept) 298.508 297.995 9.058 0.405 0.438
#> Days_c::(Intercept) 10.375 10.333 1.469 0.066 <2e-16
#>
#> (Intercept)::(Intercept)
#> Days_c::(Intercept) ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Distribution percentiles:
#>
#> 1.0% 2.5% 5.0% Median 95.0% 97.5% 99.0%
#> (Intercept)::(Intercept) 278.379 281.610 284.103 297.280 313.238 316.400 319.07
#> Days_c::(Intercept) 7.147 7.554 8.022 10.323 12.810 13.233 13.71
#>
#> === Block 1: Random effects (group-level) ===
#>
#> Summary of posterior mode (ranef.mode) across groups:
#>
#> Mean SD Min Q1 Median Q3 Max
#> (Intercept) 298.5079 36.4724 219.4041 290.8999 304.3215 314.6673 371.8535
#> Days_c 10.4463 5.0724 0.1222 7.1416 10.4194 12.8294 19.1842
#> MCMC.mean
#> (Intercept) 298.37
#> Days_c 10.46
#>
#> Per-group random effects: inspect fit$groupef.mode or fit$groupef,
#> or call summary(fit, groups = <level ids>) for selected groups.The Pr(Prior_tail) column is a directional tail
probability: how far into the prior’s tail the posterior has moved. A
value near 0.5 means the data did not move that coefficient much; the
near-zero value on Days_c says the sleep-deprivation effect
is strongly supported relative to the prior, which centered it at
zero.
The usual lme4 extractors work, though
fixef() and ranef() are lme4
generics — load lme4 or call them
namespace-qualified.
head(lme4::ranef(fit)$Subject) # per-subject coefficients, chain means
#> (Intercept) Days_c
#> 308 339.7664 19.173801
#> 309 219.3786 4.153587
#> 310 234.4354 7.161023
#> 330 302.3252 4.789262
#> 331 308.6786 6.538323
#> 332 307.4449 9.699681fixef() has three modes, and the default is not the
Bayesian one:
lme4::fixef(fit) # "mer": the reference lmer estimates
#> (Intercept) Days_c
#> 298.50789 10.46729
lme4::fixef(fit, type = "hyper") # Block 2 posterior mode
#> (Intercept).(Intercept) Days_c.(Intercept)
#> 298.50789 10.37476
lme4::fixef(fit, type = "hyper.mean")# Block 2 posterior mean from the draws
#> (Intercept).(Intercept) Days_c.(Intercept)
#> 297.99517 10.33255The default returns the reference lmer() fit’s fixed
effects so that fixef() behaves like it does on a
merMod object. For posterior quantities ask for
"hyper" or "hyper.mean", or read
fit$popef.mode and fit$popef.means
directly.
Note a naming subtlety specific to this package. What
ranef() returns is the full per-subject
coefficient vector, not lme4’s mean-zero deviation. Subject
308’s intercept of roughly 339 is their actual expected reaction time,
not an offset from the population mean of 298. Chapter 02 explains the
parameterization, and Chapter 05 catalogues every accessor.
m_lmer <- lme4::lmer(form, data = dat, REML = TRUE)
round(
cbind(
lmer = lme4::fixef(m_lmer),
lmerb = unlist(fit$popef.means)
),
3
)
#> lmer lmerb
#> (Intercept) 298.508 297.995
#> Days_c 10.467 10.333With a prior weight of one percent, the two agree closely — as they should. The Bayesian fit adds full posterior uncertainty rather than a point estimate with a standard error, and it does so without asymptotic approximation. Chapter 13 works through where the two approaches diverge and why.
||, and how the
parameterization maps to lme4.?lmerb, ?model_setup,
?Prior_Setup_GLMMvignette("Chapter-B00", package = "lmebayesCore") — the
engine roadmap.