6
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ So $\mu$ is the global mean of $X$, $\mu_j$ is the mean of each $X_j$ and $X_j$ is a specific observation. In your example, you only have one observation “pulling” its own subgroup mean, so the heirchical model isn’t necessary/beneifical. $\endgroup$
    – jbuddy_13
    Dec 1, 2023 at 10:18
  • $\begingroup$ However if you replaced $\mu_j$ with $\mu_i$ (mean of group i) and $X_j$ with $X_{i,j}$ (the jth observation of group i) then this design suddenly adds value. Now you can estimate the global mean, and each subgroup mean. $\endgroup$
    – jbuddy_13
    Dec 1, 2023 at 10:21

3 Answers 3

3
$\begingroup$

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)
$\endgroup$
4
$\begingroup$

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).

$\endgroup$
3
  • $\begingroup$ Sorry just to clear up ambiguity I am assuming your "that they're log-normally distributed" is refering to the mu_j ? $\endgroup$
    – gowerc
    Dec 1, 2023 at 10:02
  • 2
    $\begingroup$ @gowerc, this commentary is valid for any variable with a log normal prior. Except when dealing with a class of models called “conjugate priors”, Bayesian models do not have closed form solutions; this means you cannot analytically derive mean, variance, etc for a given random variable. Samplers (Hamiltonian Monte Carlo / NUTS is very common) will sample the posterior, returning “chains” (the nth element is an array indicating the sampled value of each rv at the nth step in the chain.) These sampled values can and often will diverge from the priors specified, depending on the data observed. $\endgroup$
    – jbuddy_13
    Dec 1, 2023 at 10:32
  • $\begingroup$ Yes. Generally the priors are assumptions about each variable in the model. So in your case, $\mu_j$, $\tau$ and $\sigma$ are assumed to be log-normal before seeing the data, and $X_{ij}$ and $\mu$ are assumed to be normally distributed. $\endgroup$
    – LmnICE
    Dec 1, 2023 at 10:33
1
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ Your answer is almost entirely correct, but like I said in my answer, the posterior of the $\mu_j$ might not be lognormally distributed; therefore the closed-form formula for the mean of the lognormal is not appropriate. If you estimate the posterior by way of simulation, why not just calculate the mean of the simulation samples? This is simpler and distro-agnostic. $\endgroup$
    – LmnICE
    Dec 1, 2023 at 12:20
  • 1
    $\begingroup$ I think if you calculate the mean via taking the mean of the posterior samples you only get uncertainty estimates for the sample statistic and not the population statistic. I'm not sure I fully understand your comment though; the posterior of mu_j is defined by the assumption that they are log-normally distributed with tau/sigma. The mean is then being analytically calculated under that assumption. Surely the posterior is then the values that best meet that assumption. Is this not equivalent to the assumption of the likelihood being a distribution which may actually be wrong? $\endgroup$
    – gowerc
    Dec 1, 2023 at 12:44
  • $\begingroup$ The priors of each $\mu_j$ are defined as being lognormally distributed, not the posteriors. Assuming you have data for the $X_{ij}$, its Normal priors will get "mixed" (via Bayes' Theorem) with the available data, yielding a new distribution (possibly without a closed-form formula). The other variables in the model (including the $mu_j$) will then get updated to reflect the new distribution for each $X_{ij}$. (This isn't how sampling from the posterior via Monte Carlo simulation actually works; I'm just trying to give you an intuition). $\endgroup$
    – LmnICE
    Dec 3, 2023 at 1:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.