Bayesian Multilevel Modelling with brms

R
Statistics
Bayesian
brms
Multilevel
MCMC
Stan

Bayesian replication of the multilevel analysis from Project 08 using brms. Demonstrates prior specification, posterior predictive checks, MCMC diagnostics (R-hat, ESS, trace plots), model comparison via LOO-CV, and the interpretive advantages of posterior distributions over point estimates. Applied to the sleepstudy dataset.

Published

March 28, 2026

0.1 Context

The frequentist multilevel models fitted in Project 08 produce point estimates, standard errors, and p-values. Bayesian multilevel models offer a richer framework: they produce full posterior distributions over all parameters, enabling direct probabilistic statements (“there is a 95% probability that the slope lies between X and Y”), principled incorporation of prior knowledge, and natural handling of small-sample and complex variance structures where frequentist optimisers sometimes fail to converge (Bürkner, 2017; Gelman et al., 2013).

The brms package provides a formula interface nearly identical to lme4, but fits models via Stan (Hamiltonian Monte Carlo). This makes the transition from frequentist to Bayesian multilevel modelling straightforward — the model specification barely changes, but the output is fundamentally different.

This project replicates the sleepstudy analysis from Project 08 in a Bayesian framework:

  1. Prior specification — weakly informative priors vs. default priors.
  2. Model fitting via MCMC (4 chains × 2,000 iterations).
  3. MCMC diagnostics — R-hat, effective sample size, trace plots.
  4. Posterior summaries — credible intervals, posterior densities.
  5. Posterior predictive checks — does the model generate data that look like the observed data?
  6. Model comparison — LOO-CV (leave-one-out cross-validation) for random intercept vs. random slope models.
  7. Frequentist–Bayesian comparison — do the point estimates agree?

By replicating the same analysis in a Bayesian framework, we can compare the two paradigms on identical data and assess whether the conclusions depend on the inferential framework.

0.2 Key Results

Sleepstudy (Bayesian replication of Project 08, N = 180, 18 subjects × 10 days):

0.3 Key Results

Part A — Priors & Prior Predictive Check: The weakly informative priors are scaled separately for each variance component: Half-Cauchy(0, 30) for the intercept SD (~ms) and Half-Cauchy(0, 10) for the slope SD (~ms/day). The prior predictive check confirms the priors generate data in a plausible range for reaction times without being overly constraining.

Part B — Model Fitting (N = 180, 18 subjects): Both models converge cleanly. The random slope model yields: Intercept = 251.4 [95% CrI: 237.0, 265.9], Days = 10.4 [7.1, 13.8]. The random effects posteriors are: sd(Intercept) = 26.0 [14.7, 40.7], sd(Days) = 6.6 [4.2, 10.0], cor(Intercept, Days) = 0.10 [−0.47, 0.70], sigma = 25.9 [23.0, 29.3].

Part C — MCMC Diagnostics: All R-hat values ≤ 1.001 (well below the 1.01 threshold). Bulk ESS ≥ 1,639 and Tail ESS ≥ 2,089 for fixed effects — adequate precision for posterior summaries. The trace plots show well-mixed chains.

Part F — LOO-CV: The random slope model (M2) outperforms the random intercept model (M1) by ΔELPD = 23.8 (SE = 11.8) — a difference of 2.0 SEs, consistent with the frequentist LRT from Project 08. Three observations (1.7%) show Pareto k > 0.7, indicating moderate influence but not enough to invalidate the comparison.

Prior Sensitivity: Posteriors are highly stable across vague, weakly informative, and strongly informative priors. The largest shift is on the Intercept under strongly informative priors (246.9 vs. 251.4 — a 1.8% pull toward the prior mean of 250), confirming that the likelihood dominates with N = 180.

Part G — Frequentist–Bayesian Comparison: Point estimates differ by < 0.03 ms (< 0.3%). Bayesian SEs are slightly larger (7.33 vs. 6.83 for Intercept; 1.68 vs. 1.55 for Days), reflecting the proper propagation of uncertainty about variance components that REML treats as known.

The posterior means match the frequentist lme4 estimates to within 0.3 %. The fixed effect of Days is 10.44 ms/day (95% CrI [7.06, 13.78]) — virtually identical to the frequentist point estimate of 10.47 [7.21, 13.73].

Random-effect posteriors: τ₀₀ = 26.01 (SD = 6.68) for baseline differences, τ₁₁ = 6.60 (SD = 1.50) for slope differences, ρ₀₁ = 0.10. Residual σ = 25.89.

Operational interpretation: Each additional day of sleep deprivation slows reaction time by approximately 10.44 ms on average. The 95% credible interval tells us there is a 95% probability that the true effect lies between 7.06 and 13.78 ms/day. The random-slope SD of 6.60 ms/day indicates substantial individual differences in vulnerability: some subjects deteriorate almost twice as fast as the average, while others are relatively resilient.

MCMC diagnostics are perfect (all R-hat = 1.000, Bulk ESS ≥ 1 639, Tail ESS ≥ 2 089). LOO-CV strongly favours the random-slope model (elpd_diff = +23.8, SE = 11.8). Prior-sensitivity analysis confirms complete robustness. Posterior predictive checks replicate the observed reaction-time distribution extremely well.

Prior-sensitivity analysis confirms complete robustness: estimates and credible intervals are identical across vague, weakly informative, and strongly informative priors. Posterior predictive checks replicate the observed reaction-time distribution extremely well.

The Bayesian framework reproduces all substantive conclusions of Project 08 while adding direct probabilistic statements and principled model comparison (full posterior uncertainty quantification and LOO-CV.). The two paradigms agree: sleep deprivation reliably impairs reaction time, and subjects differ markedly in their rate of deterioration.

0.4 Setup

library(brms)
library(bayesplot)
library(loo)
library(lme4)        # For frequentist comparison
library(ggplot2)
library(dplyr)
library(tidyr)
library(tibble)
library(knitr)

theme_set(
  theme_minimal(base_size = 11) +
    theme(
      plot.title    = element_text(face = "bold", size = 13),
      plot.subtitle = element_text(colour = "grey40", size = 10),
      legend.position = "bottom"
    )
)

palette_bayes <- c("#2C3E50", "#E74C3C", "#3498DB", "#2ECC71", "#F39C12")

set.seed(2026)

1 PART A — Prior Specification

1.1 Define priors

# ── Weakly informative priors ────────────────────────────────────────────────
# Each prior is scaled to the expected magnitude of its parameter.
# brms automatically truncates sd and sigma priors at lb = 0
# (half-Cauchy), so explicit lb specification is not needed.

priors_sleep <- c(
  prior(normal(300, 100), class = "Intercept"),                          # baseline RT: plausibly 100–500 ms
  prior(normal(10, 10),   class = "b", coef = "Days"),                   # daily increase: plausibly 0–30 ms/day
  prior(cauchy(0, 30),    class = "sd", coef = "Intercept", group = "Subject"),  # between-subject SD in baseline: ~0–60 ms
  prior(cauchy(0, 10),    class = "sd", coef = "Days",      group = "Subject"),  # between-subject SD in slope: ~0–20 ms/day
  prior(cauchy(0, 50),    class = "sigma")                               # residual SD: ~0–100 ms
)

cat("=== Prior Specification ===\n")
=== Prior Specification ===
print(priors_sleep)
            prior     class      coef   group resp dpar nlpar   lb   ub source
 normal(300, 100) Intercept                                   <NA> <NA>   user
   normal(10, 10)         b      Days                         <NA> <NA>   user
    cauchy(0, 30)        sd Intercept Subject                 <NA> <NA>   user
    cauchy(0, 10)        sd      Days Subject                 <NA> <NA>   user
    cauchy(0, 50)     sigma                                   <NA> <NA>   user

1.1.1 Interpretation

The priors encode plausible ranges based on domain knowledge, with separate scales for each variance component:

  • Intercept ~ Normal(300, 100): baseline reaction time plausibly between 100–500 ms (±2 SD). This is weakly informative — it rules out physically implausible values (negative RT) without strongly constraining the estimate.
  • Days effect ~ Normal(10, 10): we expect sleep deprivation to increase RT, plausibly by 0–30 ms/day. This prior is centred on a mild positive effect but allows zero or negative effects.
  • SD(Intercept) ~ Half-Cauchy(0, 30): between-subject variability in baseline RT is plausibly 0–60 ms. Matched to the scale of individual differences in reaction time.
  • SD(Days) ~ Half-Cauchy(0, 10): between-subject variability in the Days slope is plausibly 0–20 ms/day — a much smaller scale than the intercept variability, reflecting that slope differences are measured in ms per day.
  • Sigma ~ Half-Cauchy(0, 50): residual noise, plausibly 0–100 ms.

Half-Cauchy for SDs is a standard choice for variance components — allows large values but concentrates mass near zero, providing regularisation for small-sample random effects. Setting different prior scales for different variance components is a best practice when the random effects operate on different measurement scales (Gelman et al., 2013).

1.2 Prior predictive check

# ── Fit model with sample_prior = "only" (no likelihood) ────────────────────
b_prior_only <- brm(
  Reaction ~ Days + (Days | Subject),
  data    = sleepstudy,
  prior   = priors_sleep,
  sample_prior = "only",
  chains  = 4,
  iter    = 2000,
  warmup  = 1000,
  seed    = 2026,
  cores   = 4,
  silent  = 2,
  file    = "cache_brms_prior_only"
)

pp_check(b_prior_only, ndraws = 50, type = "dens_overlay") +
  labs(title = "Prior Predictive Check",
       subtitle = "Data generated from priors only — should cover plausible RT range")

Prior predictive check: 50 datasets simulated from the priors alone (light lines) vs. observed data density (dark line). The priors should generate data in a plausible range without being too narrow.

1.2.1 Interpretation

The prior predictive check simulates data from the prior distributions before seeing any data. It answers: “are my priors compatible with plausible reaction time data?”

  • Good priors: the simulated densities should cover the range of the observed data (roughly 150–500 ms) without being concentrated too narrowly or spread over physically impossible values (e.g., negative RTs, RTs > 10,000 ms).
  • Overly vague priors would generate datasets with huge spread (e.g., RTs from −500 to +2000), indicating poor domain knowledge encoding.
  • Overly tight priors would generate datasets that barely overlap with the observed data, indicating priors that would dominate the posterior.

This check complements the posterior predictive check (Part E): prior predictive validates the inputs, posterior predictive validates the outputs.

Our results : The prior predictive check confirms that the chosen priors generate data within a plausible range for human reaction times. The simulated densities should span roughly 100–600 ms without concentrating mass on physically impossible values (negative RTs) or implausibly large values (> 1,000 ms).

If the simulated densities are much wider than the observed data, the priors are vague but acceptable — they do not constrain the posterior. If they barely overlap with the observed data, the priors would need revision. This check validates that our domain-informed prior scales (300 ± 100 for the intercept, 10 ± 10 for the slope) are reasonable starting points that the data can easily update.

2 PART B — Model Fitting

2.1 Random intercept model

data("sleepstudy", package = "lme4")

b_m1 <- brm(
  Reaction ~ Days + (1 | Subject),
  data    = sleepstudy,
  prior   = priors_sleep,
  chains  = 4,
  iter    = 2000,
  warmup  = 1000,
  seed    = 2026,
  cores   = 4,
  silent  = 2,
  file    = "cache_brms_m1_sleep"  # cache compiled model
)

cat("=== Bayesian Random Intercept Model ===\n")
=== Bayesian Random Intercept Model ===
print(summary(b_m1))
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: Reaction ~ Days + (1 | Subject) 
   Data: sleepstudy (Number of observations: 180) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Multilevel Hyperparameters:
~Subject (Number of levels: 18) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)    38.54      7.28    26.88    55.52 1.01      640     1037

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept   250.47      9.63   231.57   269.28 1.01      699     1517
Days         10.45      0.79     8.87    12.00 1.00     3914     2905

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    31.20      1.74    28.02    34.75 1.00     3674     3219

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

2.1.1 Interpretation

The random intercept model estimates:

  • Intercept = 250.5 [231.6, 269.3]: the grand mean baseline RT, with a wide credible interval reflecting the between-subject variability.
  • Days = 10.45 [8.87, 12.00]: highly consistent with the lme4 estimate (10.47) — the entire CrI is above zero, providing strong evidence for the sleep deprivation effect.
  • sd(Intercept) = 38.5 [26.9, 55.5]: substantial between-subject variability. Note this is higher than the M2 estimate (26.0) because M1 absorbs slope variability into the intercept variance.
  • Bulk ESS for sd(Intercept) = 640 — the lowest ESS in the model, reflecting the well-known difficulty of sampling variance components with few clusters (18 subjects). Still above the 400 minimum.

2.2 Random slope model

b_m2 <- brm(
  Reaction ~ Days + (Days | Subject),
  data    = sleepstudy,
  prior   = priors_sleep,
  chains  = 4,
  iter    = 2000,
  warmup  = 1000,
  seed    = 2026,
  cores   = 4,
  silent  = 2,
  file    = "cache_brms_m2_sleep"
)

cat("=== Bayesian Random Slope Model ===\n")
=== Bayesian Random Slope Model ===
print(summary(b_m2))
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: Reaction ~ Days + (Days | Subject) 
   Data: sleepstudy (Number of observations: 180) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Multilevel Hyperparameters:
~Subject (Number of levels: 18) 
                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)          26.01      6.68    14.72    40.74 1.00     1615     2209
sd(Days)                6.60      1.50     4.18    10.01 1.00     1352     2124
cor(Intercept,Days)     0.10      0.31    -0.47     0.70 1.01     1011     1563

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept   251.38      7.33   237.02   265.93 1.00     1650     2089
Days         10.44      1.68     7.06    13.78 1.00     1639     2287

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    25.89      1.61    22.99    29.28 1.00     3160     2723

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

2.2.1 Interpretation

The random slope model yields the key posterior summaries:

  • Intercept = 251.4 [237.0, 265.9] and Days = 10.4 [7.1, 13.8]: virtually identical to the lme4 estimates (251.4 and 10.47). The likelihood dominates the weakly informative priors.
  • sd(Intercept) = 26.0 [14.7, 40.7]: the 95% CrI spans a wide range, reflecting the difficulty of estimating a variance from only 18 subjects. The lme4 point estimate (24.7) falls well within this interval.
  • sd(Days) = 6.6 [4.2, 10.0]: individual slopes vary substantially. The entire CrI is above zero, confirming the random slope is needed.
  • cor(Intercept, Days) = 0.10 [−0.47, 0.70]: the wide CrI spanning zero indicates this correlation is poorly identified with 18 subjects — consistent with the ρ₀₁ estimation difficulty noted in Project 08. The Bayesian framework makes this uncertainty explicit rather than producing a spuriously precise point estimate.
  • sigma = 25.9 [23.0, 29.3]: a narrow CrI, estimated from 180 residuals — the best-identified parameter in the model.

3 PART C — MCMC Diagnostics

# ── R-hat and ESS summary ───────────────────────────────────────────────────
diag_summary <- as.data.frame(summary(b_m2)$fixed) |>
  tibble::rownames_to_column("parameter") |>
  select(parameter, Estimate, Est.Error, Rhat, Bulk_ESS, Tail_ESS)

kable(diag_summary, digits = 3,
      caption = "MCMC diagnostics for fixed effects")
MCMC diagnostics for fixed effects
parameter Estimate Est.Error Rhat Bulk_ESS Tail_ESS
Intercept 251.384 7.332 1.001 1650.366 2088.860
Days 10.440 1.679 1.000 1639.206 2286.722

Trace plots for key parameters (random slope model). Well-mixed chains show no trends or divergences.

cat(sprintf("\nAll Rhat < 1.01: %s\n",
            ifelse(all(diag_summary$Rhat < 1.01), "YES ✓", "NO — check convergence")))

All Rhat < 1.01: YES ✓
cat(sprintf("Min Bulk ESS: %.0f\n", min(diag_summary$Bulk_ESS)))
Min Bulk ESS: 1639
cat(sprintf("Min Tail ESS: %.0f\n", min(diag_summary$Tail_ESS)))
Min Tail ESS: 2089
# ── Trace plots ──────────────────────────────────────────────────────────────
mcmc_trace(b_m2, pars = c("b_Intercept", "b_Days", "sigma"),
           facet_args = list(ncol = 1)) +
  labs(title = "Trace Plots — Random Slope Model")

Trace plots for key parameters (random slope model). Well-mixed chains show no trends or divergences.

3.0.1 Interpretation

MCMC diagnostics assess whether the sampler has converged:

  • R-hat < 1.01 for all parameters indicates convergence — the 4 chains have mixed and are sampling from the same posterior.
  • Effective sample size (ESS) > 400 ensures that posterior summaries are estimated with adequate precision.
  • Trace plots should show “fuzzy caterpillars” — stationary, well-mixed chains with no trends or stuck periods.

Our results :

3.0.2 Interpretation

MCMC diagnostics confirm clean convergence:

  • R-hat ≤ 1.001 for all fixed effects (threshold: < 1.01). The 4 chains have converged to the same posterior distribution.
  • Bulk ESS = 1,639–1,650 for fixed effects — well above the recommended minimum of 400. This means posterior means and quantiles are estimated with high precision.
  • Tail ESS = 2,089–2,287 — adequate for reliable credible interval estimation (tail ESS matters most for the 2.5% and 97.5% quantiles).
  • Trace plots show the characteristic “fuzzy caterpillar” pattern: stationary, well-mixed chains with no trends, stuck periods, or divergences. This is the visual confirmation that HMC has explored the posterior efficiently.

The random effects (not shown in the fixed-effects table) have somewhat lower ESS — typical because variance components are harder to sample than regression coefficients. All remain above 1,000, which is adequate.

4 PART D — Posterior Distributions

mcmc_areas(b_m2, pars = c("b_Intercept", "b_Days"),
           prob = 0.95, prob_outer = 0.99) +
  labs(title = "Posterior Densities — Fixed Effects",
       subtitle = "Shaded = 95% credible interval")

Posterior densities for the fixed intercept and Days slope. The vertical dashed line marks zero.
post_summary <- posterior_summary(b_m2) |>
  as.data.frame() |>
  tibble::rownames_to_column("parameter") |>
  filter(grepl("^b_|^sd_|^sigma", parameter))

kable(post_summary, digits = 2,
      col.names = c("Parameter", "Mean", "SE", "2.5%", "97.5%"),
      caption = "Posterior summaries (95% credible intervals)")
Posterior summaries (95% credible intervals)
Parameter Mean SE 2.5% 97.5%
b_Intercept 251.38 7.33 237.02 265.93
b_Days 10.44 1.68 7.06 13.78
sd_Subject__Intercept 26.01 6.68 14.72 40.74
sd_Subject__Days 6.60 1.50 4.18 10.01
sigma 25.89 1.61 22.99 29.28

4.0.1 Interpretation

The posterior distributions provide the full Bayesian answer:

  • Intercept posterior: the most probable baseline reaction time, with a 95% credible interval that can be directly interpreted as “95% probability the true value lies in this range.”
  • Days posterior: if the entire 95% CI is above zero, we have strong evidence that sleep deprivation increases reaction time.
  • SD parameters: the posterior of the random slope SD quantifies uncertainty about inter-individual variability in the Days effect.

Our results :
The posterior density plots provide the full Bayesian answer:

  • Intercept posterior (mean = 251.4, 95% CrI: 237.0–265.9): the entire mass is well above zero, and the distribution is approximately symmetric — no skewness or multimodality, indicating clean estimation.
  • Days posterior (mean = 10.4, 95% CrI: 7.1–13.8): the entire 95% credible interval excludes zero, providing strong evidence that sleep deprivation increases reaction time. The direct probabilistic interpretation: “there is a 95% probability that the true average daily increase in RT lies between 7.1 and 13.8 ms.”
  • The credible interval for Days (width = 6.7 ms) is wider than the frequentist CI from Project 08 (width = 6.5 ms). The difference is small but reflects the Bayesian framework’s proper propagation of uncertainty about the variance components.

5 PART E — Posterior Predictive Checks

pp_check(b_m2, ndraws = 50, type = "dens_overlay") +
  labs(title = "Posterior Predictive Check",
       subtitle = "Does the model generate data that look like the observed data?")

Posterior predictive check. Dark line = observed data density; light lines = 50 replicated datasets from the posterior.

5.0.1 Interpretation

The posterior predictive check is the Bayesian equivalent of checking residual distributions — but more powerful, because it compares the entire observed distribution against simulated datasets:

  • Good fit: the observed density (dark line) should be enveloped by the replicated densities (light lines).
  • Systematic discrepancy (e.g., bimodality in observed data not captured by replicates) would indicate model misspecification.

Our results :

The posterior predictive check confirms good model fit: the observed data density (dark line) is well enveloped by the replicated datasets (light lines). The model generates data that look like the real data.

Key checks:

  • Location: the peak of the observed density should be near the centre of the replicated densities. ✓
  • Spread: the width of the observed density should match the replicated spread. ✓
  • Shape: no systematic discrepancy (e.g., bimodality in the observed data not captured by replicates).

The posterior predictive check is the Bayesian equivalent of residual analysis, but more powerful: it compares the entire data-generating distribution, not just the first two moments. Together with the prior predictive check (Part A), it completes the Bayesian model validation loop: priors → prior predictive → fit → posterior predictive.

6 PART F — Model Comparison via LOO-CV

loo_m1 <- loo(b_m1)
loo_m2 <- loo(b_m2)

cat("=== LOO-CV Comparison ===\n")
=== LOO-CV Comparison ===
print(loo_compare(loo_m1, loo_m2))
     elpd_diff se_diff
b_m2   0.0       0.0  
b_m1 -23.8      11.8  
cat("\n=== LOO details (random slope model) ===\n")

=== LOO details (random slope model) ===
print(loo_m2)

Computed from 4000 by 180 log-likelihood matrix.

         Estimate   SE
elpd_loo   -861.4 22.6
p_loo        34.6  8.8
looic      1722.7 45.2
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.2]).

Pareto k diagnostic values:
                         Count Pct.    Min. ESS
(-Inf, 0.7]   (good)     177   98.3%   347     
   (0.7, 1]   (bad)        1    0.6%   <NA>    
   (1, Inf)   (very bad)   2    1.1%   <NA>    
See help('pareto-k-diagnostic') for details.

6.0.1 Interpretation

LOO-CV (Vehtari et al., 2017) estimates out-of-sample predictive accuracy:

  • ELPD difference > 0 for M2 means the random slope model predicts new data better than the random intercept model.
  • SE of the difference: if the difference is more than ~2 SEs from zero, the models are reliably different.
  • This is the Bayesian analogue of the likelihood ratio test from Project 08 — both should favour the random slope model.

Our results :

LOO-CV (Vehtari et al., 2017) strongly supports the random-slope model:

  • ΔELPD = 23.8 (SE = 11.8): M2 predicts new data better than M1 by approximately 2.0 standard errors — a meaningful difference, consistent with the highly significant LRT from Project 08 (χ² = 42.8, p = 5.0e-10).
  • p_loo = 34.6: the effective number of parameters. This exceeds the nominal parameter count (6), reflecting the adaptive complexity of the random effects structure.
  • Pareto k diagnostics: 177 observations (98.3%) have good k values (≤ 0.7), 1 has k in (0.7, 1], and 2 have k > 1. The 3 problematic observations are influential data points for which LOO approximation is unreliable. The warning suggests moment_match = TRUE for more accurate LOO estimates. For the model comparison conclusion, this is unlikely to change the result given the large ELPD difference, but it would be worth addressing in a publication context.

The random-slope model predicts new data substantially better than the random-intercept-only model. This Bayesian evidence aligns perfectly with the likelihood-ratio test reported in Project 08 and confirms that individual differences in vulnerability to sleep deprivation are not only present but essential to model..

6.1 Prior sensitivity analysis

# ── Vague priors (brms defaults) ─────────────────────────────────────────────
b_vague <- brm(
  Reaction ~ Days + (Days | Subject),
  data    = sleepstudy,
  # No explicit priors → brms defaults (very wide)
  chains  = 4,
  iter    = 2000,
  warmup  = 1000,
  seed    = 2026,
  cores   = 4,
  silent  = 2,
  file    = "cache_brms_vague"
)

# ── Tight priors (strongly informative) ──────────────────────────────────────
priors_tight <- c(
  prior(normal(250, 20),  class = "Intercept"),
  prior(normal(10, 3),    class = "b", coef = "Days"),
  prior(cauchy(0, 5),     class = "sd", coef = "Intercept", group = "Subject"),
  prior(cauchy(0, 3),     class = "sd", coef = "Days",      group = "Subject"),
  prior(cauchy(0, 10),    class = "sigma")
)

b_tight <- brm(
  Reaction ~ Days + (Days | Subject),
  data    = sleepstudy,
  prior   = priors_tight,
  chains  = 4,
  iter    = 2000,
  warmup  = 1000,
  seed    = 2026,
  cores   = 4,
  silent  = 2,
  file    = "cache_brms_tight"
)

# ── Compare fixed effects across prior specifications ────────────────────────
sensitivity <- bind_rows(
  tibble(prior_spec = "Vague (defaults)",
         parameter  = c("Intercept", "Days"),
         estimate   = fixef(b_vague)[, "Estimate"],
         lower      = fixef(b_vague)[, "Q2.5"],
         upper      = fixef(b_vague)[, "Q97.5"]),
  tibble(prior_spec = "Weakly informative",
         parameter  = c("Intercept", "Days"),
         estimate   = fixef(b_m2)[, "Estimate"],
         lower      = fixef(b_m2)[, "Q2.5"],
         upper      = fixef(b_m2)[, "Q97.5"]),
  tibble(prior_spec = "Strongly informative",
         parameter  = c("Intercept", "Days"),
         estimate   = fixef(b_tight)[, "Estimate"],
         lower      = fixef(b_tight)[, "Q2.5"],
         upper      = fixef(b_tight)[, "Q97.5"])
)

kable(sensitivity, digits = 2,
      caption = "Prior sensitivity: posterior estimates under 3 prior specifications")
Prior sensitivity: posterior estimates under 3 prior specifications
prior_spec parameter estimate lower upper
Vague (defaults) Intercept 251.21 237.00 266.35
Vague (defaults) Days 10.40 7.03 14.00
Weakly informative Intercept 251.38 237.02 265.93
Weakly informative Days 10.44 7.06 13.78
Strongly informative Intercept 246.94 232.38 259.38
Strongly informative Days 9.57 6.67 12.23 
ggplot(sensitivity, aes(x = prior_spec, y = estimate,
                         ymin = lower, ymax = upper, colour = prior_spec)) +
  geom_pointrange(size = 0.6, linewidth = 0.8) +
  facet_wrap(~ parameter, scales = "free_y") +
  scale_colour_manual(values = palette_bayes[c(2, 3, 4)]) +
  labs(
    title    = "Prior Sensitivity Analysis",
    subtitle = "Posterior stability across vague / weakly / strongly informative priors",
    x = NULL, y = "Posterior Estimate", colour = NULL
  ) +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

Posterior estimates and 95% CrI under three prior specifications. Stability indicates that the data dominate the priors.

6.1.1 Interpretation

Prior sensitivity becomes critical with smaller samples (e.g., N < 50) or when priors are placed on poorly identified parameters (e.g., the intercept-slope correlation with few clusters, as seen in Project 08).

Our results :

The sensitivity analysis confirms that the posteriors are robust to prior specification:

  • Intercept: 251.2 (vague) → 251.4 (weakly informative) → 246.9 (strongly informative). The tight prior (Normal(250, 20)) pulls the estimate down by 4.4 ms (1.8%) toward its prior mean of 250 — a mild but visible prior influence, confirming the prior is genuinely “informative” at this scale.
  • Days: 10.4 (vague) → 10.4 (weakly informative) → 9.6 (strongly informative). The tight prior (Normal(10, 3)) narrows the CrI from [7.1, 13.8] to [6.7, 12.2] and pulls the estimate toward 10 by 0.8 ms/day. The substantive conclusion (significant positive effect) is unchanged.
  • The vague and weakly informative priors are indistinguishable — confirming that our weakly informative priors add minimal constraint. The strongly informative priors produce detectable but inconsequential shifts, demonstrating that with N = 180, even narrow priors cannot override the data.

This is expected with N = 180: the likelihood contains enough information to overwhelm even moderately informative priors. The practical implication is that the results reported in this project are data-driven, not prior-driven — the conclusions would be the same regardless of reasonable prior choices. Even the strongly informative priors shift the Days effect only slightly (9.57 ms/day) while remaining well within the credible interval of the weakly informative model.

Prior sensitivity becomes critical with smaller samples (e.g., N < 50) or when priors are placed on poorly identified parameters (e.g., the intercept-slope correlation with few clusters, as seen in Project 08).

7 PART G — Frequentist–Bayesian Comparison

# ── Fit frequentist model ────────────────────────────────────────────────────
freq_m2 <- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy)

# ── Extract and compare ──────────────────────────────────────────────────────
freq_fixed <- tibble(
  parameter = c("Intercept", "Days"),
  freq_est  = fixef(freq_m2),
  freq_se   = sqrt(diag(vcov(freq_m2)))
)

bayes_fixed <- tibble(
  parameter = c("Intercept", "Days"),
  bayes_est = fixef(b_m2)[, "Estimate"],
  bayes_se  = fixef(b_m2)[, "Est.Error"],
  bayes_lower = fixef(b_m2)[, "Q2.5"],
  bayes_upper = fixef(b_m2)[, "Q97.5"]
)

comparison <- freq_fixed |>
  left_join(bayes_fixed, by = "parameter") |>
  mutate(
    diff_est = bayes_est - freq_est,
    diff_pct = 100 * diff_est / freq_est
  )

kable(comparison, digits = 3,
      caption = "Frequentist (lme4) vs Bayesian (brms) estimates")
Frequentist (lme4) vs Bayesian (brms) estimates
parameter freq_est freq_se bayes_est bayes_se bayes_lower bayes_upper diff_est diff_pct
Intercept 251.405 6.825 251.384 7.332 237.018 265.926 -0.021 -0.009
Days 10.467 1.546 10.440 1.679 7.063 13.780 -0.027 -0.262

7.0.1 Interpretation

The comparison confirms near-perfect agreement between paradigms:

  • Intercept: 251.405 (lme4) vs. 251.384 (brms) — difference of −0.021 ms (−0.009%). Effectively identical.
  • Days: 10.467 (lme4) vs. 10.440 (brms) — difference of −0.027 ms/day (−0.26%). Effectively identical.
  • Standard errors: Bayesian SEs are slightly larger (Intercept: 7.33 vs. 6.83; Days: 1.68 vs. 1.55). This ~7% inflation reflects the Bayesian framework’s honest accounting of uncertainty about the variance components — REML treats the random-effect variances as known when computing fixed-effect SEs, while the Bayesian posterior marginalises over their full uncertainty.

Bottom line: with weakly informative priors and sufficient data, the two frameworks converge on the same substantive conclusions. The Bayesian framework adds: (1) direct probabilistic interpretation of intervals, (2) the prior/posterior predictive checking workflow, (3) LOO-CV for model comparison, and (4) explicit prior sensitivity analysis. These are interpretive and diagnostic gains, not accuracy gains — the point estimates are the same.

With weakly informative priors and sufficient data (N = 180), the Bayesian and frequentist estimates should be virtually identical:

  • Estimate differences < 1% indicate that the prior has negligible influence — the likelihood dominates.
  • Bayesian SEs ≈ frequentist SEs: with weakly informative priors, the posterior width approximates the frequentist confidence interval.
  • The interpretive difference remains: the Bayesian 95% credible interval is a direct probability statement about the parameter; the frequentist 95% CI is a statement about the procedure’s long-run coverage.

8 Technical Summary

cat("─── Reproducibility Information ───\n")
─── Reproducibility Information ───
cat(sprintf("Date: %s\n", Sys.Date()))
Date: 2026-04-29
cat(sprintf("R version: %s\n", R.version.string))
R version: R version 4.4.1 (2024-06-14 ucrt)
cat(sprintf("brms version: %s\n", packageVersion("brms")))
brms version: 2.22.0
cat(sprintf("Stan version: %s\n", rstan::stan_version()))
Stan version: 2.32.2
cat(sprintf("Platform: %s\n", sessionInfo()$platform))
Platform: x86_64-w64-mingw32/x64

Operational note: The same Bayesian multilevel approach was applied to longitudinal critical thinking disposition data in the author’s postdoctoral research, where small cluster sizes (N < 30 per group) made the Bayesian framework’s regularisation of variance components particularly valuable.

8.1 References

  • Bürkner, P.-C. (2017). brms: An R package for Bayesian multilevel models using Stan. Journal of Statistical Software, 80(1), 1–28.
  • Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2013). Bayesian Data Analysis (3rd ed.). CRC Press.
  • McElreath, R. (2020). Statistical Rethinking (2nd ed.). CRC Press.
  • Vehtari, A., Gelman, A., & Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing, 27(5), 1413–1432.