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?
 A: 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.
A: 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).
