I believe you linked to an older post. Is this the post you meant to link?

I recently looked into this and there’s an easier way without overwriting the likelihood. It depends on the fact that brms models random effects with a non-centralized parametrization. E.g. if this is our model:

` y ~ 1 + (1|ID)`

the underlying model is

` y = b_0 + eps_resid`

` b_0 ~ a_0 + sigma_id * Norm(0, 1)`

`eps_resid ~ Norm(0,sigma)`

then you can use non-linear formulas to:

- fix sigma_id to 1
- multiply by a new parameter mysigma, which can then be predicted by any linear formula.

For example.

## Setup

```
library(cmdstanr)
library(brms)
library(tidyverse)
library(knitr)
library(bayesplot)
options(mc.cores = parallel::detectCores())
dat <- rtdists::speed_acc %>%
filter(block %in% c(1,2,3,4,5))
```

## Basic model

```
f1 <- bf(rt ~ stim_cat + (1|id))
fit1 <- brm(f1, data = dat, chains = 4)
```

## Custom model to predict random effect sigma

```
# non-linear formula
f2 <- bf(rt ~ murt + idintercept * exp(logmysigma), # specify explicit terms in a non-linear formula
murt ~ stim_cat, # predict rt mean via murt
idintercept ~ 1 + (1|id), # specify random effects (will fix via prior)
logmysigma ~ response, # model random effects variance (response is a categorical predictor here, but can be anything)
nl = T)
# set priors to fix parameters, and to avoid flat priors on nlp parameters
prior2 <- prior(student_t(3,0.6,2.5), coef='Intercept', nlpar="murt") +
prior(constant(0), class='b', nlpar='idintercept') + # center random effects at 0
prior(constant(1), class='sd', nlpar='idintercept') + # fix random effects sd to 1
prior(student_t(3,0,2.5), class='b', nlpar="logmysigma") # priors on our sigma
fit2 <- brm(f2, data = dat, chains = 4, iter = 4000, prior=prior2)
```

You get the following output:

```
Family: gaussian
Links: mu = identity; sigma = identity
Formula: rt ~ murt + idintercept * exp(logmysigma)
murt ~ stim_cat
idintercept ~ 1 + (1 | id)
logmysigma ~ response
Data: dat (Number of observations: 8160)
Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
total post-warmup draws = 8000
Multilevel Hyperparameters:
~id (Number of levels: 17)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(idintercept_Intercept) 1.00 0.00 1.00 1.00 NA NA NA
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
murt_Intercept 0.61 0.02 0.57 0.65 1.00 791 1831
murt_stim_catnonword 0.03 0.01 0.01 0.05 1.00 1742 3484
logmysigma_Intercept -2.52 0.21 -2.90 -2.09 1.00 950 1847
logmysigma_responsenonword 0.29 0.08 0.12 0.45 1.00 3830 5070
logmysigma_responseerror -1.73 2.88 -9.56 1.37 1.00 2638 1289
idintercept_Intercept 0.00 0.00 0.00 0.00 NA NA NA
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.30 0.00 0.30 0.31 1.00 4650 4822
```

In this model, the idintercept predictor is a random effects intercept with mean 0 and sd 1. We multiply it by exp(logmysigma), thus results in exp(logmysigma) as the standard deviation of the random intercept. From the results we see that indeed the random effect variance differs between conditions. You can specify any formula on logmysigma.

If I understood your post correctly, this achieves the same. Still hacky, but it does not require custom likelihoods and works for any model.