I’m trying to use brm() to assess the interaction between a non-linear smooth and the 2-way interaction between two binary factors. I think that such 3-way interactions (between a continuous and more than one categorical predictors) are not possible in mgcv. Are the possible in brms? If so, how would I specify the formula for something that in effect would be ~ s(x1) * x2 * x3 (but I don’t think that works?)

As brms uses mgcv under the hood, it shares the same limitations. The workaround I am using myself for interactions with categorical variables, is a simple yet (I believe) effective one:

data$x2x3 <- paste0(data$x2, "_", data$x3)
be ~ s(x1, by = x2x3)

The framework for s(x1) * x2 * x3 to make sense it actually already implemented in brms, but splines are not (yet?) part of it. Maybe at some point, this syntax might actually work.

Thank you for the quick clarification. I considered that option, but I didn’t see how you can then assess the interaction (as opposed to the joint consequence of the main effects and interaction)?

I guess one can compare (using loo or compare_ic?) to a model with just x3 by x1 and x3 by x2, in order to know how much the interaction matters?

And if I want to know for which values of x3 , say, a_1 - a_2 > b_1-b_2 (where a, b would be the values x2, and 1, 2 would be the values of x3 in the example) meets a certain evidence criterion, can I do that through hypothesis ()? Is there perhaps even an already implemented way of plotting a_1 - a_2 > b_1-b_2 against x3 along with uncertainty intervals?