We now present a full treatment of a simple hierarchical model based on the normal distribution, in which observed data are normally distributed with a different mean for each ‘group’ or ‘experiment,’ with known observation variance, and a normal population distribution for the group means. This model is sometimes termed the one-way normal random-effects model with known data variance and is widely applicable, being an important special case of the hierarchical normal linear model.
The data structure
Consider $J$ independent experiments, with experiment $j$ estimating the parameter $\theta_j$ from $n_j$ independent normally distributed data points $y_{ij}$, each with known error variance $\sigma^2$; that is,
$y_{ij} \mid \theta_j \sim N(\theta_j, \sigma^2), \quad i = 1, \ldots, n_j; \; j = 1, \ldots, J$
Using standard notation from the analysis of variance, we label the sample mean of each group $j$ as
$\bar{y}_{\cdot j} = \frac{1}{n_j} \sum_{i=1}^{n_j} y_{ij},$
with sampling variance $\sigma_j^2 = \frac{\sigma^2}{n_j}.$
We can then write the likelihood for each $\theta_j$ using the sufficient statistic $\bar{y}_{\cdot j}$:
$\bar{y}_{\cdot j} \mid \theta_j \sim N(\theta_j, \sigma_j^2)$
a notation that will prove useful later because of the flexibility in allowing a separate variance $\sigma_j^2$ for each group $j$.
This assumption of known variances is often adequate in practice even if not strictly true.
Constructing a prior distribution from pragmatic considerations
Rather than immediately specifying a prior for the parameter vector $\theta = (\theta_1, \ldots, \theta_J)$, consider what kinds of posterior estimates might be reasonable for $\theta$, given the data $y_{ij}$.
A simple natural approach is to estimate $\theta_j$ by its sample mean $\bar{y}_{\cdot j}$.
But suppose there are $J = 20$ experiments with only $n_j = 2$ observations per group, and the groups represent 20 nearly identical experimental conditions.
The two observations per group are insufficient for accurate estimation, so we might prefer a pooled estimate:
$\bar{y}_{\cdot \cdot} = \frac{\sum_{j=1}^J \frac{\bar{y}_{\cdot j}}{\sigma_j^2}} {\sum_{j=1}^J \frac{1}{\sigma_j^2}}.$
To decide between using individual or pooled estimates, a traditional approach is an ANOVA $F$-test:
if between-group variance is significantly larger than within-group variance, estimate groups separately; otherwise, pool them.
From classical to Bayesian interpretation
In classical random-effects ANOVA, one computes mean squares and compares them.
Let $\tau^2$ denote the variance of $\theta_1, \ldots, \theta_J$.
For a balanced design ($n_j = n$, so $\sigma_j^2 = \sigma^2 / n$), the ANOVA table is:
| Source | df | SS | MS | $E(\text{MS} \mid \sigma^2, \tau^2)$ |
|---|---|---|---|---|
| Between groups | $J – 1$ | $\sum_j (\bar{y}_{\cdot j} – \bar{y}_{\cdot \cdot})^2$ | $\text{SS}/(J-1)$ | $\tau^2 + \sigma^2$ |
| Within groups | $J(n – 1)$ | $\sum_{i,j} (y_{ij} – \bar{y}_{\cdot j})^2$ | $\text{SS}/(J(n-1))$ | $\sigma^2$ |
| Total | $Jn – 1$ | $\sum_{i,j} (y_{ij} – \bar{y}_{\cdot \cdot})^2$ |
If the ratio of between- to within-group mean squares is large, the groups differ ($\tau^2 > 0$); otherwise, pooling ($\tau^2 = 0$) is reasonable.
Bayesian interpretation of pooling
The Bayesian model provides a continuum between “no pooling” and “complete pooling.”
The posterior mean of $\theta_j$ can be expressed as a weighted combination:
$\hat{\theta}_j = \lambda_j \bar{y}_{\cdot j} + (1 – \lambda_j)\bar{y}_{\cdot \cdot}, \quad 0 \le \lambda_j \le 1$
Different priors lead to the following posterior estimates:
- No pooling:
- $\hat{\theta}_j = \bar{y}_{\cdot j}$
- if $\theta_j$ have independent uniform priors on $(-\infty, \infty)$.
- Complete pooling:
- $\hat{\theta}_j = \bar{y}_{\cdot \cdot}$
- if all $\theta_j$ are constrained to be equal with a uniform prior on the common value.
- Partial pooling:
- $\hat{\theta}_j$ as above if $\theta_j \stackrel{iid}{\sim} N(\mu, \tau^2)$.
No pooling corresponds to $\tau^2 \to \infty$ ($\lambda_j = 1$),
complete pooling to $\tau^2 = 0$ ($\lambda_j = 0$).
The hierarchical model
For conjugacy, assume $\theta_j \mid \mu, \tau \sim N(\mu, \tau^2), \quad j = 1, \ldots, J.$
We assign a noninformative uniform hyperprior to $\mu$ given $\tau$:
$p(\mu, \tau) = p(\mu \mid \tau)p(\tau) \propto p(\tau)$
The joint posterior distribution
Combining the sampling model and prior gives the joint posterior:
$p(\theta, \mu, \tau \mid y) \propto p(\mu, \tau) \prod_{j=1}^J N(\theta_j \mid \mu, \tau^2) \prod_{j=1}^J N(\bar{y}_{\cdot j} \mid \theta_j, \sigma_j^2)$
Conditional posterior of group means given hyperparameters
Because the $\theta_j$ are independent in both prior and likelihood,
their conditional posterior distributions factorize:
$\theta_j \mid \mu, \tau, y \sim N(\hat{\theta}_j, V_j)$
where
$\hat{\theta}_j = \frac{\bar{y}_{\cdot j}/\sigma_j^2 + \mu / \tau^2} {1/\sigma_j^2 + 1/\tau^2}, \qquad V_j = \frac{1}{1/\sigma_j^2 + 1/\tau^2}$
The posterior mean is a precision-weighted average of the sample mean and population mean.
Marginal posterior distribution of hyperparameters
Integrating over $\theta$, the group means have marginal distribution $\bar{y}_{\cdot j} \mid \mu, \tau \sim N(\mu, \sigma_j^2 + \tau^2)$,
so
$p(\mu, \tau \mid y) \propto p(\mu, \tau) \prod_{j=1}^J N(\bar{y}_{\cdot j} \mid \mu, \sigma_j^2 + \tau^2).$ (5.18)
We can factor this as $p(\mu, \tau \mid y) = p(\mu \mid \tau, y) p(\tau \mid y).$
Given $\tau$, the posterior for $\mu$ is normal:
$\mu \mid \tau, y \sim N(\hat{\mu}, V_\mu)$
where
$\hat{\mu} = \frac{\sum_{j=1}^J \bar{y}_{\cdot j}/(\sigma_j^2 + \tau^2)} {\sum_{j=1}^J 1/(\sigma_j^2 + \tau^2)}, \qquad V_\mu = \left(\sum_{j=1}^J \frac{1}{\sigma_j^2 + \tau^2}\right)^{-1}$
Substituting into (5.18), we get:
$p(\tau \mid y) \propto p(\tau) V_\mu^{1/2} \prod_{j=1}^J (\sigma_j^2 + \tau^2)^{-1/2} \exp\!\left[-\frac{(\bar{y}_{\cdot j} – \hat{\mu})^2}{2(\sigma_j^2 + \tau^2)}\right]$
Prior distribution for τ\tauτ
A convenient noninformative choice is $p(\tau) \propto 1$, which yields a proper posterior in this model.
By contrast, $p(\log \tau) \propto 1$ produces an improper posterior.
In practice, weakly informative priors (e.g., scaled inverse-$\chi^2$) are often preferred.
Computation
The posterior factorizes as $p(\theta, \mu, \tau \mid y) = p(\tau \mid y)p(\mu \mid \tau, y)p(\theta \mid \mu, \tau, y).$
Simulation steps:
- Draw $\tau$ from $p(\tau \mid y)$ using a numerical grid.
- Given $\tau$, draw $\mu \sim N(\hat{\mu}, V_\mu)$.
- Given $(\mu, \tau)$, draw each $\theta_j \sim N(\hat{\theta}_j, V_j)$.
Posterior predictive distributions
- For existing groups:
- Draw $(\theta, \mu, \tau)$ from the posterior and simulate
- $\tilde{y}_{ij} \mid \theta_j \sim N(\theta_j, \sigma^2)$.
- Draw $(\theta, \mu, \tau)$ from the posterior and simulate
- For new groups:
- Draw $(\mu, \tau)$ from the posterior, then
- $\tilde{\theta}_j \sim N(\mu, \tau^2)$,
- followed by $\tilde{y}_{ij} \mid \tilde{\theta}_j \sim N(\tilde{\theta}_j, \sigma^2)$
- Draw $(\mu, \tau)$ from the posterior, then
Difficulty with plug-in (non-Bayesian) estimation
Classical random-effects ANOVA estimates:
$\hat{\mu} = \bar{y}_{\cdot \cdot}, \qquad \hat{\tau}^2 = \frac{MS_B – MS_W}{n}$
Inference based on $p(\theta \mid \hat{\mu}, \hat{\tau}, y)$ ignores uncertainty in $\mu$ and $\tau$,
and $\hat{\tau}^2$ can even be negative.
Setting $\hat{\tau}^2 = 0$ when $MS_W > MS_B$ implies all $\theta_j$ are identical,
which is too strong.
A full Bayesian approach integrates over the uncertainty in $\tau$, avoiding such issues.