Multilevel model with partially pooled variance

Question

I have a model with N data points and J groups, where I want to partially pool the means and the variances of the groups.

Within group j, for data $\{y_i\}$ I'm assuming

$y[i] \sim \mathcal{N}(\theta[j], \sigma_y[j])$

For $\theta$, I have

$\theta[j] \sim \mathcal{N}(\mu_\theta, \sigma_\theta)$ for j = 1, ..., J

I'm not sure what prior to use for the $\sigma_y$ parameter. I've tried inverse gamma; i.e,

$\sigma_y^2[j] \sim IG(a, b)$ for j = 1, ..., J

but my chains aren't converging for a toy data set I've created. Even when I pass the "true" parameter values as the initial values for the chain, it still doesn't work.

How do people typically handle models like this? This isn't working in stan or in pymc, so I'm thinking there's something wrong with my probability model.

Answer 1

The general idea of your model can be handled in Stan (and possibly in PyMC and other statistical inference packages). Here are a couple of suggestions:

Hyperparameters

We (i.e. Andrew) have been suggesting a half Cauchy prior.

Reparameterization

The non-centered parameterization should help this model. The posterior geometry of $theta$ and $\sigma_y$ will look like a funnel. This is hard to sample from. Reparameterizing the model should help. See Section 20.1 in the Stan 2.8.0 manual for details.

Note: this reparameterization is dependent on the shape of the posterior distribution. The posterior distribution changes with different data, so thinking that one way to write the model will be the best way for all data isn't true.