How to calculate VPC in brms

Interfaces brms
,
1

Hello! I run the binary multilevel model using brms . I should calculate the variance partition coefficient (VPC) to represent the percentage variance explained at the group level (\sigma_μ^2). I first set up a null model with only random intercept. The code is as below.

fm.obj.v00<-bf(Relief3 ~ 1 + (1 | Group))
fit.obj.v00<-brm(formula = fm.obj.v00, data = orgin, family=bernoulli(link = "logit"), 
  control = list(adapt_delta = 0.99, max_treedepth=15), 
   iter=8000, warmup=4000, thin=4, chains=1, cores=1)
fit.obj.v00<-add_criterion(fit.obj.v00, c("loo", "waic"))

I then would like to calculate VPC based on brms output. However, I don’t know how to calculate the variance at the individual level (\sigma_ε^2 )and that at the group level (\sigma_u^2). Below is the equation for VPC. Just wondering if it’s possible to provide the code to calculate VPC based on brms output.

VPC = \sigma_u^2 / (\sigma_u^2 + \sigma_ε^2) = \sigma_u^2 / (\sigma_u^2 + \pi^2 / 3)

I have used the code as below to calculate the mean and variance of posterior samples. However, I don’t think r_mean represents the variance of residuals at the group level. If so, how can I calculate VPC based on the output of brms?

vars <- posterior_samples(fit.obj.v00, pars = "^sd_")^2
r_mean<-apply(vars, 2, mean)
r_variance<-apply(vars, 2, sd)

I have also review one topic “extracting variance components” and it shows codes as below. Can I replace “post” with “fit.obj.v00”? How to interprate the output in summary()? Does the calculation results represent the variance components or proportion of variance at the group level?

PPD <- posterior_predict(post)
vars <- apply(PPD, MARGIN = 1, FUN = var)
PPD_0 <- posterior_predict(post, re.form = ~ 0)
vars_0 <- apply(PPD_0, MARGIN = 1, FUN = var)
summary(vars - vars_0)

Another way is to use “performance” package and adopt icc() function. I run it and the code and results are as below. Just wonder if it’s right for the output of brms.

library("performance")
icc(fit.obj.v00, by_group="Group")
 #Intraclass Correlation Coefficient
 #Adjusted ICC: 0.449
 #Conditional ICC: 0.449

P.S.: Please find the relevant documents on VPC in http://seis.bris.ac.uk/~frwjb/materials/wbvpc.pdf.

Best,

MODERATOR EDIT: This post was edited for code and math formatting by @martinmodrak

1 Like
2

Hi,
generally, in the Bayesian context, you would not get a single estimate for any quantity, instead you seek a posterior distribution of a quantity. But this is quite simple to achieve - for each posterior sample of \sigma_\mu you calculate a single value of VPC. This gives you a posterior distribution for VPC, which you can then summarise as you see fit (e.g. by mean + 95% credible interval). I’ve never used VPC myself, but if I understand the formulas right, you could compute the posterior distribution as

vars <- posterior_samples(fit.obj.v00, pars = "sd_Group")^2
VPC_samples <- vars/(vars + pi^2 / 3)

Does that make sense?

Best of luck with your model!

P.S.: Note that you can use triple backticks (```) to format code blocks and that Latex-style math is supported within segments marked by dollar signs… I’ve edited your post to add this formatting.

1 Like