Computing variance ratio

Modeling
1

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,

y_{ijk} \sim N(a_i + \beta_j + \gamma_k, \sigma^2)

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?

2

I think you’d want S_e^2 to vary by province as well as C_k, that way you can compute a separate ratio for each from the posterior.

3

Thanks, Mike!

Do you mean to add an interaction term (\beta\eta)_{jk} in the model?

y_{ijk} \sim N( a_i + \beta_j + \eta_k + (\beta\eta)_{jk}, \sigma^2)

But how can I keep track of \sigma^2 for each k?

4

I think what he’s suggesting (or at least what I’d suggest) is

y_ijk ~ N(a_i + b_j + c_k, S_e^2_k)

Kent

5

Thanks, Kent! But how to exactly implement a model like y_{ijk} \sim N( a_i + \beta_j + \eta_k, \sigma_k^2)? With the R notation as adopted in brms,

y ~ a + (1 | b) + (1 | c)

how can I get a varying \sigma_k^2 at each k?

6

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].

Kent

1 Like
7

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.

1 Like