Allowing groups to have different inter-individual variances in hierarchical model

Modeling
1

I am modeling a dataset where each of several hundred individuals have ranked the same 10 objects from best to worst. Each individual is the member of a different group. The goal of the model is to estimate the underlying valence (desirability) of each object, taken to be continuous. I am not showing the part that relates rankings to continuous values, but the hierarchical structure of the model looks like:

data {
    int<lower=1> n_objects; // e.g. 10
    int<lower=1> n_individuals; // e.g. ~300
    int<lower=1, upper=n_individuals> n_groups; // e.g. 9
    int<lower=1, upper=n_groups> group_id[n_individuals]; // The integer group IDs of each individual
    int<lower=1, upper=n_objects> ranks[n_individuals, n_objects]; // The matrix of integer ranks
    }
    
parameters {
    vector[n_objects] mu_global;
    matrix[n_groups, n_objects] mu_group;
    matrix[n_individuals, n_objects] mu_ind;
    real<lower=0.0001> sigma_global;
    real<lower=0.0001> sigma_group;
    real<lower=0.0001> sigma_ind[n_groups];
    }
   
transformed parameters {
    real mu_mean = mean(mu_global); // Used for centering
}
    
model {
    mu_mean ~ normal(0, 0.001); // Used for centering
    mu_global ~ normal(0, sigma_global);
    sigma_global ~ cauchy(0, 3);
    sigma_group ~ cauchy(0, 3);
    sigma_ind ~ cauchy(0, 3);
    for (i in 1:n_groups) {
        mu_group[i] ~ normal(mu_global, sigma_group);
    }
    for (i in 1:n_individuals) {
        mu_ind[i] ~ normal(mu_group[group_id[i]], sigma_ind[group_id[i]]);
        target += some_func(mu_ind[i], ranks[i]); // A function not shown here
        }
    }

This model fits pretty well, R_hat ~ 1, etc., and I get reasonable looking histograms for all parameters and point estimates that make sense given the data. I know from the data that some groups have much higher inter-individual variance than others in how they rank the objects, but allowing sigma_ind to vary by group as I do (sigma_ind[n_groups]) makes the model much slower to fit, and in general I am surprised at how slow sampling proceeds.

Is there a better strategy for this kind of hierarchical model? I also recognize that my model may be misparameterized in other ways, so any suggestions there would also be appreciated.

2

Hi, sorry that your question slipped through.

In general this is hard to diagnose, but I would look at pair plots of some representative parameters to look for any correlations/funnels/banana-like shapes, those might provide a hint.

As a pure guess, I think the problem might be that you just don’t have enough data to reliably estimate sigma per-group. This then can show up as problems described at https://betanalpha.github.io/assets/case_studies/underdetermined_linear_regression.html.

You usually need more data to estimate variance than to estimate mean (and even more to estimate higher moments). An intermediate step might be to force some hierarchical structure on the sigmas, e.g. to have: \log \sigma \sim N(\nu, \tau) with \nu and \tau as additional paremeters.

Few technical notes:

Your approach to centering is a bit unusual (but I don’t understand your model deep enough to be sure). If each element of mu_global has mean 0, than the their overall mean is also zero, so not sure what you are trying to achieve.

The cauchy(0, 3) priors are very wide and usually allow unrealistic model configurations. This might in turn hinder sampling. For most models it is usually OK to restrict to something like normal(0, 1), but you need to understand what a reasonable prior might be your model. When in doubt, you can do prior predictive checks to see if your priors are reasonable.

Best of luck with your model.