# How to estimate multivariate posterior probability?

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"])
 0.9555
sum(pe[,"back_sexUNK"] > 0)/length(pe[,"back_sexUNK"])
 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?

``````tarsus_pos = pe[,"tarsus_sexUNK"] > 0
back_pos = pe[,"back_sexUNK"] > 0
n_samples = nrow(pe)
p_tarsus_back_pos = sum(tarsus_pos & back_pos)/n_samples
``````
2 Likes

Thanks a lot, @Guido_Biele!