In the book Bayesian Data Analysis by Gelman et al. (3rd edition, 2014), a hierarchical model (or one-way random-effects ANOVA) is presented in section 5.4 as follows,

\begin{equation}\label{eq:lme1} y_{ij} = b_0 + \lambda_i + \varepsilon_{ij}, \end{equation}

where the data $y_{ij}$ come from the $i$th measuring entity (e.g., student performance in a school district) collected under the $j$th condition (e.g., a school within the district), $b_0$ is the population mean, $\lambda_i$ is the deviation of the $i$th measuring entity from the population mean, and $\varepsilon_{ij}$ is the measuring error ($i=1, 2, \ldots, k;\ j=1, 2, \ldots, n$).

A posterior inference is derived in the book for the effect of each measuring entity $\theta_i=b_0 + \lambda_i$ based on a Gaussian assumption with a known variance $\sigma^2$ for the residuals $\varepsilon_{ij}$ and a prior distribution $G(0, \tau^2)$ for $\lambda_i$. Specifically, the mean and variance for $\theta_i$ are estimated as below:

\begin{align} {\rm mean}(\theta_i) &= \frac{\frac{n}{\sigma^2}\bar{y}_{i\cdot}+\frac{1}{\tau^2}b_0}{\frac{n}{\sigma^2}+\frac{1}{\tau^2}} \\[7pt] {\rm Var}(\theta_i) &= \frac{1}{\frac{n}{\sigma^2}+\frac{1}{\tau^2}} \end{align}

where $\bar{y}_{i\cdot}=\frac{1}{n}\sum_{j=1}^n y_{ij}$.

Even though the variance for the $\lambda_i$ is assumed to be known, I could solve the model as a mixed-effects model through, for example, function lmer() in the R package lme4, and use the estimated variances $\tau^2$ and $\sigma^2$ to obtain the posterior distribution using the formulation above. Is this a reasonable and solid approach?

I know that I could directly obtain the posterior distribution through R packages such as brms and rstanarm. However, the computational cost is too heavy in my case, and that's why I'm trying to see if the above closed form is a a reasonable approach to directly obtaining the posterior distribution by plugging the variance estimate $\hat{\sigma}^2$ from lmer, rather than going through the typical Bayesian route..

  • 1
    $\begingroup$ You might be better to split this up and ask question 1 first and then when you have a satisfactory answer ask 2 with a link back to 1, and so on. $\endgroup$
    – mdewey
    Jul 12, 2017 at 10:30
  • $\begingroup$ Having been edited to constrain the question, I think this is no longer too broad. I'm voting to leave open. $\endgroup$ Jul 12, 2017 at 15:12
  • $\begingroup$ should it be $j$ for measuring identity, $i$ for measuring condition (with $\lambda_i$ for the condition instead of identity)? Currently it is $i$-th identity and $k$-th condition. With $\lambda$ for identity. Another point... the 'estimates' means maximum posteriori probability? $\endgroup$ Jul 18, 2017 at 20:34
  • $\begingroup$ @MartijnWeterings Thanks. $k$ should be $j$, and think of $i$ and $j$ as indices for school districts and schools. What "estimates" are you referring to: mean and variance of $\theta_i$? $\endgroup$
    – bluepole
    Jul 18, 2017 at 21:37

1 Answer 1


The general concept here is updating a conjugate prior.

I'm trying to see if the above closed form is a feasible solution

I'm not sure what you mean by feasible. For the Gaussian distribution $N\left(\mu, \sigma^2\right)$, any specification of mean $\mu$ and variance $\sigma^2$ - with mean parameter unrestricted and variance parameter strictly positive $\sigma^2 > 0$ - defines a unique and valid Gaussian distribution.

I will assume you are asking if the proposed closed form solution is correct.

The closed form solution you provide is the correct posterior distribution given a Gaussian prior distribution $b_0+\lambda_i = \theta_i \sim N\left(b_0, \tau_0^2\right)$ and observed data with known observation variance $\sigma^2$.

Precision $K$ is the inverse of variance $\sigma^2$: $K = \sigma^{-2}$.

Considering the mean of the posterior distribution for $\theta_i$, the numerator in the formula $mean(\theta_i)$ you provide is the precision weighted combination of the prior distribution's mean $b_0$ and the observed data's mean $\bar{y_i}$. The denominator of the posterior mean formula simply scales the sum of these weights to one.

Considering the variance of the posterior distribution, the operation reflected in the formula $Var(\theta_0)$ you provide is the summing of the precision of the prior distribution and the precision of the observations; then inverting to obtain variance.

In summary, the formulas you provide for posterior mean and posterior variance are correct given the prior $\theta_i \sim N(b_0,\tau^2)$, and given $n$ observations with known observation variance $\sigma^2$.

If your concern has to do with specifying the prior distribution $N(b_0,\tau^2)$ parameters, note the influence of the prior mean $b_0$ in your formula for the posterior mean $mean(\theta_i)$ diminishes as the number of observations increases, and diminishes as the prior variance increases. Therefore, don't worry too much about the prior mean, just increase the prior variance as seems appropriate; and, increase the number of observations!

  • $\begingroup$ Thanks for the answer! Sorry that my original question was unclear. I was not questioning whether the close form of the posterior distribution is correct given the known variance $\sigma^2$ and number of observations, $n$. Instead, I want to know whether using estimated $\hat{\sigma}^2$ from solving LME through function lmer in the R package lme4 is a reasonable approach to directly obtaining the posterior distribution, rather than going through the typical Bayesian route. Could you comment about this? $\endgroup$
    – bluepole
    Jul 23, 2017 at 13:41
  • $\begingroup$ Ok, life just got more complicated :)You should look at Peter Hoff's notes, particularly Ch. 7 & 8. $\endgroup$
    – krkeane
    Jul 23, 2017 at 14:12
  • $\begingroup$ Apparently I'm not allowed to edit after 5 minutes. The concept is now updating a semi-conjugate prior (unknown mean and unknown variance). The correct chapters in Hoff's notes are Ch. 6 & 7; or better yet, look at his book, Hoff, Peter D. A first course in Bayesian statistical methods. Springer Science & Business Media, 2009. Hoff provides simple R code in both his notes and book. So, if you want to perform a Bayesian analysis, I would study Hoff's Ch. 6 & 7 in his notes, or similar chapters in his book (don't have the book handy today); and not use lmer $\endgroup$
    – krkeane
    Jul 23, 2017 at 14:21
  • $\begingroup$ Why is using lmer to directly estiamte $\hat{\sigma}^2$ problematic for the posterior distribution? Andrew Gelmen seems to hint that this would be fine in this paper: stat.columbia.edu/~gelman/research/published/multiple2f.pdf $\endgroup$
    – bluepole
    Jul 24, 2017 at 15:46

Your Answer

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

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