Possible identifiability issue in hierarchical model

I'm trying to fit some data using a hierarchical normal model

$y_i \sim N(\theta_i,\sigma^2)$

$\theta_i \sim N(\mu, \sigma_\theta^2)$

$(\mu,\sigma^2,\sigma_\theta^2) \sim diffuse$

I fit this model and I'm getting posteriors for $\sigma_\theta^2$ and $\sigma^2$ that are nearly identical. Is this an identifiability issue or a coincidence? There's no other information in the data to use to determine where the variability is coming from. Is there any way to still use this type of model, or is it just not useful without more data?

• It help to mention how much data you have. – John Salvatier May 9 '11 at 17:20

Your notation is a little strange (what do you mean by "diffuse"?), but I suspect that your prior on $\sigma^2_\theta$ is leading to an improper or nearly improper posterior, for one thing. See here for a detailed exposition of just this model and appropriate prior specification.

In short, yes, this model can be very useful and there probably ought to be some information about the variance parameters even in relatively small samples - but you need to be careful in how you specify and fit it.

Edit: When I wrote this answer I apparently hadn't read the OP properly (see my comment to @probabilityislogic's answer). Anyway as this model is written the parameters $\sigma, \sigma_\theta$ aren't separately identifiable as @probabilityislogic points out. I suspect that if you looked at the posterior distribution of $\sigma^2 + \sigma_\theta^2$ it would be doing something much more reasonable, and if you looked at the joint posterior of $\sigma, \sigma_\theta$ there would be a strong negative correlation.

You should go back to the original problem and try to reformulate this model - either it's not posed correctly in the OP or you're hosed, I think.

• Thanks. I had problems using improper priors (as you anticipated), but I got the results above using an inverse-gamma(1,1) prior on both the variance parameters. – user4528 May 9 '11 at 4:37
• Try to use a Uniform with a good range, instead of inverse-gamma. That's what Gelman suggests in his book with Hill. Besides that, if you want to chek your results, check if it makes sense. Do you expect, based on you data and knowledge, that variation among individuals is roughly the same as variation within individuals? – Manoel Galdino May 9 '11 at 14:15
• As @Manoel indicates, the conjugate prior isn't necessarily a good choice. Try another prior and see what happens. The paper above has BUGS/R code for the half-t prior IIRC. – JMS May 9 '11 at 16:04
• And another suggestion: In general, if you're wondering whether your model/code can perform like you hope it helps to generate some data from your prior and fit it. So here you might take $\sigma^2$ and $\sigma^2_\theta$ orders of magnitude different, generate a bunch of data and see what happens. – JMS May 9 '11 at 16:08
• One useful non-informative prior for a variance is to write $\sigma_{\theta}^{2}=g\sigma^{2}$ and then specify an F-distribution for g, with standard jeffreys prior for $p(\sigma^{2})\propto\frac{1}{\sigma^{2}}$. – probabilityislogic May 22 '11 at 4:04

Because you are dealing with normal-normal model, its not to hard to work out analytically whats going on. Now the standard argument for "diffuse" priors is usually $\frac{1}{\sigma}$ for variance parameters ("jeffreys" prior). But you will be able to see that if you were to use jeffreys prior for both parameters, you would have an improper posterior. But note that the main justification for using jeffreys prior is that it is a scale parameter. However you can show for your model, that neither parameter sets the scale of the problem.

If we consider the marginal model, with $\theta_{i}$ integrated out. It is a well-known result that if you integrate a normal with another normal, you get a normal. So we can skip the integration, and just work out the expectation and variance. We then get:

$$E(y_{i}|\mu\sigma\sigma_{\theta})=E\left[E(y_{i}|\mu\sigma\sigma_{\theta}\theta_{i})\right]=E\left[\theta_{i}|\mu\sigma\sigma_{\theta}\right]=\mu$$ $$V(y_{i}|\mu\sigma\sigma_{\theta})=E\left[V(y_{i}|\mu\sigma\sigma_{\theta}\theta_{i})\right]+V\left[E(y_{i}|\mu\sigma\sigma_{\theta}\theta_{i})\right]=\sigma^{2}+\sigma_{\theta}^{2}$$

And hence we have the marginal model:

$$(y_{i}|\mu\sigma\sigma_{\theta})\sim N(\mu,\sigma^{2}+\sigma_{\theta}^{2})$$

And this does show an identifiability problem with this model - so the data cannot distinguish between the two variances, it can only give information about their sum. You may have been able to see this intuitively. For example, we can always take $\theta_{i}=y_{i}$ for all $i$ and hence this will set $\sigma=0$. Alternatively we can set $\theta_{i}=\mu$ for all $i$ and this will set $\sigma_{\theta}=0$. Both of these scenarios will be indistinguishable by the data - in the sense that if I was to generate two data sets, one from the first case, and one from the second (but ensured that $\sigma^{2}+\sigma_{\theta}^{2}$ was the same in both cases), you would not be able to tell which data set came from which case. This suggests that it is fundamentally the sum that sets the scale and so we should apply jeffreys prior to the parameter $\tau^{2}=\sigma^{2}+\sigma_{\theta}^{2}$. Now suppose that $\tau^{2}$ was known, I would have thought a non-informative choice of prior for $\sigma^{2}$ would be uniform between $0$ and $\tau^{2}$ (for a more informative choice I would use a re-scaled beta distribution over this range). So we have the prior:

$$p(\tau^{2},\sigma^{2})\propto\frac{1}{\tau^{2}}\frac{I(0<\sigma^{2}<\tau^{2})}{\tau^{2}}$$

If we make the change of variables to $\sigma^{2},\tau^{2}\to\sigma,\sigma_{\theta}$ so that. We then get:

$$p(\sigma_{\theta},\sigma)\propto\frac{1}{(\sigma^{2}+\sigma_{\theta}^{2})^{2}}|\frac{\partial\sigma^{2}}{\partial\sigma}\frac{\partial\tau^{2}}{\partial\sigma_{\theta}}-\frac{\partial\sigma^{2}}{\partial\sigma_{\theta}}\frac{\partial\tau^{2}}{\partial\sigma}| =\frac{2\sigma\sigma_{\theta}}{(\sigma^{2}+\sigma_{\theta}^{2})^{2}}$$

Note that the non-identifiability is preserved in this prior because it is symmetric in its arguments. Another not so obvious symmetry is that if you were to integrate out either one of the variance parameters you would be left with the jeffreys prior for the other one:

$$\int_{0}^{\infty}\frac{2\sigma\sigma_{\theta}}{(\sigma^{2}+\sigma_{\theta}^{2})^{2}}d\sigma=\frac{1}{\sigma_{\theta}}$$

Hence, all you are required to input is the prior range for one of the parameters, as this will stop you from getting into trouble with improper priors. Call this $0<L_{\sigma}<\sigma<U_{\sigma}<\infty$. It is then easy to sample from the joint density using the inverse CDF method, for we have:

$$F_{\sigma}(x)=\frac{\log\left(\frac{x}{L_{\sigma}}\right)}{\log\left(\frac{U_{\sigma}}{L_{\sigma}}\right)}\implies F^{-1}_{\sigma}(p)=\frac{U_{\sigma}^{p}}{L_{\sigma}^{p-1}}$$ $$F_{\sigma_{\theta}|\sigma}(y|x)=1-\frac{x^{2}}{y^{2}+x^{2}}\implies F^{-1}_{\sigma_{\theta}|\sigma}(p|x)=x\sqrt{\frac{p}{1-p}}$$

So you sample two independent uniform random variables $q_{1b},q_{2b}$, and then your random value of $\sigma^{(b)}=U_{\sigma}^{q_{1b}}L_{\sigma}^{1-q_{1b}}$ and your random value of $\sigma^{(b)}_{\theta}=U_{\sigma}^{q_{1b}}L_{\sigma}^{1-q_{1b}}\sqrt{\frac{q_{2b}}{1-q_{2b}}}$. Combine this with the usual flat prior for $-\infty<L_{\mu}<\mu<U_{\mu}<\infty$ generated by a third random uniform variable $\mu^{(b)}=L_{\mu}+q_{3b}(U_{\mu}-L_{\mu})$ and you have all the ingredients to do monte carlo posterior simulation - note that this is much better than "Gibbs sampling" because each simulation is independent, so no need to wait for convergence (and also less need for a large number of simulations) - and you are dealing with proper priors - so divergence is impossible (however some moments may or may not exist, but all quantiles exist).

• Interesting prior! I was about to object to your claim that there is only information about the sum of the variance components but noticed that he has $i$ indexing both $y$ and $\theta$; in my head I'd imputed the hierarchical model $y_{ij}|-\sim N(\theta_j, \sigma^2)$, $\theta_j\sim N(\mu, \sigma^2_\theta)$... – JMS May 22 '11 at 4:17
• But assuming he meant what he wrote, what's the point? Even with restrictions on $\sigma$ he can't estimate both $\sigma$ and $\sigma_\theta$ separately... – JMS May 22 '11 at 4:30
• The prior range restriction effectively breaks the symmetry, and hence identifiability problem. They still will remain negatively correlated in the posterior. However it does this in a conservative and defensible way - you can usually provide a reasonable case for the choice of upper and lower limits without the need for any data. – probabilityislogic May 22 '11 at 13:02
• How does the range restriction solve the identifiability problem? I could still lop a little off the one and add it to the other, respecting interval constraints - hence, unidentified and not (separately) estimable. While there are of course situations where having portions of the model unidentified is sensible, this isn't one of them. – JMS May 22 '11 at 15:22