I’d say the problem is rather different if you want to implement predictors of the group-level parameters, or use the group-level parameters as predictors of something else. In the case where you want to implement predictors of the group-level parameters; you could view this as a nonlinear modelling problem where the intercept and slope are nonlinear parameters, and they are described by linear submodels that implement the between-subject variation and the predictors. This is a typical form of the problem as applied in pharmacometrics for example. It corresponds to the group-level values including both systematic and random variation.

I’m not sure if there is a way to specify this in the typical linear brms syntax, but you could do it with the nonlinear syntax.

Here’s a toy example (excuse the potentially poor priors, they would need to be set with care in the real case):

```
example_nlin_data <- data.frame(ID = rep(seq(1,20), each = 10), x1 = runif(200,-5,5), x2 = rep(runif(10,-5,5), each = 20))
example_nlin_data$ID_intercept <- rep(rnorm(20,0,1), each = 10)
example_nlin_data$ID_slope <- 1+(example_nlin_data$x2*2)
example_nlin_data$y1 <- example_nlin_data$ID_intercept+(example_nlin_data$x1*example_nlin_data$ID_slope) + rnorm(200,0,5)
ggplot(data = example_nlin_data, aes(x = x1, y = y1)) + geom_point(alpha = 0.6) + facet_wrap(~ID, ncol = 4) + theme_bw()
example_nlin_mod1 <- brm(bf(y1 ~ nlinintercept + (x1*nlinslope), nlinintercept ~ 1 + (1|ID), nlinslope ~ 1 + x2 + (1|ID), nl = TRUE), data = example_nlin_data, control = list(adapt_delta = 0.95), prior = set_prior('normal(0,5)', class = 'b', nlpar = 'nlinintercept') + set_prior('normal(0,5)', class = 'b', nlpar = 'nlinslope'))
summary(example_nlin_mod1)
```

So in this example the subject-level slope in the linear relationship between `y1`

and `x1`

is itself dependent on a linear function of the predictor `x2`

(which has one value per subject), with intercept=1 and slope=2.

Those are reasonably well-recovered:

```
Family: gaussian
Links: mu = identity; sigma = identity
Formula: y1 ~ nlinintercept + (x1 * nlinslope)
nlinintercept ~ 1 + (1 | ID)
nlinslope ~ 1 + x2 + (1 | ID)
Data: example_nlin_data (Number of observations: 200)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Group-Level Effects:
~ID (Number of levels: 20)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(nlinintercept_Intercept) 0.78 0.50 0.04 1.88 1.00 1604 1989
sd(nlinslope_Intercept) 0.19 0.14 0.01 0.52 1.00 1527 2055
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
nlinintercept_Intercept 0.15 0.40 -0.62 0.93 1.00 5021 2964
nlinslope_Intercept 1.05 0.13 0.79 1.32 1.00 5581 2738
nlinslope_x2 2.00 0.05 1.91 2.09 1.00 5112 2912
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 4.96 0.27 4.48 5.53 1.00 6341 2784
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).
```

For the case where you want to use the subject-level parameters as predictors for something else, I think that could be viewed as a joint modeling problem (might not be very workable in brms, but readily implemented in Stan).