Bayesian Mixed Effects Models with brms for Linguists
You’ve fitted a Bayesian model with brms. The code ran without errors. You got parameter estimates and credible intervals. But before you interpret or report these results, you need to ask:
“Did the MCMC sampler actually converge to the posterior distribution?”
If the chains haven’t converged, your estimates are unreliable — no matter how plausible they look!
MCMC (Markov Chain Monte Carlo) is a sampling algorithm that explores the posterior distribution. Think of it like:
Convergence means all chains have “forgotten” their starting points and are sampling from the same distribution.
Today we’ll learn to check convergence using:
Module 01-02: Build model + set priors
↓
Module 03: Check model fit (posterior predictive)
↓
Module 04: Test prior sensitivity
↓
Module 05: Compare models (LOO-CV for prediction)
↓
Module 06: Practical significance (ROPE, emmeans, marginaleffects)
↓
Module 07: Hypothesis comparison (Bayes Factors)
↓
Module 08 (TODAY): Convergence diagnostics
↓
Report results!Answer: ALWAYS, before interpreting any results!
Typical workflow:
For each model, we’ll ask:
We’ll also test: 6. ✅ Stability test: Does doubling iterations change our conclusions?
Let’s work through these with real examples.
library(brms) # Bayesian regression models
library(tidyverse) # Data manipulation
library(bayesplot) # MCMC diagnostic plots
library(posterior) # Working with posterior draws
library(tidybayes) # Tidy MCMC tools
library(patchwork) # Combining plots
# Set plotting theme
theme_set(theme_minimal(base_size = 14))
# Set seed for reproducibility
set.seed(42)
# Bayesplot options
color_scheme_set("viridis")We’ll use the same data structure as Modules 06 and 07: a psycholinguistic experiment with reaction times.
# Parameters
n_subjects <- 40
n_items <- 60
n_per_cell <- 3
# Generate data
rt_data <- expand_grid(
subject = factor(1:n_subjects),
item = factor(1:n_items),
condition = factor(c("A", "B", "C"))
) %>%
slice_sample(n = n_subjects * n_items * n_per_cell) %>%
# Add random effects
group_by(subject) %>%
mutate(
subj_intercept = rnorm(1, mean = 0, sd = 0.15),
subj_slope_B = rnorm(1, mean = 0, sd = 0.08),
subj_slope_C = rnorm(1, mean = 0, sd = 0.08)
) %>%
group_by(item) %>%
mutate(
item_intercept = rnorm(1, mean = 0, sd = 0.10)
) %>%
ungroup() %>%
# Generate log-RT
mutate(
condition_effect_B = if_else(condition == "B", 0.12, 0), # 12% slower
condition_effect_C = if_else(condition == "C", 0.18, 0), # 18% slower
log_rt = 6.0 + # baseline ≈ 400ms
subj_intercept +
item_intercept +
condition_effect_B +
condition_effect_C +
(condition == "B") * subj_slope_B +
(condition == "C") * subj_slope_C +
rnorm(n(), mean = 0, sd = 0.20), # residual noise
rt = exp(log_rt)
) %>%
select(subject, item, condition, log_rt, rt)
# Data summary
cat("=== DATA SUMMARY ===\n")=== DATA SUMMARY ===rt_data %>%
group_by(condition) %>%
summarise(
n = n(),
mean_rt = mean(rt) %>% round(0),
median_rt = median(rt) %>% round(0),
sd_rt = sd(rt) %>% round(0),
mean_log_rt = mean(log_rt) %>% round(3),
.groups = "drop"
) %>%
print()# A tibble: 3 × 6
condition n mean_rt median_rt sd_rt mean_log_rt
<fct> <int> <dbl> <dbl> <dbl> <dbl>
1 A 2400 414 400 117 5.99
2 B 2400 474 454 145 6.12
3 C 2400 509 480 154 6.19cat("\nCondition B - A effect (log):",
round(mean(rt_data$log_rt[rt_data$condition == "B"]) -
mean(rt_data$log_rt[rt_data$condition == "A"]), 3), "\n")
Condition B - A effect (log): 0.13 cat("Condition C - A effect (log):",
round(mean(rt_data$log_rt[rt_data$condition == "C"]) -
mean(rt_data$log_rt[rt_data$condition == "A"]), 3), "\n")Condition C - A effect (log): 0.203 p1 <- ggplot(rt_data, aes(x = condition, y = rt, fill = condition)) +
geom_violin(alpha = 0.6) +
geom_jitter(width = 0.2, alpha = 0.3, size = 0.5) +
stat_summary(fun = mean, geom = "point", size = 3, color = "red") +
stat_summary(fun.data = mean_cl_boot, geom = "errorbar",
width = 0.2, color = "red", linewidth = 1) +
labs(title = "Reaction Times by Condition",
y = "RT (ms)", x = "Condition") +
theme(legend.position = "none")
p2 <- ggplot(rt_data, aes(x = condition, y = log_rt, fill = condition)) +
geom_violin(alpha = 0.6) +
geom_jitter(width = 0.2, alpha = 0.3, size = 0.5) +
stat_summary(fun = mean, geom = "point", size = 3, color = "red") +
stat_summary(fun.data = mean_cl_boot, geom = "errorbar",
width = 0.2, color = "red", linewidth = 1) +
labs(title = "Log-Transformed RTs",
y = "log(RT)", x = "Condition") +
theme(legend.position = "none")
p1 + p2
Let’s start with a model that should converge well.
# Define priors
rt_priors <- c(
prior(normal(6, 1), class = Intercept), # Baseline around 400ms
prior(normal(0, 0.5), class = b), # Effects typically < 50% change
prior(exponential(1), class = sd), # Moderate random effects
prior(exponential(2), class = sigma), # Residual noise
prior(lkj(2), class = cor) # Slight correlation regularization
)
# Fit model with default settings
rt_model_good <- brm(
log_rt ~ condition + (1 + condition | subject) + (1 | item),
data = rt_data,
family = gaussian(),
prior = rt_priors,
iter = 2000,
warmup = 1000,
chains = 4,
cores = 4,
seed = 42,
backend = "cmdstanr",
silent = 2,
refresh = 0
)A trace plot shows parameter values across MCMC iterations for each chain.
Good trace plot (well-mixed chains): - All chains overlap completely - Look like “hairy caterpillars” - No trends or drift over time - Chains quickly forget their starting points
Bad trace plot (convergence problems): - Chains explore different regions - Systematic trends over time - One chain stuck while others move - Slow wandering (high autocorrelation)
# Extract posterior draws
posterior_good <- as_draws_array(rt_model_good)
# Plot traces for key parameters
mcmc_trace(
posterior_good,
pars = c("b_Intercept", "b_conditionB", "b_conditionC",
"sd_subject__Intercept", "sigma"),
facet_args = list(ncol = 1)
) +
labs(title = "Trace Plots: Well-Converged Model",
subtitle = "All chains mix well and explore the same space")
Good signs: - ✅ All 4 chains completely overlap - ✅ No systematic drift up or down - ✅ “Fuzzy caterpillar” appearance - ✅ Chains explore same range of values
Warning signs: - ❌ Chains separated vertically (exploring different modes) - ❌ Trends over time (non-stationarity) - ❌ One chain behaves differently - ❌ Slow, snake-like wandering
# Get all parameter names (exclude lp__ and individual random effects)
all_pars <- variables(posterior_good)
focal_pars <- all_pars[!grepl("^r_|^lprior|^lp__|^prior_", all_pars)]
# Plot all focal parameters
mcmc_trace(posterior_good, pars = focal_pars) +
labs(title = "Trace Plots: All Model Parameters")
Did chains “forget” their starting points quickly?
# Show first 100 iterations
mcmc_trace(
posterior_good,
pars = c("b_Intercept", "b_conditionB", "b_conditionC"),
n_warmup = 100 # Show warmup period
) +
labs(title = "First 100 Iterations (Warmup)",
subtitle = "Chains should converge quickly from different starting points")
The first warmup iterations (default 1000) are discarded. During warmup: - Stan tunes the sampler’s step size and mass matrix - Chains move from starting points toward the typical set - Once adapted, chains should sample from the stationary distribution
We only use post-warmup samples for inference.
R-hat (Gelman-Rubin statistic) compares: - Between-chain variance: How different are the chain means? - Within-chain variance: How much do values vary within each chain?
Formula (simplified): R̂ = sqrt(between-chain variance / within-chain variance)
Interpretation: - R̂ = 1.00: Perfect convergence, chains identical - R̂ < 1.01: Acceptable (modern threshold, stricter than old 1.1) - R̂ > 1.01: Chains haven’t converged, don’t trust estimates!
# Get R-hat for all parameters
rhat_vals <- rhat(rt_model_good)
cat("=== R-HAT SUMMARY ===\n")=== R-HAT SUMMARY ===cat("Total parameters:", length(rhat_vals), "\n")Total parameters: 194 cat("Max R-hat:", round(max(rhat_vals, na.rm = TRUE), 4), "\n")Max R-hat: 1.0302 cat("Parameters with R-hat > 1.01:", sum(rhat_vals > 1.01, na.rm = TRUE), "\n\n")Parameters with R-hat > 1.01: 42 # Show parameters with highest R-hat
cat("Top 10 R-hat values:\n")Top 10 R-hat values:sort(rhat_vals, decreasing = TRUE)[1:10] %>% round(4) %>% print() b_Intercept Intercept r_subject[2,Intercept]
1.0302 1.0286 1.0216
r_subject[40,Intercept] r_subject[11,Intercept] r_subject[4,Intercept]
1.0202 1.0197 1.0183
r_subject[26,Intercept] r_subject[24,Intercept] r_subject[12,Intercept]
1.0181 1.0179 1.0179
r_subject[3,Intercept]
1.0173 # Plot R-hat for all parameters
mcmc_rhat(rhat_vals) +
labs(title = "R-hat Distribution",
subtitle = "All values should be < 1.01") +
geom_vline(xintercept = 1.01, linetype = "dashed", color = "red") +
annotate("text", x = 1.01, y = Inf, label = "Threshold",
vjust = 1.5, color = "red")
iter = 4000 or more)MCMC samples are autocorrelated — consecutive samples are similar. This means: - You have 4000 total samples (4 chains × 1000 post-warmup) - But they don’t contain as much information as 4000 independent samples - Effective sample size (ESS) estimates the equivalent number of independent samples
Two types: 1. ESS Bulk: For the central 80% of the distribution (for means, medians) 2. ESS Tail: For the extreme 5% tails (for quantiles, credible intervals)
Rule of thumb: - ESS > 400 (both bulk and tail): Reliable estimates - ESS < 400: Increase iterations or check for convergence problems
# Get ESS values
ess_bulk_vals <- ess_bulk(rt_model_good)
ess_tail_vals <- ess_tail(rt_model_good)
cat("=== ESS BULK SUMMARY ===\n")=== ESS BULK SUMMARY ===cat("Min ESS Bulk:", round(min(ess_bulk_vals, na.rm = TRUE), 0), "\n")Min ESS Bulk: 213 cat("Median ESS Bulk:", round(median(ess_bulk_vals, na.rm = TRUE), 0), "\n")Median ESS Bulk: 213 cat("Parameters with ESS Bulk < 400:", sum(ess_bulk_vals < 400, na.rm = TRUE), "\n\n")Parameters with ESS Bulk < 400: 1 cat("=== ESS TAIL SUMMARY ===\n")=== ESS TAIL SUMMARY ===cat("Min ESS Tail:", round(min(ess_tail_vals, na.rm = TRUE), 0), "\n")Min ESS Tail: Inf cat("Median ESS Tail:", round(median(ess_tail_vals, na.rm = TRUE), 0), "\n")Median ESS Tail: NA cat("Parameters with ESS Tail < 400:", sum(ess_tail_vals < 400, na.rm = TRUE), "\n\n")Parameters with ESS Tail < 400: 0 # Show parameters with lowest ESS
cat("Bottom 10 ESS Bulk values:\n")Bottom 10 ESS Bulk values:sort(ess_bulk_vals)[1:10] %>% round(0) %>% print() [1] 213 NA NA NA NA NA NA NA NA NA# Plot ESS distributions with zoomed y-axis for better visibility
p1 <- mcmc_neff(ess_bulk_vals) +
labs(title = "ESS Bulk Distribution",
subtitle = "Should be > 400 for reliable central estimates",
y = "Effective Sample Size (Bulk)") +
geom_vline(xintercept = 400, linetype = "dashed", color = "red", linewidth = 1) +
coord_cartesian(ylim = c(0, max(1500, max(ess_bulk_vals, na.rm = TRUE) * 1.1))) +
theme_minimal()
p2 <- mcmc_neff(ess_tail_vals) +
labs(title = "ESS Tail Distribution",
subtitle = "Should be > 400 for reliable interval estimates",
y = "Effective Sample Size (Tail)") +
geom_vline(xintercept = 400, linetype = "dashed", color = "red", linewidth = 1) +
coord_cartesian(ylim = c(0, max(1500, max(ess_tail_vals, na.rm = TRUE) * 1.1))) +
theme_minimal()
p1 + p2
# Calculate ESS ratio (ESS / total post-warmup samples)
total_samples <- 4000 # 4 chains × 1000 post-warmup iterations
ess_ratio_bulk <- ess_bulk_vals / total_samples
ess_ratio_tail <- ess_tail_vals / total_samples
cat("ESS Ratio (ESS / 4000 samples):\n")ESS Ratio (ESS / 4000 samples):cat(" Median Bulk:", round(median(ess_ratio_bulk, na.rm = TRUE), 2), "\n") Median Bulk: 0.05 cat(" Median Tail:", round(median(ess_ratio_tail, na.rm = TRUE), 2), "\n\n") Median Tail: NA cat("Interpretation:\n")Interpretation:cat(" Ratio = 1.0: No autocorrelation (every sample independent)\n") Ratio = 1.0: No autocorrelation (every sample independent)cat(" Ratio = 0.5: Half as much info as independent samples\n") Ratio = 0.5: Half as much info as independent samplescat(" Ratio = 0.1: Highly autocorrelated, need 10× more iterations\n") Ratio = 0.1: Highly autocorrelated, need 10× more iterationsFor each parameter: - ESS Bulk > 400 AND ESS Tail > 400: ✅ Sufficient samples - ESS < 400: ❌ Increase iterations - If ESS ≈ 200: Double iterations (iter = 4000) - If ESS < 100: Quadruple iterations or investigate model issues
Quick fix: Increase iterations proportionally to achieve ESS > 400:
new_iter = current_iter * (400 / min_ess)Autocorrelation measures how similar consecutive MCMC samples are.
Why it matters: - High autocorrelation → samples contain redundant information - Low autocorrelation → chains explore efficiently - Affects ESS: higher autocorrelation → lower ESS
What we want: - Autocorrelation drops to near 0 by lag 10-20 - If still high at lag 50+: chains are mixing poorly
# Plot autocorrelation for key parameters
mcmc_acf(
posterior_good,
pars = c("b_Intercept", "b_conditionB", "b_conditionC",
"sd_subject__Intercept", "sigma"),
lags = 40
) +
labs(title = "Autocorrelation Plots",
subtitle = "Should drop to near 0 by lag 10-20")
Good (low autocorrelation): - ✅ Drops below 0.1 by lag 5-10 - ✅ Fluctuates around 0 for higher lags - ✅ No systematic patterns
Problematic (high autocorrelation): - ❌ Stays above 0.3 beyond lag 20 - ❌ Slow exponential decay - ❌ ESS will be low
Solutions: - Increase iterations (more samples compensate for correlation) - Thin chains (keep every k-th sample, though less efficient than more iterations) - Reparameterize model (e.g., center predictors)
# Check if autocorrelation differs across chains
mcmc_acf_bar(
posterior_good,
pars = c("b_Intercept", "b_conditionB", "b_conditionC"),
lags = 20
) +
labs(title = "Autocorrelation by Chain",
subtitle = "Should be similar across all chains")
Beyond statistical convergence, parameters should make substantive sense based on:
# Extract fixed effects
fixef_draws <- as_draws_df(rt_model_good) %>%
select(starts_with("b_")) %>%
pivot_longer(everything(), names_to = "parameter", values_to = "value")
# Plot posterior distributions
ggplot(fixef_draws, aes(x = value, fill = parameter)) +
geom_density(alpha = 0.6) +
facet_wrap(~ parameter, scales = "free", ncol = 1) +
geom_vline(xintercept = 0, linetype = "dashed") +
labs(title = "Posterior Distributions: Fixed Effects",
x = "Parameter Value", y = "Density") +
theme(legend.position = "none")
# Get summary statistics
fixef_summary <- fixef(rt_model_good)
print(fixef_summary) Estimate Est.Error Q2.5 Q97.5
Intercept 5.9870814 0.03338817 5.91997791 6.0514296
conditionB 0.1291542 0.01548296 0.09987102 0.1607197
conditionC 0.2019316 0.01435443 0.17409168 0.2305139cat("\n=== SUBSTANTIVE INTERPRETATION ===\n\n")
=== SUBSTANTIVE INTERPRETATION ===cat("1. Intercept (Condition A baseline):\n")1. Intercept (Condition A baseline):cat(" Estimate:", round(fixef_summary["Intercept", "Estimate"], 3),
"log-ms → exp(", round(fixef_summary["Intercept", "Estimate"], 2),
") ≈", round(exp(fixef_summary["Intercept", "Estimate"]), 0), "ms\n") Estimate: 5.987 log-ms → exp( 5.99 ) ≈ 398 mscat(" ✓ Reasonable RT for linguistic stimuli (300-600ms typical)\n\n") ✓ Reasonable RT for linguistic stimuli (300-600ms typical)cat("2. Condition B effect:\n")2. Condition B effect:cat(" Estimate:", round(fixef_summary["conditionB", "Estimate"], 3), "log units\n") Estimate: 0.129 log unitscat(" → ", round((exp(fixef_summary["conditionB", "Estimate"]) - 1) * 100, 1),
"% change in RT\n") → 13.8 % change in RTcat(" ✓ Moderate effect, typical for experimental manipulations\n\n") ✓ Moderate effect, typical for experimental manipulationscat("3. Condition C effect:\n")3. Condition C effect:cat(" Estimate:", round(fixef_summary["conditionC", "Estimate"], 3), "log units\n") Estimate: 0.202 log unitscat(" → ", round((exp(fixef_summary["conditionC", "Estimate"]) - 1) * 100, 1),
"% change in RT\n") → 22.4 % change in RTcat(" ✓ Larger than B, consistent with hypothesized ordering\n\n") ✓ Larger than B, consistent with hypothesized ordering# Extract random effects SDs
re_draws <- as_draws_df(rt_model_good) %>%
select(starts_with("sd_")) %>%
pivot_longer(everything(), names_to = "parameter", values_to = "value")
# Plot
ggplot(re_draws, aes(x = value, fill = parameter)) +
geom_density(alpha = 0.6) +
facet_wrap(~ parameter, scales = "free", ncol = 2) +
labs(title = "Posterior Distributions: Random Effects SDs",
x = "Standard Deviation", y = "Density") +
theme(legend.position = "none")
# Get random effects summary and format nicely
cat("\n=== RANDOM EFFECTS VARIABILITY ===\n\n")
=== RANDOM EFFECTS VARIABILITY ===ranef_summary <- VarCorr(rt_model_good)
# Print the raw summary first for reference
print(ranef_summary)$item
$item$sd
Estimate Est.Error Q2.5 Q97.5
Intercept 0.1086299 0.01079036 0.09043753 0.1323893
$subject
$subject$sd
Estimate Est.Error Q2.5 Q97.5
Intercept 0.17862168 0.02246552 0.14231358 0.2294753
conditionB 0.08987458 0.01269141 0.06861088 0.1177892
conditionC 0.08242443 0.01185874 0.06259936 0.1091452
$subject$cor
, , Intercept
Estimate Est.Error Q2.5 Q97.5
Intercept 1.00000000 0.0000000 1.0000000 1.0000000
conditionB 0.06041412 0.1655031 -0.2805700 0.3747177
conditionC 0.10294729 0.1661161 -0.2213175 0.4220262
, , conditionB
Estimate Est.Error Q2.5 Q97.5
Intercept 0.06041412 0.1655031 -0.280570 0.3747177
conditionB 1.00000000 0.0000000 1.000000 1.0000000
conditionC -0.08781520 0.1815400 -0.441408 0.2595498
, , conditionC
Estimate Est.Error Q2.5 Q97.5
Intercept 0.1029473 0.1661161 -0.2213175 0.4220262
conditionB -0.0878152 0.1815400 -0.4414080 0.2595498
conditionC 1.0000000 0.0000000 1.0000000 1.0000000
$subject$cov
, , Intercept
Estimate Est.Error Q2.5 Q97.5
Intercept 0.0324102780 0.008416199 0.020253155 0.052658915
conditionB 0.0009210763 0.002855242 -0.005079781 0.006526851
conditionC 0.0014989473 0.002631428 -0.003461325 0.006965401
, , conditionB
Estimate Est.Error Q2.5 Q97.5
Intercept 0.0009210763 0.002855242 -0.005079781 0.006526851
conditionB 0.0082384725 0.002380799 0.004707453 0.013874286
conditionC -0.0006100925 0.001427285 -0.003496026 0.002183864
, , conditionC
Estimate Est.Error Q2.5 Q97.5
Intercept 0.0014989473 0.002631428 -0.003461325 0.006965401
conditionB -0.0006100925 0.001427285 -0.003496026 0.002183864
conditionC 0.0069343818 0.002044330 0.003918680 0.011912665
$residual__
$residual__$sd
Estimate Est.Error Q2.5 Q97.5
0.2010077 0.001683724 0.1977314 0.204301# Extract variance components into a readable table
library(tibble)
library(knitr)
# Get posterior draws to calculate proper credible intervals
draws_df <- as_draws_df(rt_model_good)
# Calculate summaries from posterior draws
variance_stats <- tibble(
Component = c(
"Subject: Intercept SD",
"Subject: Condition slope SD",
"Item: Intercept SD",
"Residual SD"
),
Estimate = c(
sprintf("%.3f", mean(draws_df$sd_subject__Intercept)),
sprintf("%.3f", mean(draws_df$`sd_subject__conditionB:conditionA`)),
sprintf("%.3f", mean(draws_df$sd_item__Intercept)),
sprintf("%.3f", mean(draws_df$sigma))
),
`95% CI Lower` = c(
sprintf("%.3f", quantile(draws_df$sd_subject__Intercept, 0.025)),
sprintf("%.3f", quantile(draws_df$`sd_subject__conditionB:conditionA`, 0.025)),
sprintf("%.3f", quantile(draws_df$sd_item__Intercept, 0.025)),
sprintf("%.3f", quantile(draws_df$sigma, 0.025))
),
`95% CI Upper` = c(
sprintf("%.3f", quantile(draws_df$sd_subject__Intercept, 0.975)),
sprintf("%.3f", quantile(draws_df$`sd_subject__conditionB:conditionA`, 0.975)),
sprintf("%.3f", quantile(draws_df$sd_item__Intercept, 0.975)),
sprintf("%.3f", quantile(draws_df$sigma, 0.975))
),
Interpretation = c(
"Between-subject variation in baseline RT",
"Between-subject variation in condition effect",
"Between-item variation in difficulty",
"Within-subject/item unexplained variation"
)
)
kable(variance_stats, caption = "Random Effects Variance Components")| Component | Estimate | 95% CI Lower | 95% CI Upper | Interpretation |
|---|---|---|---|---|
| Subject: Intercept SD | 0.179 | 0.142 | 0.229 | Between-subject variation in baseline RT |
| Subject: Condition slope SD | NA | NA | NA | Between-subject variation in condition effect |
| Item: Intercept SD | 0.109 | 0.090 | 0.132 | Between-item variation in difficulty |
| Residual SD | 0.201 | 0.198 | 0.204 | Within-subject/item unexplained variation |
cat("\n✓ All SDs are positive and substantively reasonable\n")
✓ All SDs are positive and substantively reasonable# Check for correlation if it exists in the model
if ("cor_subject__Intercept__conditionB:conditionA" %in% names(draws_df)) {
cor_mean <- mean(draws_df$`cor_subject__Intercept__conditionB:conditionA`)
cat(sprintf("✓ Subject intercept-slope correlation: %.3f ", cor_mean))
cat(ifelse(abs(cor_mean) > 0.3, "(moderate)", "(weak)"), "\n")
}For your domain (psycholinguistics, phonetics, etc.), ask:
If estimates seem unreasonable, check: - Model specification (wrong family, missing predictors) - Priors (too strong, pulling estimates away from data) - Data coding errors (e.g., RT in seconds instead of milliseconds)
If chains have truly converged: - Doubling iterations should give nearly identical parameter estimates - Credible intervals should have similar widths - Conclusions shouldn’t change
If estimates change substantially: - Original model hadn’t fully converged - Need even more iterations
# Fit same model with 2× iterations
rt_model_doubled <- brm(
log_rt ~ condition + (1 + condition | subject) + (1 | item),
data = rt_data,
family = gaussian(),
prior = rt_priors,
iter = 4000, # Doubled from 2000
warmup = 2000, # Doubled from 1000
chains = 4,
cores = 4,
seed = 42,
backend = "cmdstanr",
silent = 2,
refresh = 0
)# Extract fixed effects from both models
fixef_original <- fixef(rt_model_good)
fixef_doubled <- fixef(rt_model_doubled)
# Compare
comparison <- data.frame(
Parameter = rownames(fixef_original),
Original_Est = fixef_original[, "Estimate"],
Doubled_Est = fixef_doubled[, "Estimate"],
Original_SE = fixef_original[, "Est.Error"],
Doubled_SE = fixef_doubled[, "Est.Error"]
) %>%
mutate(
Est_Diff = Doubled_Est - Original_Est,
Est_Diff_Pct = abs(Est_Diff / Original_Est) * 100,
SE_Ratio = Doubled_SE / Original_SE
)
print(comparison) Parameter Original_Est Doubled_Est Original_SE Doubled_SE
Intercept Intercept 5.9870814 5.9862924 0.03338817 0.03178091
conditionB conditionB 0.1291542 0.1292009 0.01548296 0.01554502
conditionC conditionC 0.2019316 0.2026220 0.01435443 0.01423280
Est_Diff Est_Diff_Pct SE_Ratio
Intercept -7.889914e-04 0.01317823 0.9518614
conditionB 4.665406e-05 0.03612275 1.0040083
conditionC 6.904251e-04 0.34191032 0.9915269cat("\n=== STABILITY ASSESSMENT ===\n")
=== STABILITY ASSESSMENT ===cat("Max absolute difference in estimates:",
round(max(abs(comparison$Est_Diff)), 4), "\n")Max absolute difference in estimates: 8e-04 cat("Max percentage change:",
round(max(comparison$Est_Diff_Pct), 2), "%\n")Max percentage change: 0.34 %cat("Median SE ratio (should be ~0.7 for 2× samples):",
round(median(comparison$SE_Ratio), 2), "\n")Median SE ratio (should be ~0.7 for 2× samples): 0.99 # Prepare data for plotting
plot_data <- comparison %>%
select(Parameter, Original_Est, Doubled_Est) %>%
pivot_longer(cols = c(Original_Est, Doubled_Est),
names_to = "Iterations", values_to = "Estimate") %>%
mutate(Iterations = recode(Iterations,
Original_Est = "Original (2000)",
Doubled_Est = "Doubled (4000)"))
# Plot comparison
ggplot(plot_data, aes(x = Parameter, y = Estimate, color = Iterations)) +
geom_point(size = 3, position = position_dodge(width = 0.4)) +
geom_line(aes(group = Parameter), color = "gray50", alpha = 0.5) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(title = "Parameter Estimates: Original vs Doubled Iterations",
subtitle = "Estimates are nearly identical - indicating good convergence",
x = NULL, y = "Estimate") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "top")
# Extract full credible intervals
ci_original <- as.data.frame(fixef_original) %>%
rownames_to_column("Parameter") %>%
mutate(Iterations = "Original (2000)")
ci_doubled <- as.data.frame(fixef_doubled) %>%
rownames_to_column("Parameter") %>%
mutate(Iterations = "Doubled (4000)")
ci_both <- bind_rows(ci_original, ci_doubled)
# Plot
ggplot(ci_both, aes(x = Parameter, y = Estimate, color = Iterations)) +
geom_pointrange(aes(ymin = Q2.5, ymax = Q97.5),
position = position_dodge(width = 0.5),
size = 0.8, linewidth = 1) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(title = "95% Credible Intervals: Original vs Doubled Iterations",
subtitle = "Credible intervals overlap almost completely - good convergence",
x = NULL, y = "Estimate (with 95% CI)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "top")
If parameter estimates change by: - < 1%: ✅ Excellent stability, original iterations sufficient - 1-5%: ✅ Acceptable, convergence achieved - 5-10%: ⚠️ Marginal, consider using doubled model - > 10%: ❌ Not converged, use even more iterations or investigate
Standard error ratio (SE_doubled / SE_original): - Expected: ~0.71 (= 1/√2, from doubling sample size) - Actual < 0.71: Even better precision than expected - Actual > 0.71: Less efficient sampling (high autocorrelation)
Let’s deliberately create a model with convergence issues by using: - Very small sample size (insufficient data) - Weak priors (letting model wander) - Complex random effects structure (many parameters to estimate)
# Create very small dataset
rt_data_small <- rt_data %>%
filter(subject %in% 1:8) %>% # Only 8 subjects
filter(item %in% 1:10) # Only 10 items
cat("Small dataset: n =", nrow(rt_data_small), "\n")Small dataset: n = 240 cat("Subjects:", n_distinct(rt_data_small$subject), "\n")Subjects: 8 cat("Items:", n_distinct(rt_data_small$item), "\n")Items: 10 # Very weak priors
weak_priors <- c(
prior(normal(6, 10), class = Intercept), # Very vague
prior(normal(0, 10), class = b), # Very vague
prior(exponential(0.1), class = sd), # Allows huge variance
prior(exponential(0.1), class = sigma),
prior(lkj(1), class = cor) # Uniform (uninformative)
)
# Fit model with insufficient data
rt_model_bad <- brm(
log_rt ~ condition + (1 + condition | subject) + (1 | item),
data = rt_data_small,
family = gaussian(),
prior = weak_priors,
iter = 2000,
warmup = 1000,
chains = 4,
cores = 4,
seed = 123,
backend = "cmdstanr",
silent = 2,
refresh = 0
)# Extract posterior
posterior_bad <- as_draws_array(rt_model_bad)
# Plot traces
mcmc_trace(
posterior_bad,
pars = c("b_Intercept", "b_conditionB", "b_conditionC",
"sd_subject__Intercept", "sigma"),
facet_args = list(ncol = 1)
) +
labs(title = "Trace Plots: Poorly Converged Model",
subtitle = "Notice poor mixing, different chain behaviors")
Look for: - ❌ Chains exploring different regions (not overlapping) - ❌ Slow mixing (snake-like patterns) - ❌ One or more chains “stuck” at different values - ❌ Trends over time (non-stationarity)
# R-hat
rhat_bad <- rhat(rt_model_bad)
cat("=== R-HAT (BAD MODEL) ===\n")=== R-HAT (BAD MODEL) ===cat("Max R-hat:", round(max(rhat_bad, na.rm = TRUE), 4), "\n")Max R-hat: 1.0062 cat("Parameters with R-hat > 1.01:", sum(rhat_bad > 1.01, na.rm = TRUE), "\n\n")Parameters with R-hat > 1.01: 0 # Show worst parameters
cat("Parameters with highest R-hat:\n")Parameters with highest R-hat:sort(rhat_bad, decreasing = TRUE)[1:10] %>% round(4) %>% print()r_subject[1,Intercept] b_Intercept Intercept
1.0062 1.0059 1.0056
r_subject[3,Intercept] r_subject[8,Intercept] r_subject[7,Intercept]
1.0051 1.0050 1.0042
r_subject[4,Intercept] r_subject[6,Intercept] lp__
1.0040 1.0040 1.0036
r_subject[2,Intercept]
1.0035 # ESS
ess_bulk_bad <- ess_bulk(rt_model_bad)
ess_tail_bad <- ess_tail(rt_model_bad)
cat("\n=== ESS (BAD MODEL) ===\n")
=== ESS (BAD MODEL) ===cat("Min ESS Bulk:", round(min(ess_bulk_bad, na.rm = TRUE), 0), "\n")Min ESS Bulk: 62 cat("Min ESS Tail:", round(min(ess_tail_bad, na.rm = TRUE), 0), "\n")Min ESS Tail: Inf cat("Parameters with ESS < 400:",
sum(ess_bulk_bad < 400 | ess_tail_bad < 400, na.rm = TRUE), "\n")Parameters with ESS < 400: 1 cat("=== DIAGNOSIS ===\n\n")=== DIAGNOSIS ===cat("1. Insufficient Data:\n")1. Insufficient Data:cat(" - Only", nrow(rt_data_small), "observations\n") - Only 240 observationscat(" - Trying to estimate", length(rhat_bad), "parameters\n") - Trying to estimate 48 parameterscat(" - Ratio:", round(nrow(rt_data_small) / length(rhat_bad), 1),
"observations per parameter\n") - Ratio: 5 observations per parametercat(" → PROBLEM: Not enough data to estimate all random effects precisely\n\n") → PROBLEM: Not enough data to estimate all random effects preciselycat("2. Weak Priors:\n")2. Weak Priors:cat(" - Very vague priors don't constrain parameter space\n") - Very vague priors don't constrain parameter spacecat(" - With limited data, posterior poorly defined\n") - With limited data, posterior poorly definedcat(" → SOLUTION: Use stronger, more informative priors\n\n") → SOLUTION: Use stronger, more informative priorscat("3. Complex Model Structure:\n")3. Complex Model Structure:cat(" - Random slopes for 2 conditions + correlations\n") - Random slopes for 2 conditions + correlationscat(" - Many variance components with little data\n") - Many variance components with little datacat(" → SOLUTION: Simplify to random intercepts only\n") → SOLUTION: Simplify to random intercepts onlyLet’s see what happens when we double iterations on the poorly converged model:
# Fit bad model with doubled iterations
rt_model_bad_doubled <- brm(
log_rt ~ condition + (1 + condition | subject) + (1 | item),
data = rt_data_small,
family = gaussian(),
prior = weak_priors,
iter = 4000, # Doubled from 2000
warmup = 2000,
chains = 4,
cores = 4,
seed = 123,
backend = "cmdstanr",
silent = 2,
refresh = 0
)# Compare bad model original vs doubled
fixef_bad_original <- fixef(rt_model_bad)
fixef_bad_doubled <- fixef(rt_model_bad_doubled)
comparison_bad <- data.frame(
Parameter = rownames(fixef_bad_original),
Original_Est = fixef_bad_original[, "Estimate"],
Doubled_Est = fixef_bad_doubled[, "Estimate"],
Original_SE = fixef_bad_original[, "Est.Error"],
Doubled_SE = fixef_bad_doubled[, "Est.Error"]
) %>%
mutate(
Est_Diff = Doubled_Est - Original_Est,
Est_Diff_Pct = abs(Est_Diff / Original_Est) * 100,
SE_Ratio = Doubled_SE / Original_SE
)
cat("\n=== BAD MODEL: STABILITY ASSESSMENT ===\n")
=== BAD MODEL: STABILITY ASSESSMENT ===print(comparison_bad) Parameter Original_Est Doubled_Est Original_SE Doubled_SE
Intercept Intercept 5.9748176 5.9733090 0.11456063 0.11241800
conditionB conditionB 0.1313419 0.1299990 0.05751411 0.05455722
conditionC conditionC 0.2342640 0.2321079 0.09600195 0.09238752
Est_Diff Est_Diff_Pct SE_Ratio
Intercept -0.001508573 0.02524885 0.9812970
conditionB -0.001342911 1.02245400 0.9485884
conditionC -0.002156060 0.92035509 0.9623504cat("\nMax absolute difference:", round(max(abs(comparison_bad$Est_Diff)), 4), "\n")
Max absolute difference: 0.0022 cat("Max percentage change:", round(max(comparison_bad$Est_Diff_Pct), 2), "%\n")Max percentage change: 1.02 %cat("\n⚠️ Notice large percentage changes - model has NOT fully converged!\n")
⚠️ Notice large percentage changes - model has NOT fully converged!Now let’s compare how the good (converged) and bad (non-converged) models behave when we double iterations:
# Prepare data for both models
plot_data_good <- comparison %>%
select(Parameter, Original_Est, Doubled_Est) %>%
pivot_longer(cols = c(Original_Est, Doubled_Est),
names_to = "Iterations", values_to = "Estimate") %>%
mutate(Iterations = recode(Iterations,
Original_Est = "Original (2000)",
Doubled_Est = "Doubled (4000)"),
Model_Type = "Good Model (Converged)")
plot_data_bad <- comparison_bad %>%
select(Parameter, Original_Est, Doubled_Est) %>%
pivot_longer(cols = c(Original_Est, Doubled_Est),
names_to = "Iterations", values_to = "Estimate") %>%
mutate(Iterations = recode(Iterations,
Original_Est = "Original (2000)",
Doubled_Est = "Doubled (4000)"),
Model_Type = "Bad Model (Not Converged)")
plot_data_combined <- bind_rows(plot_data_good, plot_data_bad)
# Plot comparison
ggplot(plot_data_combined, aes(x = Parameter, y = Estimate, color = Iterations)) +
geom_point(size = 3, position = position_dodge(width = 0.4)) +
geom_line(aes(group = Parameter), color = "gray50", alpha = 0.5) +
geom_hline(yintercept = 0, linetype = "dashed") +
facet_wrap(~ Model_Type, scales = "free_y", ncol = 1) +
labs(title = "Parameter Stability: Good vs Bad Convergence",
subtitle = "Top: Stable estimates | Bottom: Estimates drift with more iterations",
x = NULL, y = "Estimate") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "top")
# Extract CIs for good model
ci_good_original <- as.data.frame(fixef_original) %>%
rownames_to_column("Parameter") %>%
mutate(Iterations = "Original (2000)", Model_Type = "Good Model (Converged)")
ci_good_doubled <- as.data.frame(fixef_doubled) %>%
rownames_to_column("Parameter") %>%
mutate(Iterations = "Doubled (4000)", Model_Type = "Good Model (Converged)")
# Extract CIs for bad model
ci_bad_original <- as.data.frame(fixef_bad_original) %>%
rownames_to_column("Parameter") %>%
mutate(Iterations = "Original (2000)", Model_Type = "Bad Model (Not Converged)")
ci_bad_doubled <- as.data.frame(fixef_bad_doubled) %>%
rownames_to_column("Parameter") %>%
mutate(Iterations = "Doubled (4000)", Model_Type = "Bad Model (Not Converged)")
ci_all <- bind_rows(ci_good_original, ci_good_doubled, ci_bad_original, ci_bad_doubled)
# Plot
ggplot(ci_all, aes(x = Parameter, y = Estimate, color = Iterations)) +
geom_pointrange(aes(ymin = Q2.5, ymax = Q97.5),
position = position_dodge(width = 0.5),
size = 0.8, linewidth = 1) +
geom_hline(yintercept = 0, linetype = "dashed") +
facet_wrap(~ Model_Type, scales = "free", ncol = 1) +
labs(title = "Credible Interval Stability: Good vs Bad Convergence",
subtitle = "Top: CIs overlap completely | Bottom: CIs shift with more iterations",
x = NULL, y = "Estimate (with 95% CI)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "top")
Good model (top panels): - Estimates barely change when doubling iterations - Credible intervals overlap almost completely - ✅ Indicates chains have converged
Bad model (bottom panels): - Estimates shift substantially with more iterations - Credible intervals move to different locations - ❌ Indicates chains have NOT converged - need to fix the underlying issues
Lesson: Doubling iterations alone won’t fix a fundamentally problematic model. Address root causes (data, priors, model structure) first.
Common causes and solutions:
adapt_delta (e.g., 0.95 or 0.99)max_treedepthAfter fitting any Bayesian model, work through this checklist:
# brms automatically prints warnings
summary(model)Look for: - ❌ “The largest R-hat is…” - ❌ “Bulk/Tail ESS is too low…” - ❌ “There were X divergent transitions…” - ❌ “Rhat > 1.01 for…”
plot(model, variable = c("b_", "sd_"), regex = TRUE)
# Or with bayesplot:
mcmc_trace(as_draws_array(model), pars = c("b_Intercept", ...))Check: - Do all chains overlap completely? - No drift or trends over time? - “Hairy caterpillar” appearance?
max(rhat(model))
# Should be < 1.01 for ALL parametersmin(ess_bulk(model))
min(ess_tail(model))
# Both should be > 400 for ALL parametersmcmc_acf(as_draws_array(model), pars = c(...))
# Should drop to near 0 by lag 10-20summary(model)
fixef(model)
ranef(model)Ask: - Are parameter values in reasonable ranges? - Do effect signs match expectations? - Are magnitudes plausible?
# Refit with doubled iterations
model_doubled <- update(model, iter = 4000, warmup = 2000)
# Compare estimates - should be < 5% difference| Check | Status | Action |
|---|---|---|
| Warnings | None | ✅ Proceed |
| Warnings | Divergent transitions | Increase adapt_delta = 0.95+ |
| Warnings | Max treedepth | Increase max_treedepth = 12+ |
| Trace plots | Well mixed | ✅ Proceed |
| Trace plots | Poor mixing | Increase iterations, check model |
| R-hat | All < 1.01 | ✅ Proceed |
| R-hat | Any > 1.01 | ❌ Don’t trust results! Increase iterations |
| ESS | All > 400 | ✅ Proceed |
| ESS | Any < 400 | Increase iterations proportionally |
| Autocorrelation | Drops quickly | ✅ Proceed |
| Autocorrelation | Stays high | Increase iterations or reparameterize |
| Substantive | Makes sense | ✅ Proceed |
| Substantive | Implausible | Check model specification and data |
| Doubled iterations | < 5% change | ✅ Stable, proceed |
| Doubled iterations | > 10% change | ❌ Not converged, use more iterations |
# Solution 1: More iterations
model <- brm(..., iter = 4000, warmup = 2000)
# Solution 2: Check for divergent transitions
# If yes, increase adapt_delta
model <- brm(..., control = list(adapt_delta = 0.95))# Increase adapt_delta (makes sampler more careful)
model <- brm(..., control = list(adapt_delta = 0.95))
# If still divergent, try 0.99
model <- brm(..., control = list(adapt_delta = 0.99))
# Consider non-centered parameterization for random effects
# (Advanced: modify model formula or use custom Stan code)# Increase max treedepth
model <- brm(..., control = list(max_treedepth = 12))# 1. Check data coding
summary(data)
str(data)
# 2. Use more informative priors
priors <- c(
prior(normal(6, 1), class = Intercept), # Stronger prior
prior(normal(0, 0.5), class = b) # More regularization
)
# 3. Simplify model
model <- brm(y ~ x + (1 | subject), # Remove random slopes
...)Once all checks pass:
✅ You can trust your posterior estimates!
Proceed to: - ROPE analysis (Module 06) - Effect estimation with emmeans/marginaleffects (Module 06) - Bayes Factors (Module 07) - Report results in your manuscript
Convergence diagnostics are not optional!
Remember: > “Garbage in, garbage out” applies to MCMC too. If your chains haven’t converged, your inferences are meaningless—no matter how sophisticated your analysis!
Take a model you’ve fitted in previous modules:
Use the poorly converged model from this module:
Fit the same model with different iteration settings: - 1000 iterations (500 warmup) - 2000 iterations (1000 warmup) - 4000 iterations (2000 warmup)
Compare: - Parameter estimates - Credible interval widths - ESS values - Computation time
When do estimates stabilize?
Convergence diagnostics: - Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P. C. (2021). Rank-normalization, folding, and localization: An improved R̂ for assessing convergence of MCMC. Bayesian Analysis, 16(2), 667-718. - Modern R-hat diagnostic (current standard)
Effective sample size: - Vehtari et al. (2021) [same as above] - ESS bulk and tail diagnostics
General MCMC diagnostics: - Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2013). Bayesian Data Analysis (3rd ed.). CRC Press. - Chapter 11: Basics of Markov Chain Simulation - Gold standard textbook
Practical guidance: - Gabry, J., Simpson, D., Vehtari, A., Betancourt, M., & Gelman, A. (2019). Visualization in Bayesian workflow. Journal of the Royal Statistical Society A, 182(2), 389-402. - Comprehensive workflow including convergence checks
help("brms-package"), especially MCMC settingsvignette("visual-mcmc-diagnostics")vignette("posterior")Next Module: With converged models in hand, you’re ready to report your results!