Let me use the example code about multivariate modeling in Paul Bürkner’s blog post to demonstrate my question:

```
library(MCMCglmm)
library(brms)
data("BTdata", package = "MCMCglmm")
fit1 <- brm(
mvbind(tarsus, back) ~ sex + hatchdate + (1|p|fosternest) + (1|q|dam),
data = BTdata, chains = 2, cores = 2
)
```

I can obtain the posterior probability of an effect being positive for each of the two response variables with the following

```
pe <- fixef(fit1, summary = FALSE)
sum(pe[,"tarsus_sexUNK"] > 0)/length(pe[,"tarsus_sexUNK"])
[1] 0.9555
sum(pe[,"back_sexUNK"] > 0)/length(pe[,"back_sexUNK"])
[1] 0.834
```

Suppose that I want to estimate the posterior probability of the effects for both response variables being positive (e.g., `pe[,"tarsus_sexUNK"]`

and `pe[,"back_sexUNK"]`

). How can I achieve this?