Suppose that I have the following model for each grouping unit k (e.g., province):

y_{ijk} \sim N (a_{ik} + \beta_{jk}, \sigma_k^2), k=1,2,..., K

where a_{ik} are the population-level effects while the group-level effects are \beta_{jk} \sim N(0, \tau_k^2). I’m interested in the variance ratio under the above model,

r_k = \tau_k^2/( \tau_k^2 + \sigma_k^2)

However, the above computations are done for each k separately, leading to a multiplicity problem. So, I would like to compute the variance ratio for each of those entities by adopting partial pooling through incorporating all the units into one Bayesian multilevel model with index,

with the assumptions \beta_j \sim N(0, \tau^2) and \gamma_k \sim N(0, \eta^2).

However, I’m not so sure if the above model is reasonable. In addition, I cannot conceptually figure out a way to formulate the variance ratio for each k with the posterior samples from the Bayesian model. Any suggestions?

Here’s an outline of how to do it. Check the documentation for a discussion of centering options. This also isn’t vectorized. A fully vectorized version might run substantially faster. I also haven’t debugged it, so I can’t guarantee the syntax will work, but the idea should be right.

model {
for (i in 1:I) {
for (j in 1:J) {
for (k in 1:K) {
y[i,j,k] = normal(mu[i,j,k], S_e[k])
}
}
}
}

You’ll need to define mu[i,j,k] in the transformed parameters block and set appropriate priors for a[i], b[j], c[k], and S_e[k].

You can predict sigma with brms via brms distributional syntax: bf(y ~ ..., sigma ~ ...). For more details see ?brmsformula and the brms_distributional and brms_multilevel vignettes.