Bayesian Hierarchical Model - Exact Conjugate Solution?

Question

I was hoping to get some help. In understand how to compute an exact numerical solution (http://www.cs.berkeley.edu/~jordan/courses/260-spring10/lectures/lecture5.pdf) for the following Bayesian model: $$ \tau \sim Ga(\alpha, \beta)$$ $$\mu \sim N(m ,p)$$ $$Y_i \sim N(\mu, \tau)$$ where the data is $Y_i$.

I am performing a meta-analysis where my data are $Y_i$ from each study, and $\tau_i$ from each study. The problem is described by a hierarchical model: $$ \tau \sim Ga(\alpha, \beta)$$ $$\mu \sim N(m ,p)$$ $$\theta_i \sim N(\mu, \tau)$$ $$Y_i \sim N(\theta_i, \tau_i)$$ Is there any way to solve this model analytically? Any help/suggestions would be greatly appreciated.

Answer 1

I'm assuming that by "exact numerical solution" and "analytically" you mean that you want $$p(\mu,\theta_i,\tau|Y_1,\ldots,Y_n, \tau_1,\ldots,\tau_n)$$ to be a known distribution. If this is what you want, then NO there is no way to solve for this posterior analytically.

Let me use $\sigma^2$ instead of your $\tau$, so that we have $$\theta_i \stackrel{iid}{\sim} N(\mu,\sigma^2)$$. Some suggestions that would get you closer to an analytic solution are

• Let $Y_i\stackrel{iid}{\sim} N(\theta_i,\sigma^2\tau_i)$.
• Let $\mu \sim N(m,\sigma^2p)$
• Let $\sigma^2 \sim IG(\alpha,\beta)$.

Now you should be able to integrate out $\theta_1,\ldots,\theta_n,\mu$ to obtain $Y\sim N(m, \sigma^2 S)$ where $Y=(Y_1,\ldots,Y_n)$and $S$ is a known matrix. Then, you should be able to get $\mu|Y,\tau_1,\ldots,\tau_n,\sigma^2$ followed by $\theta_1,\ldots,\theta_n|Y,\tau_1,\ldots,\tau_n,\mu,\sigma^2$. Although this wouldn't get you to an analytic solution, you can at least use Monte Carlo (rather than Markov chain Monte Carlo) to obtain draws from the full posterior.

Answer 2

Unfortunately the algebraic structure of the normal likelihood does not allow us to separate the $\tau_i$ and $\tau$ in a way that allows for exact conjugate relationships to be used here. The exponential family conjugate relationships are a direct consequence of the sum/product properties of exponentials.. to see the problem look at the log likelihood of the data: $$ \text{LL}(\text{data}) = \text{constant} + \frac{1}{2}\sum_i \log(\tau_i) + \frac{1}{2}\sum_i \tau_i (Y_i - \theta_i)^2. $$ There is no way to combine terms involving $\theta_i$ with the prior for $\theta_i$, $$ \log(p(\theta_i)) = \text{constant} + \frac{1}{2} \log(\tau) + \frac{1}{2} \tau (\theta_i - \mu)^2. $$

To combine these two (and the other distributions involved naturally) and get a posterior distribution as a function of the parameters, we would normally combine the squared sum at the end of each of these through completing the square (see section 3 here). This is not possible here because each of the $\tau_i$ infront of each $(Y_i - \theta_i)^2$ are not the same. By not being able to complete the square we are not able to put the parameters in a normal distribution form, so there is by definition no conjugacy. This is without even investigating the additional levels of depth necessitated by the prior (it's really just one Normal-Gamma prior) on $\tau$ and $\mu$.

Further, this non-conjugacy result would hold if you use the convolution of normals mentioned by SeanEaster due to the way the variances are added, making them impossible to wrangle into a Normal-Gamma posterior.

---

An alternative approach to working with a similar model specification is given by jaradniemi