Simple intercept only hierarchical model with two groups: deadly slow, poor convergence

Modeling
, ,
4

Hierarchical models in Stan (and Bayes generally) can have difficulties with sampling and convergence when specified directly, instead it’s recommended to use a non-centered parameterisation which constructs the parameter as a function of a standard normal. There’s more background for this in the User’s Guide.

For parameters which don’t have a lower bound, you can implement this very simply using the offset and multiplier syntax without changing anything else in the model:

real<offset=prior_mu_grand,multiplier=prior_sigma_grand> mu_grand;

However for parameters with a lower bound (e.g., sigma_sp), you need to use the ‘manual’ approach which specifies a standard normal parameter with an appropriate lower bound:

data {
    real prior_sigma_sp_mu;
    real prior_sigma_sp_sigma;
}
parameters {                 
    real<lower=-prior_sigma_sp_mu/prior_sigma_sp_sigma> sigma_sp_raw;
}
transformed parameters {
  // Implies sigma_sp<lower=0> ~ normal(prior_sigma_sp_mu, prior_sigma_sp_sigma)
  real sigma_sp = prior_sigma_sp_mu + sigma_sp_raw * prior_sigma_sp_sigma;
}
model {
    sigma_sp_raw ~ std_normal();
}

See this post for more background on how these constraints are derived.

Additionally, with stan models in general you will want to use:

vector[N] parameter;

Instead of:

real parameter[N];

As it is more efficient to work with vectors than arrays of real. This also allows you to vectorise your likelihood, so you can instead write:

yTraiti ~ normal(mu_grand + muSp[species] + muStudy[study], sigmaTrait_y);

Without the need for the loop or the ymu variable