Code
#install.packages("remotes")
#remotes::install_github("DevPsyLab/petersenlab")Bayesian Analysis
Adapted from brms workshop by Paul-Christian Bürkner
1 Preamble
1.1 Install Libraries
1.2 Load Libraries
1.3 Simulate Data
Code
set.seed(52242)
sampleSize <- 1000
id <- rep(1:100, each = 10)
X <- rnorm(sampleSize)
M <- 0.5*X + rnorm(sampleSize)
Y <- 0.7*M + rnorm(sampleSize)
nonlinear_interaction <- 1.5 * sin(X) * M
Y2 <- 0.7 * M + nonlinear_interaction + rnorm(sampleSize, sd = 1)
X[sample(1:length(X), size = 10)] <- NA
M[sample(1:length(M), size = 10)] <- NA
Y[sample(1:length(Y), size = 10)] <- NA
Y2[sample(1:length(Y2), size = 10)] <- NA
mydata <- data.frame(
id = id,
X = X,
Y = Y,
Y2 = Y2,
M = M)1.4 Load Data
Code
data("sleepstudy", package = "lme4")1.5 Prepare Data
Code
conditions <- make_conditions(sleepstudy, "Subject")2 brms
2.1 Post-Processing Methods
Code
methods(class = "brmsfit") [1] add_criterion add_ic as_draws_array
[4] as_draws_df as_draws_list as_draws_matrix
[7] as_draws_rvars as_draws as.array
[10] as.data.frame as.matrix as.mcmc
[13] autocor bayes_factor bayes_R2
[16] bayesfactor_models bayesfactor_parameters bayesfactor_restricted
[19] bci bridge_sampler check_model
[22] check_predictions check_prior check_priors
[25] ci clean_parameters coef
[28] conditional_effects conditional_smooths control_params
[31] default_prior describe_posterior describe_prior
[34] diagnostic_draws diagnostic_posterior effective_sample
[37] equivalence_test estimate_density estimate_grouplevel
[40] eta_squared_posterior eti expose_functions
[43] family find_algorithm find_auxiliary
[46] find_formula find_parameters find_predictors
[49] find_response find_weights fitted
[52] fixef formula get_data
[55] get_dispersion get_mixed_info get_modelmatrix
[58] get_parameters get_predicted get_priors
[61] get_residuals get_simulated get_varcov
[64] get_variance get_weights getCall
[67] hdi hypothesis inits
[70] kfold link_function link_inverse
[73] log_lik log_posterior logLik
[76] loo_compare loo_epred loo_linpred
[79] loo_model_weights loo_moment_match loo_predict
[82] loo_predictive_interval loo_R2 loo_subsample
[85] loo LOO map_estimate
[88] marginal_effects marginal_smooths mcmc_plot
[91] mcse mediation model_info
[94] model_parameters model_performance model_to_priors
[97] model_weights model.frame n_parameters
[100] nchains ndraws neff_ratio
[103] ngrps niterations nobs
[106] nsamples nuts_params nvariables
[109] p_direction p_map p_rope
[112] p_significance p_value pairs
[115] parnames performance_logloss plot
[118] point_estimate post_prob posterior_average
[121] posterior_epred posterior_interval posterior_linpred
[124] posterior_predict posterior_samples posterior_smooths
[127] posterior_summary pp_average pp_check
[130] pp_mixture predict predictive_error
[133] predictive_interval prepare_predictions print
[136] prior_draws prior_summary psis
[139] r2_loo_posterior r2_posterior r2
[142] ranef reloo report_effectsize
[145] report_info report_intercept report_model
[148] report_parameters report_performance report_priors
[151] report_random report_statistics report_table
[154] report_text report residuals
[157] restructure rhat rope
[160] sexit_thresholds si simulate_prior
[163] spi stancode standard_error
[166] standardize standata stanplot
[169] summary unupdate update
[172] VarCorr variables vcov
[175] waic WAIC weighted_posteriors
see '?methods' for accessing help and source code3 Multilevel Models
3.1 Fit the Models
3.1.1 Complete Pooling
Code
fit_sleep1 <- brm(
Reaction ~ 1 + Days,
data = sleepstudy,
seed = 52242)3.1.2 Random Intercepts
Code
fit_sleep2 <- brm(
Reaction ~ 1 + Days + (1 | Subject),
data = sleepstudy,
seed = 52242)3.1.3 Random Intercepts and Slopes
Code
fit_sleep3 <- brm(
Reaction ~ 1 + Days + (1 + Days | Subject),
data = sleepstudy,
seed = 52242
)3.2 Summarize Results
For convergence, Rhat values should not be above 1.00.
Code
summary(fit_sleep3) Family: gaussian
Links: mu = identity
Formula: Reaction ~ 1 + Days + (1 + 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.76 6.70 15.60 41.96 1.00 1836 2290
sd(Days) 6.58 1.55 4.17 10.23 1.00 1265 1938
cor(Intercept,Days) 0.08 0.30 -0.49 0.66 1.00 1041 1736
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 251.38 7.39 236.80 266.10 1.00 1912 2086
Days 10.43 1.76 7.03 13.97 1.00 1294 1865
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 25.96 1.55 23.13 29.23 1.00 3535 2807
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).3.3 Model Priors
Code
prior_summary(fit_sleep3)3.4 Model Parameters
Code
variables(fit_sleep3) [1] "b_Intercept" "b_Days"
[3] "sd_Subject__Intercept" "sd_Subject__Days"
[5] "cor_Subject__Intercept__Days" "sigma"
[7] "Intercept" "r_Subject[308,Intercept]"
[9] "r_Subject[309,Intercept]" "r_Subject[310,Intercept]"
[11] "r_Subject[330,Intercept]" "r_Subject[331,Intercept]"
[13] "r_Subject[332,Intercept]" "r_Subject[333,Intercept]"
[15] "r_Subject[334,Intercept]" "r_Subject[335,Intercept]"
[17] "r_Subject[337,Intercept]" "r_Subject[349,Intercept]"
[19] "r_Subject[350,Intercept]" "r_Subject[351,Intercept]"
[21] "r_Subject[352,Intercept]" "r_Subject[369,Intercept]"
[23] "r_Subject[370,Intercept]" "r_Subject[371,Intercept]"
[25] "r_Subject[372,Intercept]" "r_Subject[308,Days]"
[27] "r_Subject[309,Days]" "r_Subject[310,Days]"
[29] "r_Subject[330,Days]" "r_Subject[331,Days]"
[31] "r_Subject[332,Days]" "r_Subject[333,Days]"
[33] "r_Subject[334,Days]" "r_Subject[335,Days]"
[35] "r_Subject[337,Days]" "r_Subject[349,Days]"
[37] "r_Subject[350,Days]" "r_Subject[351,Days]"
[39] "r_Subject[352,Days]" "r_Subject[369,Days]"
[41] "r_Subject[370,Days]" "r_Subject[371,Days]"
[43] "r_Subject[372,Days]" "lprior"
[45] "lp__" 3.5 Model Coefficients
Code
coef(fit_sleep3)$Subject
, , Intercept
Estimate Est.Error Q2.5 Q97.5
308 253.8036 13.42340 227.2528 279.4458
309 211.5798 13.53907 184.5701 237.5476
310 212.9776 13.44786 186.3989 238.6100
330 274.5239 13.43628 248.1815 301.9403
331 273.0398 13.23918 248.1694 299.4477
332 260.2840 12.23798 236.3152 284.9809
333 267.8093 12.62690 244.2342 293.6989
334 244.3311 12.62822 219.1096 269.0711
335 250.7611 13.16495 224.1951 276.3869
337 286.1682 13.38914 259.9519 312.2581
349 226.5570 12.82480 200.4294 251.4542
350 238.6786 13.13686 212.3307 263.6997
351 255.7790 12.53802 230.2905 280.3204
352 272.2132 12.54778 247.9377 297.7906
369 254.4293 12.50311 230.2014 278.3378
370 226.7179 13.30644 200.0374 251.3793
371 252.3647 12.07777 229.1029 275.8134
372 263.5458 12.38755 239.3714 287.7161
, , Days
Estimate Est.Error Q2.5 Q97.5
308 19.6319154 2.542400 14.91127989 24.921229
309 1.7458810 2.561539 -3.16997162 6.728358
310 4.9887767 2.536991 0.03313305 10.034349
330 5.7537881 2.555568 0.62257977 10.503648
331 7.5170261 2.469967 2.42232385 12.196072
332 10.2271380 2.325028 5.69008473 14.750818
333 10.2869556 2.337880 5.51782984 14.730549
334 11.5339766 2.405073 6.85791261 16.290102
335 -0.2252846 2.582433 -5.24670583 4.873086
337 19.1100826 2.537751 14.24044379 24.060945
349 11.5705306 2.403886 7.01645714 16.487190
350 17.0331885 2.543700 12.24508630 22.034614
351 7.4875668 2.395344 2.89592486 12.176265
352 14.0213129 2.421787 9.20693205 18.780621
369 11.3270352 2.379922 6.54155107 15.929727
370 15.1243321 2.534860 10.24501333 20.291155
371 9.4030963 2.306262 4.87040862 14.017797
372 11.7617904 2.361364 7.07438198 16.4750363.6 Standardized Coefficients
Code
print(effectsize::standardize_parameters(
fit_sleep3,
method = "posthoc"),
digits = 2)# Standardization method: posthoc
Component | Parameter | Std. Median | 95% CI
------------------------------------------------------
conditional | (Intercept) | 0.00 | [0.00, 0.00]
conditional | Days | 0.18 | [0.12, 0.25]
sigma | sigma | 0.46 | [0.41, 0.52]3.7 Model Fit
Code
performance::model_performance(fit_sleep3)3.8 R-Squared
Code
performance::r2(fit_sleep3)# Bayesian R2 with Compatibility Interval
Conditional R2: 0.793 (95% CI [0.755, 0.822])
Marginal R2: 0.284 (95% CI [0.145, 0.402])3.9 Plots
3.9.1 Trace Plots
3.9.2 Visualize Predictions
3.9.2.1 Sample-Level
Code
plot(conditional_effects(fit_sleep1), points = TRUE)
3.9.2.2 Person-Level
Code
# re_formula = NULL ensures that group-level effects are included
ce2 <- conditional_effects(
fit_sleep3,
conditions = conditions,
re_formula = NULL)
plot(ce2, ncol = 6, points = TRUE)
3.9.3 Check Model Fit
3.9.3.1 Posterior Predictive Check
Evaluate how closely the posterior predictions match the observed values. If they do not match the general pattern of the observed values, a different response distribution may be necessary.
3.10 Fitted Values
Code
fitted(fit_sleep3) Estimate Est.Error Q2.5 Q97.5
[1,] 253.8036 13.423405 227.2528 279.4458
[2,] 273.4355 11.520427 251.1066 295.5886
[3,] 293.0674 9.908596 273.5416 312.1998
[4,] 312.6994 8.750308 295.6254 329.6635
[5,] 332.3313 8.239121 316.1146 348.5817
[6,] 351.9632 8.492698 335.3675 368.36323.11 Residuals
Code
residuals(fit_sleep3) Estimate Est.Error Q2.5 Q97.5
[1,] -4.372661 29.28893 -62.62331 52.58033
[2,] -14.630370 28.37156 -69.95154 39.94352
[3,] -42.373194 27.64747 -95.88675 11.66276
[4,] 8.848878 27.81129 -46.68467 62.51537
[5,] 24.240488 26.95768 -26.35339 78.13108
[6,] 62.036557 27.19291 7.32166 114.397613.12 Compare Models
elpd values: higher is better looic values: lower is better
elpd_diff values that are greater than ~2 standard errors of the elpd_diff values indicate a significantly better model (i.e., if elpd_diff value is greater than 2 times the se_diff value).
Code
loo(fit_sleep1, fit_sleep2, fit_sleep3)Output of model 'fit_sleep1':
Computed from 4000 by 180 log-likelihood matrix.
Estimate SE
elpd_loo -953.1 10.5
p_loo 3.0 0.5
looic 1906.3 21.0
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.9, 1.2]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.
Output of model 'fit_sleep2':
Computed from 4000 by 180 log-likelihood matrix.
Estimate SE
elpd_loo -884.7 14.3
p_loo 19.2 3.3
looic 1769.4 28.7
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.8]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.
Output of model 'fit_sleep3':
Computed from 4000 by 180 log-likelihood matrix.
Estimate SE
elpd_loo -861.2 22.1
p_loo 34.0 8.2
looic 1722.4 44.1
------
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% 49
(0.7, 1] (bad) 3 1.7% <NA>
(1, Inf) (very bad) 0 0.0% <NA>
See help('pareto-k-diagnostic') for details.
Model comparisons:
model elpd_diff se_diff p_worse diag_diff diag_elpd
fit_sleep3 0.0 0.0 NA 3 k_psis > 0.7
fit_sleep2 -23.5 11.4 0.98
fit_sleep1 -91.9 20.7 1.00 Output of model 'fit_sleep1':
Computed from 4000 by 180 log-likelihood matrix.
Estimate SE
elpd_loo -953.1 10.5
p_loo 3.0 0.5
looic 1906.3 21.0
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.9, 1.2]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.
Output of model 'fit_sleep2':
Computed from 4000 by 180 log-likelihood matrix.
Estimate SE
elpd_loo -884.7 14.3
p_loo 19.2 3.3
looic 1769.4 28.7
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.8]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.
Output of model 'fit_sleep3':
Computed from 4000 by 180 log-likelihood matrix.
Estimate SE
elpd_loo -861.2 22.1
p_loo 34.0 8.2
looic 1722.4 44.1
------
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% 49
(0.7, 1] (bad) 3 1.7% <NA>
(1, Inf) (very bad) 0 0.0% <NA>
See help('pareto-k-diagnostic') for details.
Model comparisons:
model elpd_diff se_diff p_worse diag_diff diag_elpd elpd_loo
fit_sleep3 0.0 0.0 NA 3 k_psis > 0.7 -861.2
fit_sleep2 -23.5 11.4 0.98 -884.7
fit_sleep1 -91.9 20.7 1.00 -953.1
se_elpd_loo p_loo se_p_loo looic se_looic
22.1 34.0 8.2 1722.4 44.1
14.3 19.2 3.3 1769.4 28.7
10.5 3.0 0.5 1906.3 21.03.12.1 Compute Model Weights
Code
model_weights(fit_sleep1, fit_sleep2, fit_sleep3, weights = "loo") fit_sleep1 fit_sleep2 fit_sleep3
1.192945e-40 6.120337e-11 1.000000e+00 Code
round(model_weights(fit_sleep1, fit_sleep2, fit_sleep3, weights = "loo"))fit_sleep1 fit_sleep2 fit_sleep3
0 0 1 4 Multilevel Mediation
The syntax below estimates random intercepts (which allows each participant to have a different intercept) to account for nested data within the same participant.
Code
bayesianMediationSyntax <-
bf(M ~ X + (1 |i| id)) +
bf(Y ~ X + M + (1 |i| id)) +
set_rescor(FALSE) # don't add a residual correlation between M and YCode
bayesianMediationModel <- brm(
bayesianMediationSyntax,
data = mydata,
seed = 52242
)Code
summary(bayesianMediationModel) Family: MV(gaussian, gaussian)
Links: mu = identity
mu = identity
Formula: M ~ X + (1 | i | id)
Y ~ X + M + (1 | i | id)
Data: mydata (Number of observations: 970)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Multilevel Hyperparameters:
~id (Number of levels: 100)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(M_Intercept) 0.07 0.05 0.00 0.18 1.00 1566
sd(Y_Intercept) 0.10 0.06 0.01 0.22 1.00 1149
cor(M_Intercept,Y_Intercept) -0.03 0.56 -0.95 0.94 1.00 2070
Tail_ESS
sd(M_Intercept) 2217
sd(Y_Intercept) 1670
cor(M_Intercept,Y_Intercept) 2596
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
M_Intercept 0.00 0.03 -0.06 0.07 1.00 6911 2510
Y_Intercept 0.06 0.04 -0.01 0.13 1.00 5359 2443
M_X 0.51 0.03 0.44 0.57 1.00 7437 2650
Y_X 0.03 0.04 -0.04 0.11 1.00 6972 3170
Y_M 0.69 0.03 0.63 0.76 1.00 6418 2803
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_M 1.01 0.02 0.97 1.06 1.00 7560 2803
sigma_Y 1.01 0.02 0.96 1.06 1.00 5275 2747
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).Code
#print(effectsize::standardize_parameters( # errors out perhaps due to joint model
# bayesianMediationModel,
# method = "posthoc"),
# digits = 2)
#performance::model_performance(bayesianMediationModel) # errors out perhaps due to joint model
performance::r2(bayesianMediationModel) M.R2_Bayes M.R2_Bayes_marginal Y.R2_Bayes Y.R2_Bayes_marginal
0.2032413 0.1978159 0.3915377 0.3842677
attr(,"class")
[1] "r2_bayes_mv"
attr(,"SE")
M.R2_Bayes M.R2_Bayes_marginal Y.R2_Bayes Y.R2_Bayes_marginal
0.02066753 0.02032979 0.02079362 0.01954297
attr(,"CI")
M.R2_Bayes.CI M.R2_Bayes.CI_low
0.9500000 0.1604359
M.R2_Bayes.CI_high M.R2_Bayes_marginal.CI
0.2438259 0.9500000
M.R2_Bayes_marginal.CI_low M.R2_Bayes_marginal.CI_high
0.1580773 0.2385614
Y.R2_Bayes.CI Y.R2_Bayes.CI_low
0.9500000 0.3491706
Y.R2_Bayes.CI_high Y.R2_Bayes_marginal.CI
0.4299376 0.9500000
Y.R2_Bayes_marginal.CI_low Y.R2_Bayes_marginal.CI_high
0.3446632 0.4220177
attr(,"ci_method")
[1] "HDI"
attr(,"robust")
[1] TRUECode
hypothesis(
bayesianMediationModel,
"b_M_X * b_Y_M = 0", # indirect effect = a path * b path
class = NULL,
seed = 52242
)Hypothesis Tests for class :
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
1 (b_M_X*b_Y_M) = 0 0.35 0.03 0.3 0.41 NA NA
Star
1 *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.Code
mediation(bayesianMediationModel)5 Multilevel Moderation
Code
bayesianModerationModel <- brm(
bayesianModerationSyntax,
data = mydata,
seed = 52242
)Code
summary(bayesianModerationModel) Family: gaussian
Links: mu = identity
Formula: Y2 ~ s(X) + s(M) + t2(X, M) + (1 | id)
Data: mydata (Number of observations: 970)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Smoothing Spline Hyperparameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sds(sX_1) 1.65 1.16 0.20 4.84 1.24 14 38
sds(sM_1) 0.64 0.67 0.01 2.27 1.49 8 13
sds(t2XM_1) 1.55 1.25 0.12 4.79 1.22 15 71
sds(t2XM_2) 1.66 1.29 0.12 4.77 1.12 23 95
sds(t2XM_3) 32.97 9.09 18.49 57.70 1.31 11 19
Multilevel Hyperparameters:
~id (Number of levels: 100)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.17 0.06 0.06 0.28 1.35 9 12
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.02 0.46 -0.67 1.17 1.53 7 12
sX_1 -1244.86 1903.69 -5453.95 1415.59 2.08 5 12
sM_1 -405.74 1583.39 -3656.27 1897.55 2.06 5 15
t2XM_1 -57.93 201.10 -478.76 228.67 2.07 5 15
t2XM_2 147.79 221.74 -165.05 630.96 2.09 5 12
t2XM_3 6.13 8.55 -6.93 24.11 2.05 5 20
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 1.02 0.02 0.98 1.06 1.16 24 115
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).Code
print(effectsize::standardize_parameters(
bayesianModerationModel,
method = "posthoc"),
digits = 2)# Standardization method: posthoc
Component | Parameter | Std. Median | 95% CI
-------------------------------------------------------------
conditional | (Intercept) | 0.00 | [ -0.00, 0.00]
conditional | sX_1 | -341.64 | [-3282.84, 852.07]
conditional | sM_1 | -67.42 | [-2200.78, 1142.17]
conditional | t2XM_1 | -10.18 | [ -288.18, 137.64]
conditional | t2XM_2 | 38.62 | [ -99.34, 379.78]
conditional | t2XM_3 | 3.22 | [ -4.17, 14.51]
sigma | sigma | 0.62 | [ 0.59, 0.64]Code
performance::model_performance(bayesianModerationModel)Code
performance::r2(bayesianModerationModel)# Bayesian R2 with Compatibility Interval
Conditional R2: 0.627 (95% CI [0.601, 0.650])
Marginal R2: 0.616 (95% CI [0.593, 0.640])6 Setting Priors
6.1 Parameters for Which to Set Priors
Code
get_prior(
Reaction ~ 1 + Days + (1 + Days | Subject),
data = sleepstudy)6.2 Define Priors
6.3 Fit the Model
Fit the model with these priors, and sample from these priors:
Code
fit_sleep4 <- brm(
Reaction ~ 1 + Days + (1 + Days | Subject),
data = sleepstudy,
prior = bprior,
sample_prior = TRUE,
seed = 52242
)7 Evaluate a Hypothesis
Code
# Evid.Ratio is the ratio of P(Days > 7) / P(Days <= 7)
(hyp1 <- hypothesis(fit_sleep4, "Days < 7"))Hypothesis Tests for class b:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (Days)-(7) < 0 2.84 1.59 0.18 5.4 0.04 0.04
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.Code
plot(hyp1)
Code
# Evid.Ratio is the Bayes Factor of the posterior
# vs the prior that Days = 10 is TRUE (Savage-Dickey Ratio)
(hyp2 <- hypothesis(fit_sleep4, "Days = 10"))Hypothesis Tests for class b:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
1 (Days)-(10) = 0 -0.16 1.59 -3.32 2.75 5.17 0.84
Star
1
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.Code
plot(hyp2)
8 Parallel Processing
8.1 Between-Chain Parallelization
Code
fit_sleep3 <- brm(
Reaction ~ 1 + Days + (1 + Days | Subject),
data = sleepstudy,
cores = 4, # between-chain parallelization
seed = 52242
)8.2 Within-Chain Parallelization
https://paul-buerkner.github.io/brms/articles/brms_threading.html (archived at https://perma.cc/NCG3-KV4G)
9 Missing Data Handling
9.1 Multiply Imputed Datasets from mice
Code
?brm_multipleCode
imp <- mice::mice(
mydata,
m = 5,
print = FALSE)Code
fit_imp <- brm_multiple(
bayesianMediationSyntax,
data = imp,
chains = 2)Code
summary(fit_imp) Family: MV(gaussian, gaussian)
Links: mu = identity
mu = identity
Formula: M ~ X + (1 | i | id)
Y ~ X + M + (1 | i | id)
Data: imp (Number of observations: 1000)
Draws: 10 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 10000
Multilevel Hyperparameters:
~id (Number of levels: 100)
Estimate Est.Error l-95% CI u-95% CI
sd(M_Intercept) 0.07 0.05 0.00 0.17
sd(Y_Intercept) 0.11 0.06 0.01 0.23
cor(M_Intercept,Y_Intercept) -0.01 0.56 -0.95 0.94
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI
M_Intercept 0.01 0.03 -0.06 0.07
Y_Intercept 0.06 0.03 -0.01 0.12
M_X 0.51 0.03 0.45 0.58
Y_X 0.03 0.04 -0.04 0.10
Y_M 0.69 0.03 0.63 0.76
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI
sigma_M 1.01 0.02 0.96 1.05
sigma_Y 1.01 0.02 0.96 1.05
Draws were sampled using sampling(NUTS). Overall Rhat and ESS estimates
are not informative for brm_multiple models and are hence not displayed.
Please see ?brm_multiple for how to assess convergence of such models.Code
#print(effectsize::standardize_parameters( # errors out perhaps due to joint model
# fit_imp,
# method = "posthoc"),
# digits = 2)
#performance::model_performance(fit_imp) # errors out perhaps due to joint model
performance::r2(fit_imp) M.R2_Bayes M.R2_Bayes_marginal Y.R2_Bayes Y.R2_Bayes_marginal
0.2051318 0.2002290 0.3925896 0.3846888
attr(,"class")
[1] "r2_bayes_mv"
attr(,"SE")
M.R2_Bayes M.R2_Bayes_marginal Y.R2_Bayes Y.R2_Bayes_marginal
0.02046519 0.01979063 0.01996989 0.01910813
attr(,"CI")
M.R2_Bayes.CI M.R2_Bayes.CI_low
0.9500000 0.1653215
M.R2_Bayes.CI_high M.R2_Bayes_marginal.CI
0.2458624 0.9500000
M.R2_Bayes_marginal.CI_low M.R2_Bayes_marginal.CI_high
0.1621308 0.2400681
Y.R2_Bayes.CI Y.R2_Bayes.CI_low
0.9500000 0.3508358
Y.R2_Bayes.CI_high Y.R2_Bayes_marginal.CI
0.4290966 0.9500000
Y.R2_Bayes_marginal.CI_low Y.R2_Bayes_marginal.CI_high
0.3477870 0.4224221
attr(,"ci_method")
[1] "HDI"
attr(,"robust")
[1] TRUECode
hypothesis(
fit_imp,
"b_M_X * b_Y_M = 0", # indirect effect = a path * b path
class = NULL,
seed = 52242
)Hypothesis Tests for class :
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
1 (b_M_X*b_Y_M) = 0 0.36 0.03 0.3 0.41 NA NA
Star
1 *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.9.2 Imputation on the Fly During Model Fitting
https://paul-buerkner.github.io/brms/articles/brms_missings.html (archived at https://perma.cc/4Y9L-USQR)
Code
?miCode
bayesianRegressionModel <- brm(
bayesianRegressionImputationSyntax,
data = mydata,
seed = 52242
)Code
summary(bayesianRegressionModel) Family: MV(gaussian, gaussian, gaussian)
Links: mu = identity
mu = identity
mu = identity
Formula: X | mi() ~ (1 | i | id)
M | mi() ~ mi(X) + (1 | i | id)
Y | mi() ~ mi(X) + mi(M) + (1 | i | id)
Data: mydata (Number of observations: 1000)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Multilevel Hyperparameters:
~id (Number of levels: 100)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(X_Intercept) 0.10 0.06 0.00 0.22 1.00 1101
sd(M_Intercept) 0.07 0.05 0.00 0.17 1.00 1452
sd(Y_Intercept) 0.11 0.06 0.01 0.22 1.00 1015
cor(X_Intercept,M_Intercept) 0.05 0.49 -0.85 0.87 1.00 3110
cor(X_Intercept,Y_Intercept) 0.07 0.47 -0.81 0.89 1.00 1975
cor(M_Intercept,Y_Intercept) -0.01 0.49 -0.87 0.87 1.00 1829
Tail_ESS
sd(X_Intercept) 1630
sd(M_Intercept) 2033
sd(Y_Intercept) 1768
cor(X_Intercept,M_Intercept) 2910
cor(X_Intercept,Y_Intercept) 2714
cor(M_Intercept,Y_Intercept) 2519
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
X_Intercept -0.02 0.03 -0.09 0.04 1.00 6569 2628
M_Intercept 0.01 0.03 -0.05 0.07 1.00 8761 2815
Y_Intercept 0.06 0.04 -0.01 0.13 1.00 6306 3129
M_miX 0.51 0.03 0.45 0.58 1.00 8105 2943
Y_miX 0.03 0.04 -0.04 0.10 1.00 5638 3094
Y_miM 0.69 0.03 0.63 0.75 1.00 6367 3592
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_X 0.98 0.02 0.94 1.03 1.00 5886 2722
sigma_M 1.01 0.02 0.96 1.06 1.00 7455 2771
sigma_Y 1.01 0.02 0.96 1.05 1.00 5585 2795
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).Code
#print(effectsize::standardize_parameters( # errors out perhaps due to joint model
# bayesianRegressionModel,
# method = "posthoc"),
# digits = 2)
#performance::model_performance(bayesianRegressionModel) # errors out perhaps due to joint model
performance::r2(bayesianRegressionModel) X.R2_Bayes X.R2_Bayes_marginal M.R2_Bayes M.R2_Bayes_marginal
NA NA NA NA
Y.R2_Bayes Y.R2_Bayes_marginal
NA NA
attr(,"class")
[1] "r2_bayes_mv"
attr(,"SE")
X.R2_Bayes X.R2_Bayes_marginal M.R2_Bayes M.R2_Bayes_marginal
NA NA NA NA
Y.R2_Bayes Y.R2_Bayes_marginal
NA NA
attr(,"CI")
X.R2_Bayes.CI X.R2_Bayes.CI_low
0.95 NA
X.R2_Bayes.CI_high X.R2_Bayes_marginal.CI
NA 0.95
X.R2_Bayes_marginal.CI_low X.R2_Bayes_marginal.CI_high
NA NA
M.R2_Bayes.CI M.R2_Bayes.CI_low
0.95 NA
M.R2_Bayes.CI_high M.R2_Bayes_marginal.CI
NA 0.95
M.R2_Bayes_marginal.CI_low M.R2_Bayes_marginal.CI_high
NA NA
Y.R2_Bayes.CI Y.R2_Bayes.CI_low
0.95 NA
Y.R2_Bayes.CI_high Y.R2_Bayes_marginal.CI
NA 0.95
Y.R2_Bayes_marginal.CI_low Y.R2_Bayes_marginal.CI_high
NA NA
attr(,"ci_method")
[1] "HDI"
attr(,"robust")
[1] TRUECode
hypothesis(
bayesianRegressionModel,
"bsp_M_miX * bsp_Y_miM = 0", # indirect effect = a path * b path
class = NULL,
seed = 52242
)Hypothesis Tests for class :
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (bsp_M_miX*bsp_Y_... = 0 0.35 0.03 0.3 0.41 NA
Post.Prob Star
1 NA *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.10 Session Info
Code
R version 4.6.1 (2026-06-24)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.4 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
locale:
[1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
[4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
[7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
time zone: UTC
tzcode source: system (glibc)
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] future_1.75.0 see_0.14.1 report_0.6.4
[4] parameters_0.29.2 performance_0.17.1 modelbased_0.16.0
[7] insight_1.5.2 effectsize_1.0.3 datawizard_1.3.1
[10] correlation_0.8.8 easystats_0.7.6 mice_3.19.0
[13] bayestestR_0.18.1 brms_2.23.0 Rcpp_1.1.2
[16] rstan_2.32.7 StanHeaders_2.32.10 lme4_2.0-6
[19] Matrix_1.7-5
loaded via a namespace (and not attached):
[1] Rdpack_2.6.6 gridExtra_2.3.1 inline_0.3.21
[4] rlang_1.3.0 magrittr_2.0.5 otel_0.2.0
[7] matrixStats_1.5.0 compiler_4.6.1 mgcv_1.9-4
[10] loo_2.10.1 callr_3.8.0 vctrs_0.7.3
[13] reshape2_1.4.5 stringr_1.6.0 pkgconfig_2.0.3
[16] shape_1.4.6.1 fastmap_1.2.0 backports_1.5.1
[19] labeling_0.4.3 rmarkdown_2.31 ps_1.9.3
[22] nloptr_2.2.1 purrr_1.2.2 xfun_0.60
[25] glmnet_5.0 jomo_2.7-6 jsonlite_2.0.0
[28] pan_2.0 broom_1.0.13 parallel_4.6.1
[31] R6_2.6.1 stringi_1.8.9 RColorBrewer_1.1-3
[34] parallelly_1.48.0 boot_1.3-32 rpart_4.1.27
[37] iterators_1.0.14 knitr_1.51 future.apply_1.20.2
[40] bayesplot_1.15.0 splines_4.6.1 nnet_7.3-20
[43] tidyselect_1.2.1 rstudioapi_0.19.0 abind_1.4-8
[46] yaml_2.3.12 codetools_0.2-20 processx_3.9.0
[49] listenv_1.0.0 pkgbuild_1.4.8 lattice_0.22-9
[52] tibble_3.3.1 plyr_1.8.9 withr_3.0.3
[55] bridgesampling_1.2-1 S7_0.2.2 posterior_1.7.0
[58] coda_0.19-4.1 evaluate_1.0.5 survival_3.8-6
[61] RcppParallel_6.2.0 pillar_1.11.1 tensorA_0.36.2.1
[64] checkmate_2.3.4 foreach_1.5.2 stats4_4.6.1
[67] reformulas_0.4.4 distributional_0.8.1 generics_0.1.4
[70] ggplot2_4.0.3 rstantools_2.7.0 scales_1.4.0
[73] minqa_1.2.8 globals_0.19.1 glue_1.8.1
[76] tools_4.6.1 mvtnorm_1.4-2 grid_4.6.1
[79] tidyr_1.3.2 rbibutils_2.4.1 QuickJSR_1.10.0
[82] nlme_3.1-169 cli_3.6.6 Brobdingnag_1.2-9
[85] dplyr_1.2.1 gtable_0.3.6 digest_0.6.39
[88] htmlwidgets_1.6.4 farver_2.1.2 htmltools_0.5.9
[91] lifecycle_1.0.5 mitml_0.4-5 MASS_7.3-65 