Multiple-outcome causal model in brms

I would like to fit a multivariate causal model in brms, according to the following hypothesized DAG:


Through the model, I want to understand the total effect of the treatment a on both outcomes y1 and y2, while taking into account a mediator m and another covariate b. I further want to understand the correlation between the outcomes y1 and y2.

This was my attempt so far:

model <- brm(
  data = data,
  bf(y1 ~ 1 + a + b + m) +
    bf(y2 ~ 1 + a + b + m) +
    bf(m ~ 1 + a) +

This estimates the residual correlation between y1 and y2, but also between y1 and m and y2 and me, which I am not interested in. Is it possible to exclude those correlations somehow from the covariance matrix?

I am not an expert in causal modelling, so my problem may be ill-posed. Any thoughts on how to approach this problem are greatly appreciated.

  • Operating System: Windows 10
  • brms Version: 2.17.0