Allowing interaction of smooth by categorical predictor to vary across levels of random effect

I am trying to figure out the correct specification of random effects (specifically an interaction varying across levels of a random effect) when moving from a linear mixed model to a GAMM. I am working on brms which uses lme4/mgcv syntax respectively.

Suppose we have experimental plots in a 15x2X10 design, where each plot has 1-15 plant species, and is assigned to one of 2 levels of a categorical predictor (control vs treatment). Then these 30 combinations are replicated across 10 different locations.

We want to test if the relationship between the response and the "continuous" predictor differs between the two treatments and since the different locations may have different conditions we want to allow all that (intercepts, slopes, interaction) to vary by location (group).

For a linear mixed model the formula, using lme4 syntax, would be (I think) that of the 4th formula here:

y ~ 1 + cont + cat + cat:cont + (1 + cont + cat + cat:cont|group)


Now for a GAMM, if I understand correctly the information here and here, the formula, using mgcv syntax and excluding the cat:cont|group bit would be:

y ~ 1
+ s(cont, by = cat)
+ cat
+ s(group, bs='re')
+ s(cont, group, bs='re')
+ s(cat, group, bs='re')


The main difference is that random slopes and intercepts get individual terms and are uncorrelated (a limitation that I am not sure it is true for brms as well)

So where I am stuck is how to allow the s(cont, by = cat) term to vary by group (which I think would be the equivalent of the cat:cont|group part of the random effect term in the lmm). I have tried adding the term s(interaction(cat, cont), group, bs='re')) but that throws an error.

Here is a reproducible example:

library(tidyverse)
library(brms)

d = data.frame(y = NA,
cat = as.factor(rep(rep(letters[1:2], 15), 10)),
cont = rep(rep(1:15, each = 2), 10),
group = as.factor(rep(1:10, each = 30)))

# linear relationship
d$$y[d$$cat == "a"] = 25 + 4*d$$cont[d$$cat == "a"]
d$$y[d$$cat == "b"] = 25 + 2*d$$cont[d$$cat == "b"]
d$$y = d$$y + rnorm(300, 0, 10)

d %>% ggplot(aes(cont, y, color = cat)) +
#geom_point() +
geom_smooth(method = "lm")

m = brm(bf(y ~ 1
+ cont
+ cat
+ cat:cont
+ (1
+ cont
+ cat
+ cat:cont|group)),
family = gaussian(),
chains = 3,
iter = 2000,
warmup = 500,
cores = 3,
#backend = "cmdstanr",
data = d)
summary(m)
plot(conditional_effects(m, "cont:cat", re_formula = NULL))

# non-linear relationship
d$$y[d$$cat == "a"] = 25 + 4*d$$cont[d$$cat == "a"] -.13*d$$cont[d$$cat == "a"]^2
d$$y[d$$cat == "b"] = 25 + 2*d$$cont[d$$cat == "b"] -.06*d$$cont[d$$cat == "b"]^2
d$$y = d$$y + rnorm(300, 0, 10)

d %>% ggplot(aes(cont, y, color = cat)) +
#geom_point() +
geom_smooth()

m = brm(bf(y ~ 1
+ s(cont, by = cat)
+ cat
+ s(group, bs='re')
+ s(cont, group, bs='re')
+ s(cat, group, bs='re')
# + s(interaction(cat, cont), group, bs='re') # this does not work
),
family = gaussian(),
chains = 3,
iter = 2000,
warmup = 500,
cores = 3,
#backend = "cmdstanr",
data = d)
summary(m)
plot(conditional_effects(m, "cont:cat", re_formula = NULL))


s(cont, group, bs = "fs", by = cat)

Where you have a random smooth of cont in the groups, allowing different smoothing parameters for each level of cat.
Re: s(cont, group, bs = "re"), yes, these are us correlated ranefs, but as for brms, you shouldn't use the ranef smooth to include random effects. Using (cont | group) would get correlated ranefs I believe.