How to interpret the population parameters of a Bayesian Hierarchical model?

Question

This is almost certainly a fatal misunderstanding of mine / knowledge gap but I am confused as to how to interpret the population parameters of a Bayesian Hierarchical model.

This is incredibly artificial but to demonstrate the point let's say we wanted to fit the following model + priors:

$$ \begin{align} X_{ij} &\sim N(\mu_j, \sigma) \\ \mu_j &\sim \text{LogNormal}(\mu, \tau) \\ \mu &\sim \text{N}(2, 1) \\ \sigma &\sim \text{LogNormal}(0.5, 0.5) \\ \tau &\sim \text{LogNormal}(0.5, 0.5) \\ \end{align} $$

Now lets say our "question of interest" that we want to answer is "what is the population-mean of the $\mu_j$'s ?".

Naively I would assume that because they follow a Log-Normal distribution that the mean of the distribution is then $exp(\mu + \tau^2/2)$ (wiki) where we just sub in the posterior samples of $\mu$ and $\tau$ to get an uncertainty distribution.

But this doesn't feel right to me...

For example the posterior of $\tau$ is not necesarily a Log-Normal distribution anymore (that was just our prior), so why would the posterior population distribution of $\mu_j$ still be a Log-Normal distribution ? But if its not a Log-Normal distribution what do the values of $\mu$ and $\tau$ even mean anymore then; do they even have an interpretation if the underlying distributions are different in the posterior ?

Any guidance as to where my mental model has gone wrong would be greatly appreciated.

Answer 1

I think I have found an explanation for why $\psi(\mu, \tau) := \exp(\mu + \tau^2/2)$ can be useful to answer the question "What is the population-mean of the $\mu_j$'s?":

Imagine you have an overall population that you are sampling your groups $j$ from, and then you are sampling units $i$ from these groups in turn. Say $j \in \cal{J}$, where $\cal{J}$ is the whole population of $j$'s, and you have drawn $J$ of them in your sample, so you observe data $X_{ij}, j = 1, \dotsc, J$, and for each of them units $i = 1, \dotsc, n_j$, say. Let's call this overall data $X$.

Now you calculate your posterior distributions via MCMC etc., in particular obtaining $M$ samples $\mu_m, \tau_m$, $m=1, \dotsc, M$, from the posterior given $X$.

And finally you are interested in statements about a new group $j = J+1$ from the same population $\cal{J}$, for which you don't have any data $X_{ij}$. It makes sense to assume that this group is exchangeable with the other groups that you have observed $X_{ij}$ for.

The important fact now is that the full conditional distribution of the new group's mean $\mu_{J+1}$ is again $\text{LogNormal}(\mu, \tau)$, i.e. for the density we have:

$f(\mu_{J+1} | X, \mu_1, \dotsc, \mu_J, \mu, \sigma, \tau) = f(\mu_{J+1} | \mu, \tau)$

and this is the density of the $\text{LogNormal}(\mu, \tau)$ distribution. Therefore, in particular we have

$\mathbb{E}(\mu_{J+1} | X, \mu_1, \dotsc, \mu_J, \mu, \sigma, \tau) = \mathbb{E}(\mu_{J+1} | \mu, \tau) = \psi(\mu, \tau)$

as defined above. So, given we know everything else, including the data $X$ from the other $J$ groups, $\psi$ is the mean for $\mu_{J+1}$.

Now that is not yet the marginal posterior expectation yet, because we still condition on $\mu$ and $\tau$. However, based on the law of iterated expectations we have

$\mathbb{E}(\mu_{J+1} | X) = \mathbb{E}[\mathbb{E}(\mu_{J+1} | \mu, \tau) | X] = \mathbb{E}[\psi(\mu, \tau) | X]$

where the outer expectation is with regards to the posterior of $\mu$ and $\tau$, and we can approximate this marginal mean easily by the Monte Carlo estimate

$\hat\psi = \frac{1}{M}\sum_{m = 1}^{M} \psi(\mu_m, \tau_m)$

based on samples $\mu_m, \tau_m$ from the posterior given $X$.

Two notes at the end:

1. It is important that even in the prior model, the marginal prior distribution for the $\mu_j$'s is actually not log normal. Only the full conditional prior distribution, i.e. conditioning on all other parameters, is log normal. This then carries over to the posterior predictive distribution for $\mu_{J+1}$, as shown above.
2. The overall concept here is applied a lot in the meta-analytic prior context (https://onlinelibrary.wiley.com/doi/10.1111/biom.12242)

Answer 2

the posterior of $\tau$ is not necesarily a Log-Normal distribution anymore (that was just our prior), so why would the posterior population distribution of $\mu_j$ still be a Log-Normal distribution ?

Your intuition is correct here. The posteriors might not be Log-Normal distributions, even if their priors are. In fact, the posteriors might not resemble any closed-form distribution.

But if its not a Log-Normal distribution what do the values of $\mu$ and $\tau$ even mean anymore then ?

$\mu$ and $\tau$ mean the same as they meant before. Except then you didn't have data, so you assumed that they followed Log-Normal distributions (based on prior experience, or maybe domain knowledge); now you have data, so you incorporate that into your estimate. You didn't discard your previous assumption (that they're log-normally distributed), but you don't discard the data that you have now either: you "mix" the two together (via Bayes' Theorem).

Answer 3

EDIT - I want to add that I no longer think my answer below is correct (i'm still unsure as to what I think the answer is) but am leaving my answer up as it already has comments.

(posting an answer to my own question) I'm not sure if this is correct but after thinking about this for ages (and based on the answer by @LmnICE) the following seems to make some intuiative sense to me, please comment/correct me if I have this wrong

---

I think the confusion is over what the posterior actually represents. Being explicit our posterior as noted in the OP can be expressed as:

$$ P(\mu_j, \mu, \tau, \sigma \ \mid X) \propto P(X \mid \mu_j, \sigma)P(\mu_j \mid \mu, \tau)P(\tau)P(\sigma)P(\mu) $$

That is to say that a sample of $\mu_j, \mu, \tau, \sigma$ from our posterior distribution is equivalent to a sample from the model which assumes these relationships between the parameters.

Or put another way, the posterior hasn't changed any of these assumed relationships; in fact the posterior is dependent on them. E.g. a sample of $\mu$ from the posterior is equivalent to a sample of $\mu$ where we are assuming it is the log-centrality parameter of a log-normal distribution for $\mu_j$.

I think the confusion I had is that we often refer to the posterior as being an updated / changed distribution but this is merely referencing the fact that $P(\theta \mid X) \ne P(\theta)$. However under this model the data doesn't change the relationship between parameters e.g. $P(\mu_j \mid \mu, \tau)$ is still $P(\mu_j \mid \mu, \tau)$, all that is changed is what values we think those parameters are.

So to answer "question of interest" I would argue that as our model is assuming a log-normal releationship that the population-mean of $\mu_j$ is indeed $exp(\mu + \tau^2/2)$. Likewise then the interpretation of the values of $\mu$ and $\tau$ that come out of the posterior are the values that most likely fit this assumed model given the available data.