I just would like to know whether there is a way to fit an interaction term between two response variables in a multivariate model in brms?
model_1 <- brm(
mvbind(y1, y2, y1:y2) ~ 1 + x1 + x2 + (1 |p| g1) + (1 |q| g2),
family = “gaussian”, chains = 4, cores = 4,
iter = 2000, warmup = 500, thin = 1