Random Effects with Formula Syntax

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

František Bartoš

29th of October 2026

BayesTools can generate JAGS (Plummer, 2003) mixed models directly from formulas with random-effect terms. Its random-coefficient syntax follows the conventions used by lme4 (Bates et al., 2015), and prior_random() allows users to specify the corresponding prior distributions, e.g., random-effect standard deviation priors, correlation priors, variance-allocation priors.

This vignette illustrates how to use the functionality on example datasets distributed with lme4. The first example uses the sleepstudy data and shows random intercepts, correlated and independent random slopes, and Bayes factor tests for the random slope and the correlation, comparing lme4, brms (Bürkner, 2017), and BayesTools. The second example continues with the sleepstudy data and shows a heterogeneous compound-symmetry covariance over the study days and a test of its correlation against a specific value. The third example uses the Pastes data and shows nested random intercepts. The fourth example uses the Penicillin data and shows crossed random effects with independent SD priors and a variance allocation prior.

Formula Semantics

One major difference from standard R formulas is that BayesTools formulas keeps intercept control and factor contrasts separate. In a fixed effects formula, BayesTools’ 0 + factor or factor - 1 fixes the intercept at zero but preserves the factor contrast specified by prior_factor(). For example, an independent factor prior creates one coefficient per factor level, while treatment and mean-difference priors retain their own contrast despite the same no-intercept formula. The reason behind this choice is that prior specification, i.e., factor contrast, should not change when other coefficients are modified.

Random effects specification differs between two types of random structures:

Structure syntax	Left-side role	Marginal SDs	Correlation
Plain bar, `us()` / `un()`	coefficient formula	one per generated coefficient (one per mean-difference/orthonormal factor term)	unrestricted LKJ correlation matrix
Double bar, `diag()`	coefficient formula	one per generated coefficient (one per mean-difference/orthonormal factor term)	independent
`id()`	coefficient formula	one shared SD	independent
`cs()`	one or more discrete indices	one shared SD	one common pairwise correlation
`hcs()`	one or more discrete indices	one per index level	one common pairwise correlation
`ar1()` / `ar()`	one discrete ordered index	one shared SD	AR(1), `cor^abs(i - j)`
`har()`	one discrete ordered index	one per index level	AR(1), `cor^abs(i - j)`
`car()`	one numeric-distance coordinate	one shared SD	continuous AR(1), `cor^abs(t_i - t_j)`

The first three structures operate on coefficient columns. Their 1, 0, and -1 terms include or remove the random intercept, and their left side can contain continuous slopes, factor slopes, and interactions. A random factor reuses its already specified fixed effect contrast by default. Set random_block(contrasts = c(factor_name = "contrast")) when the random contrast should be explicit or differ from the fixed contrast. Therefore us(0 + group | study) with an independent block contrast means one correlated random coefficient per group level and no random intercept; 0 + group by itself does not force by-group indicators (as standard R formula does). A mean-difference or orthonormal factor term is one exchangeable unit: its contrast columns are an arbitrary basis of the level deviations, so they share one SD, sd(group), and are uncorrelated with each other. Under us() the unit keeps its correlations with the other coefficients, which are reported for each level, cor(intercept, group[a]), together with the multiple correlation cor_total(intercept, group); a us() block may contain at most one such term with three or more levels.

The structured CS/HCS and AR1/HAR families operate on discrete index levels instead. They reject 1, 0, -1, and block contrast overrides. They accept factor, character, numeric/integer, or logical columns. Existing factor level order is retained; otherwise sorted unique values define the levels and AR(1) order. CS/HCS can combine several indices with +; AR1/HAR require exactly one. CAR requires one finite numeric/integer coordinate, or an ordered factor with numeric labels, and uses actual distances.

Correlation prior distributions

Correlation priors are completed only after the formula determines the structure and its dimension. An omitted US/UN correlation prior becomes LKJ(1). An omitted scalar correlation prior is uniform over the complete admissible raw-correlation interval: (-1 / (K - 1), 1) for CS/HCS with K levels, (-1, 1) for AR1/HAR, and (0, 1) for CAR. The correlation priors are specified via random_covariance(cor = ..., cor_scale = ...). The default cor_scale = "fisher_z" places the prior on Fisher’s z; "cor" places it on the raw correlation and "logit" first rescales the correlation from its admissible interval to (0, 1) and then applies the logit transformation.

Standard deviation prior distributions

Random-effect SD priors have no default and must always be supplied, either once for all blocks in prior_random(sd = ...) or per block in random_block(sd = ...).

suppressPackageStartupMessages({
  library(BayesTools)
  library(lme4)
})

Random Intercepts and Slopes

The sleepstudy data contain reaction times of 18 subjects measured over ten days of sleep deprivation. We illustrate BayesTools features in a few simple steps: 1) we specify a random intercept model, 2) extend the model with random slope, 3) compare the models using a Bayes factor, 4) specify a model without the random intercept-slope correlation, and 5) test for presence vs. absence of the correlation via a Bayes factor. For comparison we fit corresponding models with lme4 (Bates et al., 2015) and brms (Bürkner, 2017).

We standardize the response and the predictor to simplify prior specification for the example. This allows us to use standard normal priors for the centered intercept and the standardized effect of Days, half-normal prior for the residual SD using the scaled outcome, a half-normal prior for random-effect SD, and an LKJ(1) prior for the correlation between the random intercept and slope. For brms and lme4 we use the explicitly standardized variables defined below, BayesTools allows us to standardize the predictor directly in the model fitting interface.

data("sleepstudy", package = "lme4")
sleepstudy$Subject    <- factor(sleepstudy$Subject)
sleepstudy$Reaction_z <- as.numeric(scale(sleepstudy$Reaction))
sleepstudy$Days_s     <- as.numeric(scale(sleepstudy$Days))
head(sleepstudy)
#>   Reaction Days Subject Reaction_z     Days_s
#> 1 249.5600    0     308 -0.8689681 -1.5623409
#> 2 258.7047    1     308 -0.7066229 -1.2151540
#> 3 250.8006    2     308 -0.8469438 -0.8679672
#> 4 321.4398    3     308  0.4071084 -0.5207803
#> 5 356.8519    4     308  1.0357766 -0.1735934
#> 6 414.6901    5     308  2.0625736  0.1735934

Random Intercepts

In lme4, a random intercept model is specified using the (1 | Subject) syntax in the lme4::lmer function.

lme4_intercept <- lme4::lmer(
  Reaction_z ~ Days_s + (1 | Subject),
  data = sleepstudy, REML = FALSE
)

lme4_intercept
#> Linear mixed model fit by maximum likelihood  ['lmerMod']
#> Formula: Reaction_z ~ Days_s + (1 | Subject)
#>    Data: sleepstudy
#>       AIC       BIC    logLik -2*log(L)  df.resid 
#>  350.8448  363.6166 -171.4224  342.8448       176 
#> Random effects:
#>  Groups   Name        Std.Dev.
#>  Subject  (Intercept) 0.6393  
#>  Residual             0.5485  
#> Number of obs: 180, groups:  Subject, 18
#> Fixed Effects:
#> (Intercept)       Days_s  
#>   1.937e-16    5.352e-01

brms uses the same syntax. Here we explicitly specify prior distributions for the individual parameters. We call brms::prior() with its namespace because BayesTools also exports prior() function.

brms_priors <- c(
  brms::prior(normal(0, 1), class = Intercept), # intercept
  brms::prior(normal(0, 1), class = b),         # Days_s coefficient
  brms::prior(normal(0, 1), class = sigma),     # residual sd
  brms::prior(normal(0, 1), class = sd)         # random intercept sd
)

brms bounds SDs at zero, so the normal priors on sigma and sd become half-normal priors. The model is estimated via brms::brm function and save_pars allows us to retain all parameters for bridge sampling used later.

brms_intercept <- brms::brm(
  Reaction_z ~ Days_s + (1 | Subject),
  data = sleepstudy, prior = brms_priors, save_pars = brms::save_pars(all = TRUE),
  chains = 2, iter = 3000, warmup = 1000, seed = 1, refresh = 0, silent = 2
)

brms_intercept
#>  Family: gaussian 
#>   Links: mu = identity 
#> Formula: Reaction_z ~ Days_s + (1 | Subject) 
#>    Data: sleepstudy (Number of observations: 180) 
#>   Draws: 2 chains, each with iter = 3000; 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)     0.70      0.13     0.49     1.01 1.00      857     1443
#> 
#> Regression Coefficients:
#>           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept     0.00      0.17    -0.32     0.33 1.00      608     1091
#> Days_s        0.54      0.04     0.46     0.62 1.00     3881     2787
#> 
#> Further Distributional Parameters:
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma     0.55      0.03     0.50     0.62 1.00     2656     2799
#> 
#> 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).

BayesTools first requires separate specification of the JAGS likelihood and the directly supplied variables (Reaction_z and N).

sleep_likelihood <-
"model{
  for(i in 1:N){
    Reaction_z[i] ~ dnorm(mu[i], pow(sigma, -2))
  }
}
"
sleep_data <- list(Reaction_z = sleepstudy$Reaction_z, N = nrow(sleepstudy))

The mu parameter is specified via a formula that becomes assigned to fitting function’s formula_list argument later (sigma parameter is specified directly in the prior_list of the same function).

sleep_formula_list <- list(mu = ~ Days + (1 | Subject))

The BayesTools formula interface requires specifying prior distributions for the fixed and random effects parameters separately. Prior distributions for the fixed effect intercept and predictor enter the fitting function via formula_prior_list argument, here specified in the named sleep_fixed_priors list. Prior distributions for the random effects intercept enter the fitting function via formula_random_prior_list argument, here specified in the named sleep_sd_priors list using the prior_random function. prior_random() assigns priors per random-effect block which takes its name from the grouping variable Subject.

unit_prior    <- prior(distribution = "normal", parameters = list(mean = 0, sd = 1))
unit_sd_prior <- prior(distribution = "normal", parameters = list(mean = 0, sd = 1), truncation = list(0, Inf))

sleep_prior_list   <- list(sigma = unit_sd_prior)
sleep_fixed_priors <- list(mu = list(intercept = unit_prior, Days = unit_prior))
sleep_sd_priors    <- list(mu = prior_random(Subject = random_block(sd = unit_sd_prior)))

The model is estimated via JAGS_fit function. formula_scale_list standardizes Days predictor before fitting, so the Days prior refers to the same standardized effect as Days_s in lme4 and brms (but allows us to obtain non-standardized effect later).

fit_sleep_intercept <- JAGS_fit(
  model_syntax = sleep_likelihood, data = sleep_data, prior_list = sleep_prior_list,
  formula_list = sleep_formula_list, formula_data_list = list(mu = sleepstudy),
  formula_prior_list = sleep_fixed_priors, formula_random_prior_list = sleep_sd_priors,
  formula_scale_list = list(mu = list(Days = TRUE)),
  chains = 2, adapt = 1000, burnin = 1000, sample = 5000, seed = 1
)

JAGS_estimates_table(fit_sleep_intercept)
#>                     Mean    SD  0.025   0.5 0.975 error(MCMC) error(MCMC)/SD  ESS R-hat
#> (mu) intercept     0.004 0.174 -0.332 0.002 0.361     0.01029          0.059  292 1.007
#> (mu) Days          0.534 0.042  0.452 0.534 0.614     0.00042          0.010 9793 1.000
#> sigma              0.554 0.031  0.496 0.553 0.619     0.00043          0.014 5311 1.000
#> (mu) sd(intercept) 0.702 0.131  0.494 0.687 1.006     0.00538          0.041  597 1.000

Table below shows close agreement in the parameter estimates from all three packages. The lme4 maximum likelihood estimates of the random-effect SD slightly differs from the Bayesian approaches because their posterior distributions are wide and skewed with only 18 subjects.

compare_sleep_estimates(bayes_model = fit_sleep_intercept, brms_model = brms_intercept, lme4_model = lme4_intercept)
#>       parameter  lme4  brms BayesTools
#> 1     intercept 0.000 0.002      0.004
#> 2          Days 0.535 0.535      0.534
#> 3         sigma 0.548 0.553      0.554
#> 4 sd(intercept) 0.639 0.697      0.702

Random Slopes

In lme4, adding Days_s to the left side of the random-effect term, (1 + Days_s | Subject), gives each subject its own Days_s effect. The single bar | also estimates the correlation between the random intercepts and slopes.

lme4_slope <- lme4::lmer(
  Reaction_z ~ Days_s + (1 + Days_s | Subject),
  data = sleepstudy, REML = FALSE
)

lme4_slope
#> Linear mixed model fit by maximum likelihood  ['lmerMod']
#> Formula: Reaction_z ~ Days_s + (1 + Days_s | Subject)
#>    Data: sleepstudy
#>       AIC       BIC    logLik -2*log(L)  df.resid 
#>  312.7055  331.8632 -150.3527  300.7055       174 
#> Random effects:
#>  Groups   Name        Std.Dev. Corr 
#>  Subject  (Intercept) 0.6467        
#>           Days_s      0.2923   0.76 
#>  Residual             0.4543        
#> Number of obs: 180, groups:  Subject, 18
#> Fixed Effects:
#> (Intercept)       Days_s  
#>      0.0000       0.5352

brms uses the same syntax. The sd prior now applies to both random-effect SDs (brms::prior(..., class = sd, coef = Days_s) would set a separate prior for the Days_s SD), and we add a prior distribution for the correlation.

brms_slope_priors <- c(
  brms_priors,
  brms::prior(lkj(1), class = cor)              # random intercept-slope correlation
)

brms_slope <- brms::brm(
  Reaction_z ~ Days_s + (1 + Days_s | Subject),
  data = sleepstudy, prior = brms_slope_priors, save_pars = brms::save_pars(all = TRUE),
  chains = 2, iter = 3000, warmup = 1000, seed = 1, refresh = 0, silent = 2
)

brms_slope
#>  Family: gaussian 
#>   Links: mu = identity 
#> Formula: Reaction_z ~ Days_s + (1 + Days_s | Subject) 
#>    Data: sleepstudy (Number of observations: 180) 
#>   Draws: 2 chains, each with iter = 3000; 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)             0.70      0.13     0.50     1.01 1.00     1056     1996
#> sd(Days_s)                0.32      0.07     0.21     0.49 1.00     1440     1992
#> cor(Intercept,Days_s)     0.67      0.17     0.28     0.91 1.00     1450     2048
#> 
#> Regression Coefficients:
#>           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept    -0.01      0.16    -0.33     0.31 1.00      626     1106
#> Days_s        0.53      0.09     0.35     0.70 1.00     1046     1349
#> 
#> Further Distributional Parameters:
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma     0.46      0.03     0.41     0.51 1.00     3459     2757
#> 
#> 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).

BayesTools uses the same random-effect term in the mu formula. The likelihood, the data, and the prior distributions for sigma and the fixed effects remain unchanged.

sleep_slope_formula_list <- list(mu = ~ Days + (1 + Days | Subject))

The random_block() function adds the correlation prior next to the SD prior. The SD prior applies to both random-effect SDs unless the terms argument overrides it for a specific term; for example, terms = list(Days = ...) sets a separate prior for the Days SD, which we use for a hypothesis test via a spike-and-slab prior later.

sleep_slope_priors <- list(mu = prior_random(Subject = random_block(
  sd  = unit_sd_prior, 
  cor = prior_lkj(eta = 1)
)))

fit_sleep_slope <- JAGS_fit(
  model_syntax = sleep_likelihood, data = sleep_data, prior_list = sleep_prior_list,
  formula_list = sleep_slope_formula_list, formula_data_list = list(mu = sleepstudy),
  formula_prior_list = sleep_fixed_priors, formula_random_prior_list = sleep_slope_priors,
  formula_scale_list = list(mu = list(Days = TRUE)),
  chains = 2, adapt = 1000, burnin = 1000, sample = 5000, seed = 1
)

JAGS_estimates_table(fit_sleep_slope)
#>                           Mean    SD  0.025   0.5 0.975 error(MCMC) error(MCMC)/SD  ESS R-hat
#> (mu) intercept           0.026 0.166 -0.298 0.024 0.351     0.01219          0.074  184 1.014
#> (mu) Days                0.541 0.087  0.364 0.540 0.710     0.00367          0.042  579 1.004
#> sigma                    0.459 0.028  0.409 0.458 0.518     0.00046          0.016 3756 1.000
#> (mu) sd(intercept)       0.704 0.132  0.496 0.688 1.003     0.00684          0.052  373 1.002
#> (mu) sd(Days)            0.326 0.073  0.209 0.318 0.493     0.00202          0.028 1302 1.001
#> (mu) cor(intercept,Days) 0.669 0.173  0.247 0.700 0.918     0.00574          0.033  909 1.002

Table below again shows close agreement between the brms and BayesTools posterior means.

compare_sleep_estimates(bayes_model = fit_sleep_slope, brms_model = brms_slope, lme4_model = lme4_slope)
#>             parameter  lme4   brms BayesTools
#> 1           intercept 0.000 -0.007      0.026
#> 2                Days 0.535  0.530      0.541
#> 3               sigma 0.454  0.459      0.459
#> 4       sd(intercept) 0.647  0.701      0.704
#> 5            sd(Days) 0.292  0.322      0.326
#> 6 cor(intercept,Days) 0.759  0.670      0.669

Because formula_scale_list standardized Days inside BayesTools, transform_scaled = TRUE reports the estimates on the original Days scale. We compare them with lme4 fitted to the original Days. The intercept, its SD, and the correlation then refer to the first day rather than the mean day, while the prior distributions still act on the standardized scale.

lme4_slope_days <- lme4::lmer(
  Reaction_z ~ Days + (1 + Days | Subject),
  data = sleepstudy, REML = FALSE
)
compare_sleep_estimates(bayes_model = fit_sleep_slope, lme4_model = lme4_slope_days, transform_scaled = TRUE)
#>             parameter   lme4 BayesTools
#> 1           intercept -0.836     -0.819
#> 2                Days  0.186      0.188
#> 3               sigma  0.454      0.459
#> 4       sd(intercept)  0.422      0.521
#> 5            sd(Days)  0.101      0.113
#> 6 cor(intercept,Days)  0.081     -0.037

Testing the Random Slope

We can also test whether subjects differ in how their reaction times change over the days. This can be done by using a Bayes factor to compare the random-slope model with the random-intercept model. In brms we can obtain Bayes factors via bridge sampling (Gronau et al., 2020) to compare the two model fits. In BayesTools we can also use product-space approach (Lodewyckx et al., 2011) and obtain the Bayes factor from posterior samples by using a spike-and-slab prior distribution. While the product-space approach is less precise, it allows us to scale up the inference to examine multiple predictors simultaneously and automatically returns Bayesian model-averaged estimates.

First, we obtain the Bayes factor via bridge sampling for both packages. brms models are passed directly to the brms::bridge_sampler function.

marglik_brms_intercept <- brms::bridge_sampler(brms_intercept, silent = TRUE)
marglik_brms_slope     <- brms::bridge_sampler(brms_slope, silent = TRUE)

BayesTools models are passed to the JAGS_bridgesampling function, which additionally requires the log-likelihood of the data as a function of the parameters. BayesTools automatically prepares the parameters (mu and sigma) and evaluates the prior densities.

log_posterior <- function(parameters, data){
  sum(dnorm(
    x    = data[["Reaction_z"]],
    mean = parameters[["mu"]],
    sd   = parameters[["sigma"]],
    log  = TRUE
  ))
}

marglik_sleep_intercept <- JAGS_bridgesampling(
  fit_sleep_intercept, log_posterior = log_posterior, data = sleep_data,
  prior_list = sleep_prior_list, seed = 1)
marglik_sleep_slope <- JAGS_bridgesampling(
  fit_sleep_slope, log_posterior = log_posterior, data = sleep_data,
  prior_list = sleep_prior_list, seed = 1)

Table below shows the log Bayes factors from both packages, which agree up to the bridge sampling error and show extreme evidence for the presence of the random slope.

log_BF <- function(marglik_1, marglik_0) bridgesampling::bf(marglik_1, marglik_0, log = TRUE)$bf

data.frame(
  method = c("brms, bridge sampling", "BayesTools, bridge sampling"),
  log_BF = c(log_BF(marglik_brms_slope,  marglik_brms_intercept),
             log_BF(marglik_sleep_slope, marglik_sleep_intercept))
)
#>                        method   log_BF
#> 1       brms, bridge sampling 17.34698
#> 2 BayesTools, bridge sampling 17.37809

Second, BayesTools allows us to obtain the Bayes factor using a spike-and-slab prior (Kuo & Mallick, 1998) within a single model fit. An inclusion indicator switches the Days SD between zero (the spike) and the half-normal prior (the slab); with the SD at zero, the model reduces to the random-intercept model. The terms argument assigns a prior to each SD of the block: the half-normal prior to the intercept SD and the spike-and-slab prior to the Days SD. Because terms covers both SDs, the block needs no sd argument. random_monitor(latent = FALSE) reduces the size of the stored fit by omitting the subject-level random effects samples.

sleep_spike_priors <- list(
  mu = prior_random(
    Subject = random_block(
      cor   = prior_lkj(eta = 1),
      terms = list(
        intercept = unit_sd_prior,
        Days      = prior_spike_and_slab(prior_parameter = unit_sd_prior)
      ),
      monitor = random_monitor(latent = FALSE)
    )
  )
)

fit_sleep_slope_spike <- JAGS_fit(
  model_syntax = sleep_likelihood, data = sleep_data, prior_list = sleep_prior_list,
  formula_list = sleep_slope_formula_list, formula_data_list = list(mu = sleepstudy),
  formula_prior_list = sleep_fixed_priors, formula_random_prior_list = sleep_spike_priors,
  formula_scale_list = list(mu = list(Days = TRUE)),
  chains = 2, adapt = 1000, burnin = 1000, sample = 5000, seed = 1
)

The JAGS_inference_table function summarizes the inclusion of the random slope.

JAGS_inference_table(fit_sleep_slope_spike)
#>                        Prior prob. Post. prob. Inclusion BF
#> (mu) Subject: sd(Days)       0.500       1.000    >9999.000

All posterior draws include the random slope, so the inclusion Bayes factor is bounded by the number of draws. The spike-and-slab prior agrees with conclusions from bridge sampling (extreme evidence), but it cannot precisely quantify such large Bayes factors.

Removing the Random Intercept-Slope Correlation

In lme4, the double bar || removes the correlation between the random intercepts and slopes.

lme4_independent <- lme4::lmer(
  Reaction_z ~ Days_s + (1 + Days_s || Subject),
  data = sleepstudy, REML = FALSE
)

lme4_independent
#> Linear mixed model fit by maximum likelihood  ['lmerMod']
#> Formula: Reaction_z ~ Days_s + ((1 | Subject) + (0 + Days_s | Subject))
#>    Data: sleepstudy
#>       AIC       BIC    logLik -2*log(L)  df.resid 
#>  321.2096  337.1744 -155.6048  311.2096       175 
#> Random effects:
#>  Groups    Name        Std.Dev.
#>  Subject   (Intercept) 0.6467  
#>  Subject.1 Days_s      0.2923  
#>  Residual              0.4543  
#> Number of obs: 180, groups:  Subject, 18
#> Fixed Effects:
#> (Intercept)       Days_s  
#>  -3.497e-16    5.352e-01

brms uses the same syntax. Without the correlation, the prior distributions from the random-intercept model suffice; the sd prior applies to both random-effect SDs.

brms_independent <- brms::brm(
  Reaction_z ~ Days_s + (1 + Days_s || Subject),
  data = sleepstudy, prior = brms_priors, save_pars = brms::save_pars(all = TRUE),
  chains = 2, iter = 3000, warmup = 1000, seed = 1, refresh = 0, silent = 2
)

brms_independent
#>  Family: gaussian 
#>   Links: mu = identity 
#> Formula: Reaction_z ~ Days_s + (1 + Days_s || Subject) 
#>    Data: sleepstudy (Number of observations: 180) 
#>   Draws: 2 chains, each with iter = 3000; 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)     0.70      0.13     0.50     0.99 1.00      790     1155
#> sd(Days_s)        0.33      0.08     0.21     0.51 1.00     1018     1733
#> 
#> Regression Coefficients:
#>           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept    -0.00      0.17    -0.34     0.33 1.00      517      779
#> Days_s        0.53      0.08     0.36     0.69 1.00      854     1391
#> 
#> Further Distributional Parameters:
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma     0.46      0.03     0.41     0.52 1.00     2552     3030
#> 
#> 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).

BayesTools uses the same double bar (or the equivalent diag() covariance tag) in the mu formula. The random-effect prior distributions from the random-intercept model, sleep_sd_priors, again suffice.

sleep_independent_formula_list <- list(mu = ~ Days + (1 + Days || Subject))

fit_sleep_independent <- JAGS_fit(
  model_syntax = sleep_likelihood, data = sleep_data, prior_list = sleep_prior_list,
  formula_list = sleep_independent_formula_list, formula_data_list = list(mu = sleepstudy),
  formula_prior_list = sleep_fixed_priors, formula_random_prior_list = sleep_sd_priors,
  formula_scale_list = list(mu = list(Days = TRUE)),
  chains = 2, adapt = 1000, burnin = 1000, sample = 5000, seed = 1
)

JAGS_estimates_table(fit_sleep_independent)
#>                      Mean    SD  0.025    0.5 0.975 error(MCMC) error(MCMC)/SD  ESS R-hat
#> (mu) intercept     -0.007 0.169 -0.351 -0.002 0.303     0.01170          0.069  206 1.012
#> (mu) Days           0.532 0.082  0.372  0.531 0.692     0.00268          0.033  943 1.003
#> sigma               0.458 0.027  0.409  0.456 0.515     0.00042          0.016 4163 1.000
#> (mu) sd(intercept)  0.706 0.140  0.499  0.684 1.037     0.00840          0.060  312 1.009
#> (mu) sd(Days)       0.327 0.073  0.212  0.317 0.494     0.00192          0.026 1459 1.002

Table below shows close agreement in the parameter estimates from all three packages.

compare_sleep_estimates(bayes_model = fit_sleep_independent, brms_model = brms_independent, lme4_model = lme4_independent)
#>       parameter  lme4  brms BayesTools
#> 1     intercept 0.000 0.000     -0.007
#> 2          Days 0.535 0.531      0.532
#> 3         sigma 0.454 0.459      0.458
#> 4 sd(intercept) 0.647 0.700      0.706
#> 5      sd(Days) 0.292 0.329      0.327

Testing the Correlation

The Bayes factor comparing the correlated and the independent random-slope models tests the correlation. Because Days is standardized, the correlation refers to the random intercept at the mean day. BayesTools does not provide a spike-and-slab prior for LKJ correlations, so we use bridge sampling only (later example shows how to test scalar correlations).

marglik_brms_independent  <- brms::bridge_sampler(brms_independent, silent = TRUE)
marglik_sleep_independent <- JAGS_bridgesampling(
  fit_sleep_independent, log_posterior = log_posterior, data = sleep_data,
  prior_list = sleep_prior_list, seed = 1)

Table below shows the Bayes factors from both packages, which agree up to the bridge sampling error and show strong evidence for the correlation.

data.frame(
  method = c("brms", "BayesTools"),
  BF = c(bridgesampling::bf(marglik_brms_slope,  marglik_brms_independent)$bf,
         bridgesampling::bf(marglik_sleep_slope, marglik_sleep_independent)$bf)
)
#>       method       BF
#> 1       brms 19.43482
#> 2 BayesTools 22.08361

Structured Random-Effect Correlations

The random intercept and slope model implies a specific covariance between a subject’s measurements on different days: subjects differ in their level and in a linear trend, so the SD of their deviations changes over the days. Structured covariances model the covariance between days directly, without assuming that each subject follows a straight line. They use the levels of a discrete index instead of coefficient columns: in addition to the fixed Days slope, each subject receives one random effect per day, and the covariance structure determines how these day-specific effects are correlated. A random intercept alone corresponds to compound symmetry, in which all pairs of days share one correlation and all days one SD. Heterogeneous compound symmetry (HCS) keeps the common correlation but estimates a separate SD for each day. We fit the HCS model, compare its day-specific SDs with those implied by the random-slope model, and test whether its correlation equals a specific value.

sleepstudy$Days_f <- factor(sleepstudy$Days)

Heterogeneous Compound Symmetry

In lme4, the cs() term allows us to specify a compound-symmetry covariance. It estimates a separate SD for each column by default when using the 0 + Days_f model-matrix notation

lme4_hcs <- lme4::lmer(
  Reaction_z ~ Days_s + cs(0 + Days_f | Subject),
  data = sleepstudy, REML = FALSE
)
lme4_hcs
#> Linear mixed model fit by maximum likelihood  ['lmerMod']
#> Formula: Reaction_z ~ Days_s + cs(0 + Days_f | Subject)
#>    Data: sleepstudy
#>       AIC       BIC    logLik -2*log(L)  df.resid 
#>  415.5300  600.7215 -149.7650  299.5300       122 
#> Random effects:
#>  Groups   Name    Std.Dev. Corr     
#>  Subject  Days_f0 0.3777   0.82 (cs)
#>           Days_f1 0.4395            
#>           Days_f2 0.3569            
#>           Days_f3 0.5820            
#>           Days_f4 0.6554            
#>           Days_f5 0.8402            
#>           Days_f6 1.0430            
#>           Days_f7 0.7822            
#>           Days_f8 0.9856            
#>           Days_f9 1.1359            
#>  Residual         0.3691            
#> Number of obs: 180, groups:  Subject, 18
#> Fixed Effects:
#> (Intercept)       Days_s  
#>     -0.1068       0.4906

BayesTools uses an explicit hcs() term in the mu formula. The left side automatically expands the factor into one column per level (the 0 + notation does not apply here). The likelihood, the data, and the prior distributions for sigma and the fixed effects remain unchanged from the previous example.

sleep_hcs_formula_list <- list(mu = ~ Days + hcs(Days_f | Subject))

The random_block() function sets the SD prior which applies to every day-specific SD independently. Without a cor argument, the correlation receives the default uniform prior over its admissible interval, (-1/9, 1) for an HCS correlation with 10 days. random_monitor() allows us to stores only the correlation posterior samples.

day_random_priors <- list(mu = prior_random(Subject = random_block(
  sd      = unit_sd_prior,
  monitor = random_monitor(latent = FALSE, coefficients = FALSE, correlation = TRUE)
)))

fit_sleep_hcs <- JAGS_fit(
  model_syntax = sleep_likelihood, data = sleep_data, prior_list = sleep_prior_list,
  formula_list = sleep_hcs_formula_list, formula_data_list = list(mu = sleepstudy),
  formula_prior_list = sleep_fixed_priors, formula_random_prior_list = day_random_priors,
  formula_scale_list = list(mu = list(Days = TRUE)),
  chains = 2, adapt = 1000, burnin = 2000, sample = 5000, seed = 31
)

JAGS_estimates_table(fit_sleep_hcs)
#>                      Mean    SD  0.025    0.5 0.975 error(MCMC) error(MCMC)/SD  ESS R-hat
#> (mu) intercept     -0.074 0.139 -0.362 -0.070 0.190     0.00966          0.069  208 1.000
#> (mu) Days           0.502 0.060  0.377  0.504 0.615     0.00302          0.050  397 1.002
#> sigma               0.396 0.056  0.282  0.400 0.500     0.00387          0.069  210 1.007
#> (mu) sd(Days_f[0])  0.393 0.145  0.126  0.390 0.699     0.00555          0.038  682 1.003
#> (mu) sd(Days_f[1])  0.466 0.146  0.202  0.459 0.777     0.00456          0.031 1020 1.006
#> (mu) sd(Days_f[2])  0.391 0.131  0.149  0.385 0.668     0.00324          0.025 1645 1.001
#> (mu) sd(Days_f[3])  0.621 0.146  0.365  0.609 0.940     0.00434          0.030 1125 1.000
#> (mu) sd(Days_f[4])  0.709 0.151  0.450  0.700 1.045     0.00486          0.032 1075 1.001
#> (mu) sd(Days_f[5])  0.897 0.178  0.596  0.881 1.288     0.00660          0.037  734 1.003
#> (mu) sd(Days_f[6])  1.097 0.214  0.736  1.079 1.570     0.00805          0.038  708 1.003
#> (mu) sd(Days_f[7])  0.838 0.172  0.537  0.822 1.217     0.00576          0.033  911 1.002
#> (mu) sd(Days_f[8])  1.055 0.194  0.724  1.038 1.483     0.00739          0.038  692 1.000
#> (mu) sd(Days_f[9])  1.194 0.221  0.823  1.171 1.673     0.00944          0.043  552 1.004
#> (mu) cor            0.842 0.077  0.670  0.850 0.970     0.00561          0.073  192 1.007

Table below compares the estimates with the lme4 cs() fit.

compare_structured_estimates(lme4_model = lme4_hcs, bayes_model = fit_sleep_hcs)
#>        parameter   lme4 BayesTools
#> 1      intercept -0.107     -0.074
#> 2           Days  0.491      0.502
#> 3          sigma  0.369      0.396
#> 4  sd(Days_f[0])  0.378      0.393
#> 5  sd(Days_f[1])  0.440      0.466
#> 6  sd(Days_f[2])  0.357      0.391
#> 7  sd(Days_f[3])  0.582      0.621
#> 8  sd(Days_f[4])  0.655      0.709
#> 9  sd(Days_f[5])  0.840      0.897
#> 10 sd(Days_f[6])  1.043      1.097
#> 11 sd(Days_f[7])  0.782      0.838
#> 12 sd(Days_f[8])  0.986      1.055
#> 13 sd(Days_f[9])  1.136      1.194
#> 14           cor  0.824      0.842

The day-specific SDs generally increase over the study, implying similar pattern to the one modeled by the random-slope model previously. HCS recovers the increase without assuming a linear trend per subject, and its day-specific SDs also capture departures from the smooth pattern implied by the random-slope model, such as the smaller SD on day 7.

Testing the Correlation

The structured covariance constructions allow us to also test the correlation value using a Savage-Dickey density ratio (Wagenmakers et al., 2010) and the product-space approach. Since the day-specific random effects are so strongly correlated that all posterior draws of such a test would include the correlation (as in the random-slope test), we test whether the HCS correlation equals 0.9 (against a hypothesis with a uniform prior over the admissible range).

First, the Savage-Dickey density ratio allows us to obtains the Bayes factor directly from the fitted HCS model: it divides the posterior density of the correlation at 0.9 by its prior density.

hypothesis_BF(fit_sleep_hcs, hypothesis = "(mu) cor = 0.9", BF01 = TRUE, seed = 1)
#>     Alternative:          Null:  1/BF error%(BF)
#>  (mu) cor != 0.9 (mu) cor = 0.9 5.218         NA

Second, a mixture prior places half of the prior probability on a correlation of 0.9 (the spike) and half on the same uniform distribution (the slab) and tests the hypothesis within a single fit using the product space approach. random_covariance(cor_scale = "cor") specifies that the correlation prior is on the raw correlation scale.

cor_test_prior <- prior_mixture(
  prior_list = list(
    prior(distribution = "spike",   parameters = list(location = 0.9)),
    prior(distribution = "uniform", parameters = list(a = -1/9, b = 1))
  ),is_null = c(TRUE, FALSE)
)
day_test_priors <- list(mu = prior_random(Subject = random_block(
  sd         = unit_sd_prior,
  covariance = random_covariance(cor = cor_test_prior, cor_scale = "cor"),
  monitor    = random_monitor(latent = FALSE, coefficients = FALSE, correlation = TRUE)
)))

fit_sleep_hcs_test <- JAGS_fit(
  model_syntax = sleep_likelihood, data = sleep_data, prior_list = sleep_prior_list,
  formula_list = sleep_hcs_formula_list, formula_data_list = list(mu = sleepstudy),
  formula_prior_list = sleep_fixed_priors, formula_random_prior_list = day_test_priors,
  formula_scale_list = list(mu = list(Days = TRUE)),
  chains = 2, adapt = 1000, burnin = 2000, sample = 5000, seed = 31
)

The JAGS_inference_table function compares the slab with the spike, so an inclusion Bayes factor below 1 favors a correlation of 0.9.

JAGS_inference_table(fit_sleep_hcs_test)
#>                   Prior prob. Post. prob. Inclusion BF
#> (mu) Subject: cor       0.500       0.194        0.241

Both approaches come to a similar conclusion. The data favor a correlation of 0.9 over the uniform alternative by a Bayes factor of about 5. The Savage-Dickey ratio reuses the estimated model but relies on a density estimate at the tested value; the mixture prior needs a separate fit but again allows us to extend the approach to testing several parameters at once.

Nested Random Effects

The Pastes data contain the strength of a chemical paste measured twice in each of three casks from each of ten delivery batches (60 measurements). The casks are nested within the batches, so a cask’s random effect has to be defined within its batch. We again standardize the response so the unit-scale priors from the previous examples apply.

data("Pastes", package = "lme4")
Pastes$strength_z <- as.numeric(scale(Pastes$strength))
head(Pastes)
#>   strength batch cask sample strength_z
#> 1     62.8     A    a    A:a 0.84822230
#> 2     62.6     A    a    A:a 0.78645854
#> 3     60.1     A    b    A:b 0.01441154
#> 4     62.3     A    b    A:b 0.69381290
#> 5     62.7     A    c    A:c 0.81734042
#> 6     63.1     A    c    A:c 0.94086794

In lme4, the slash notation (1 | batch/cask) expands to (1 | batch) + (1 | batch:cask): one random intercept for each batch and one for each cask within its batch.

lme4_nested <- lme4::lmer(
  strength_z ~ 1 + (1 | batch/cask),
  data = Pastes, REML = FALSE
)
lme4::VarCorr(lme4_nested)
#>  Groups     Name        Std.Dev.
#>  cask:batch (Intercept) 0.89683 
#>  batch      (Intercept) 0.33818 
#>  Residual               0.25428

BayesTools accepts the same slash notation and creates one random-effect block per nesting level, named batch and cask_batch. The likelihood, the data, and the intercept prior follow the previous examples, with strength_z as the response.

nested_likelihood <-
"model{
  for(i in 1:N){
    strength_z[i] ~ dnorm(mu[i], pow(sigma, -2))
  }
}
"

nested_formula <- list(
  mu = ~ 1 + (1 | batch/cask)
)

nested_data <- list(strength_z = Pastes$strength_z, N = nrow(Pastes))
nested_fixed_priors <- list(intercept = unit_prior)

Instead of separate SD priors for the two blocks, we specify a variance allocation prior with random_variance_allocation(). It places one prior on the total random-effect SD, sd_total and a Dirichlet prior on the proportions of the total variance assigned to each block, var_prop. With alpha = c(1, 1), every split between the batch and cask levels is equally likely a priori. This is convenient when the overall amount of heterogeneity is easier to anticipate than its split between the levels, and the prior distribution on the total SD stays the same however many blocks share it.

name = "nested_re" labels the allocation and prefixes its rows (nested_re: sd_total). terms lists the random-effect blocks it scales: the values "batch" and "cask_batch" are the blocks generated by (1 | batch/cask), and the names batch and cask label the components in the output (var_prop(batch), var_prop(cask)). Other specifications are possible: independent SD priors for each block, prior_random(batch = random_block(sd = ...), cask_batch = random_block(sd = ...)), as shown for crossed random effects below.

With only ten batches, the intercept and the batch effects are strongly correlated a posteriori, which slows down the default (noncentered) sampler. random_block(parameterization = "mean_centered") samples the batch effects around the fixed intercept instead, which removes this correlation without changing the model or its priors. The block takes its SD from the allocation, so it needs no sd prior.

nested_allocation_priors <- list(
  mu = prior_random(
    batch      = random_block(parameterization = "mean_centered"),
    allocation = random_variance_allocation(
      name    = "nested_re",
      terms   = c(batch = "batch", cask = "cask_batch"),
      sd      = unit_sd_prior,
      weights = prior("dirichlet", parameters = list(alpha = c(1, 1)))
    ),
    monitor = random_monitor(latent = FALSE, coefficients = FALSE, correlation = FALSE)
  )
)

fit_nested <- JAGS_fit(
  model_syntax = nested_likelihood, data = nested_data, prior_list = list(sigma = unit_sd_prior),
  formula_list = nested_formula, formula_data_list = list(mu = Pastes),
  formula_prior_list = list(mu = nested_fixed_priors),
  formula_random_prior_list = nested_allocation_priors,
  chains = 4, adapt = 1000, burnin = 1000, sample = 4000, seed = 41
)

JAGS_estimates_table(fit_nested, random_effects_summary = "standard")
#>                                   Mean    SD  0.025    0.5 0.975 error(MCMC) error(MCMC)/SD  ESS R-hat
#> (mu) intercept                  -0.002 0.242 -0.484 -0.001 0.479     0.00906          0.037  720 1.009
#> sigma                            0.266 0.036  0.208  0.262 0.346     0.00062          0.017 3300 1.001
#> (mu) nested_re: sd_total         1.055 0.172  0.792  1.029 1.466     0.00921          0.054  353 1.005
#> (mu) nested_re: var_prop(batch)  0.251 0.164  0.011  0.232 0.607     0.00912          0.056  359 1.008
#> (mu) nested_re: var_prop(cask)   0.749 0.164  0.393  0.768 0.989     0.00912          0.056  359 1.008

The nested lme4 and BayesTools estimates are compared in the table below.

compare_nested_estimates(lme4_model = lme4_nested, bayes_model = fit_nested)
#>                    parameter  lme4 BayesTools
#> 1                  intercept 0.000     -0.002
#> 2                      sigma 0.254      0.266
#> 3       batch: sd(intercept) 0.338      0.505
#> 4  cask_batch: sd(intercept) 0.897      0.899
#> 5        nested_re: sd_total 0.958      1.055
#> 6 nested_re: var_prop(batch) 0.124      0.251
#> 7  nested_re: var_prop(cask) 0.876      0.749

Crossed Random Effects

The Penicillin data contain the diameter of the growth-inhibition zone of six penicillin samples, each tested once on each of 24 agar plates (144 measurements). Every sample appears on every plate, so the plate and sample effects are crossed: a measurement shares its plate with the other five samples and its sample with the other 23 plates. We again standardize the response so the unit-scale priors from the previous examples apply.

data("Penicillin", package = "lme4")
Penicillin$diameter_z <- as.numeric(scale(Penicillin$diameter))
head(Penicillin)
#>   diameter plate sample  diameter_z
#> 1       27     a      A  1.98311734
#> 2       23     a      B  0.01367667
#> 3       26     a      C  1.49075717
#> 4       23     a      D  0.01367667
#> 5       23     a      E  0.01367667
#> 6       21     a      F -0.97104366

In lme4, crossed random effects are separate terms, (1 | plate) + (1 | sample): one random intercept for each plate and one for each sample.

lme4_crossed <- lme4::lmer(
  diameter_z ~ 1 + (1 | plate) + (1 | sample),
  data = Penicillin, REML = FALSE
)
lme4::VarCorr(lme4_crossed)
#>  Groups   Name        Std.Dev.
#>  plate    (Intercept) 0.41633 
#>  sample   (Intercept) 0.87180 
#>  Residual             0.27076

BayesTools uses the same formula; the two random-effect blocks take the names of their grouping variables, plate and sample, which the priors below refer to. The likelihood, the data, and the intercept prior follow the previous examples, with diameter_z as the response.

crossed_likelihood <-
"model{
  for(i in 1:N){
    diameter_z[i] ~ dnorm(mu[i], pow(sigma, -2))
  }
}
"

crossed_formula <- list(mu = ~ 1 + (1 | plate) + (1 | sample))
crossed_data    <- list(diameter_z = Penicillin$diameter_z, N = nrow(Penicillin))
crossed_fixed_priors <- list(intercept = unit_prior)

Independent SD Priors

The most direct specification gives each block its own half-normal SD prior through random_block(), addressed by the block names from the formula. As in the nested example, the sample block is mean-centered: its six effects are strongly correlated with the intercept a posteriori.

crossed_independent_priors <- list(
  mu = prior_random(
    plate  = random_block(sd = unit_sd_prior, monitor = random_monitor(latent = FALSE)),
    sample = random_block(sd = unit_sd_prior, monitor = random_monitor(latent = FALSE),
                          parameterization = "mean_centered")
  )
)

fit_crossed_independent <- JAGS_fit(
  model_syntax = crossed_likelihood, data = crossed_data, prior_list = list(sigma = unit_sd_prior),
  formula_list = crossed_formula, formula_data_list = list(mu = Penicillin),
  formula_prior_list = list(mu = crossed_fixed_priors),
  formula_random_prior_list = crossed_independent_priors,
  chains = 4, adapt = 2000, burnin = 2000, sample = 4000, seed = 2
)

JAGS_estimates_table(fit_crossed_independent, random_effects_summary = "standard")
#>                             Mean    SD  0.025   0.5 0.975 error(MCMC) error(MCMC)/SD  ESS R-hat
#> (mu) intercept             0.010 0.396 -0.780 0.011 0.808     0.00400          0.010 9946 1.000
#> sigma                      0.274 0.018  0.241 0.273 0.313     0.00022          0.012 6905 1.000
#> (mu) plate: sd(intercept)  0.441 0.076  0.319 0.431 0.619     0.00267          0.035  870 1.002
#> (mu) sample: sd(intercept) 1.043 0.318  0.592 0.985 1.817     0.00431          0.014 5581 1.001

Variance Allocation Prior

Alternatively, as in the nested example, random_variance_allocation() places one prior on the total SD of the two blocks, total_re: sd_total, and a uniform Dirichlet prior on the proportions of its variance assigned to plate and sample. The sample block keeps the mean-centered parameterization while its SD now comes from the allocation specification.

crossed_allocation_priors <- list(
  mu = prior_random(
    sample     = random_block(parameterization = "mean_centered"),
    allocation = random_variance_allocation(
      name    = "total_re",
      terms   = c(plate = "plate", sample = "sample"),
      sd      = unit_sd_prior,
      weights = prior("dirichlet", parameters = list(alpha = c(1, 1)))
    ),
    monitor = random_monitor(latent = FALSE)
  )
)

fit_crossed_allocation <- JAGS_fit(
  model_syntax = crossed_likelihood, data = crossed_data, prior_list = list(sigma = unit_sd_prior),
  formula_list = crossed_formula, formula_data_list = list(mu = Penicillin),
  formula_prior_list = list(mu = crossed_fixed_priors),
  formula_random_prior_list = crossed_allocation_priors,
  chains = 4, adapt = 2000, burnin = 2000, sample = 4000, seed = 3
)

JAGS_estimates_table(fit_crossed_allocation, random_effects_summary = "standard")
#>                                   Mean    SD  0.025    0.5 0.975 error(MCMC) error(MCMC)/SD  ESS R-hat
#> (mu) intercept                  -0.001 0.365 -0.730 -0.001 0.727     0.00443          0.012 7070 1.001
#> sigma                            0.274 0.018  0.240  0.273 0.313     0.00022          0.012 6757 1.000
#> (mu) total_re: sd_total          1.023 0.233  0.701  0.980 1.627     0.01226          0.053  390 1.029
#> (mu) total_re: var_prop(plate)   0.216 0.096  0.064  0.204 0.432     0.00452          0.047  459 1.014
#> (mu) total_re: var_prop(sample)  0.784 0.096  0.568  0.796 0.936     0.00452          0.047  459 1.014

The two BayesTools fits differ only in the prior on the random-effect SDs: the independent priors treat the two blocks separately, while the allocation prior sets the scale of their total and estimates how it splits. The table below compares both with lme4.

compare_crossed_estimates(
  lme4_model = lme4_crossed,
  bayes_independent = fit_crossed_independent, 
  bayes_allocation = fit_crossed_allocation)
#>               parameter  lme4 BayesTools_independent BayesTools_allocation
#> 1             intercept 0.000                  0.010                -0.001
#> 2                 sigma 0.271                  0.274                 0.274
#> 3  plate: sd(intercept) 0.416                  0.441                 0.447
#> 4 sample: sd(intercept) 0.872                  1.043                 0.912

References

Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1–48. https://doi.org/10.18637/jss.v067.i01

Bürkner, P.-C. (2017). brms: An R package for Bayesian multilevel models using Stan. Journal of Statistical Software, 80(1), 1–28. https://doi.org/10.18637/jss.v080.i01

Gronau, Q. F., Singmann, H., & Wagenmakers, E.-J. (2020). bridgesampling: An R package for estimating normalizing constants. Journal of Statistical Software, 92(10), 1–29. https://doi.org/10.18637/jss.v092.i10

Kuo, L., & Mallick, B. (1998). Variable selection for regression models. Sankhyā: The Indian Journal of Statistics, Series B, 65–81.

Lodewyckx, T., Kim, W., Lee, M. D., Tuerlinckx, F., Kuppens, P., & Wagenmakers, E.-J. (2011). A tutorial on Bayes factor estimation with the product space method. Journal of Mathematical Psychology, 55(5), 331–347. https://doi.org/10.1016/j.jmp.2011.06.001

Plummer, M. (2003). JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling. Proceedings of the 3rd International Workshop on Distributed Statistical Computing, 124, 1–10.

Wagenmakers, E.-J., Lodewyckx, T., Kuriyal, H., & Grasman, R. (2010). Bayesian hypothesis testing for psychologists: A tutorial on the Savage-Dickey method. Cognitive Psychology, 60(3), 158–189. https://doi.org/10.1016/j.cogpsych.2009.12.001

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.