# Bayesian Hierarchical Model - Exact Conjugate Solution?

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.

• By the mixture property of the normal distribution, $Y_i | \theta_i$ is given by $N(\mu, \frac{1}{\tau} + \frac{1}{\tau_i})$. Is $\tau_i$ observed in each study, and assumed drawn from the same gamma distribution? Commented Jun 1, 2014 at 18:40
• Yes $\tau_i$ is observed in each study. It is not assumed to be drawn from the same gamma Commented Jun 1, 2014 at 19:18
• Do you have any other prior beliefs for $\tau_i$, or any assumptions about how it should be modeled? I can't quite suss out why there are two sources of variance. Commented Jun 1, 2014 at 20:18
• $\tau$ is the between studies precision which is unknown. $\tau_i$ is the within-study precision which is known from the data. Commented Jun 2, 2014 at 0:25
• That clarifying information should probably go in your question Commented Jun 2, 2014 at 1:33

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.

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