How to estimate multivariate posterior probability?

Modeling
#1

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?

#2 
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
#3

Thanks a lot, @Guido_Biele!