Hierarchical and Empirical Bayes

Intuition Beginner

Imagine you must estimate the same kind of quantity for many groups at once — the success rate of a treatment in each of eight hospitals, say. You could pool everything into one number and pretend the hospitals are identical, or you could estimate each hospital alone and ignore the others. Both extremes waste information. The hospitals are not identical, but they are similar, and a sensible estimate for a hospital with little data should lean on what the other hospitals reveal about the typical rate.

A hierarchical model captures this middle ground. It says each group has its own true value, and those values are themselves a sample from a shared population of groups. You put a prior on each group's value, and a second prior — a hyperprior — on the population the groups are drawn from. The data then teach you two things at once: where each group sits, and what the whole population looks like.

The payoff is borrowing strength. A group with noisy or sparse data gets pulled toward the population average, because the model trusts the crowd when one group is uncertain. A group with plenty of clean data barely moves, because its own evidence speaks for itself. This selective pulling — strong for the uncertain, gentle for the confident — is called partial pooling, and it usually beats both extremes.

Visual Beginner

Picture three layers stacked like a family tree. At the top sits the population description (its center and spread). Below it hang the individual group values, each one a child of the population. At the bottom sit the data points, each a noisy reading of its group's value.

Layer	What lives here	Role
Top (hyperprior)	population center and spread	describes how groups vary
Middle (group values)	one true value per group	borrows from the population
Bottom (data)	noisy measurements	speaks for its own group

Read the arrows downward and you have the model: population shapes groups, groups shape data. Read them upward and you have the inference: data inform the groups, and pooled across groups they also reveal the population, which then pulls the uncertain groups inward.

Worked example Beginner

Two clinics measure the average recovery time for a procedure. The true clinic averages are unknown. Each clinic's measurement has noise of standard deviation $\sigma =4$ days. Past experience says clinic averages scatter around a population center $\mu =30$ days with spread $\tau =3$ days.

Clinic A reports a measured average of $24$ days; clinic B reports $36$ days.

Step 1. Find the pull factor. The shrinkage factor measures how much weight the data deserve against the population. It is $B={\tau }^2/({\tau }^2+{\sigma }^2)=9/(9+16)=9/25=0.36$.

Step 2. Shrink clinic A. The hierarchical estimate moves the data toward the population center: $30+0.36\times (24-30)=30+0.36\times (-6)=30-2.16=27.84$ days.

Step 3. Shrink clinic B. Same rule: $30+0.36\times (36-30)=30+2.16=32.16$ days.

Step 4. Compare. Clinic A's raw $24$ became $27.84$; clinic B's raw $36$ became $32.16$. Both were pulled about a third of the way to the population center $30$.

What this tells us. The noise ($\sigma =4$) is larger than the population spread ($\tau =3$), so the model distrusts any single clinic's reading and pulls it firmly toward $30$. If the clinics had been measured very precisely (small $\sigma$), the factor $B$ would be near $1$ and the estimates would barely move. The hierarchy automatically tunes the pull to how much each level knows.

Check your understanding Beginner

Formal definition Intermediate+

A hierarchical (multilevel) model specifies a joint distribution through a chain of conditional distributions. For $J$ groups with group parameters $\theta =({\theta }_1,\dots ,{\theta }_J)$, observed data $y=(y_1,\dots ,y_J)$, and a population-level hyperparameter $\varphi$, the model is

$$\varphi \sim p(\varphi ),{\theta }_j\mid \varphi \stackrel{\text{ind}}{\sim }p({\theta }_j\mid \varphi ),y_j\mid {\theta }_j\stackrel{\text{ind}}{\sim }p(y_j\mid {\theta }_j),$$

so that the joint density factors as

$$p(\varphi ,\theta ,y)=p(\varphi )\prod _{j=1}^{J}p({\theta }_j\mid \varphi )p(y_j\mid {\theta }_j).$$

The defining structural feature is conditional independence: given $\varphi$ the group parameters are independent, and given ${\theta }_j$ the datum $y_j$ is independent of everything else. The hyperprior $p(\varphi )$ is the prior on the population from which the ${\theta }_j$ are drawn; $p({\theta }_j\mid \varphi )$ is the population distribution (sometimes called the prior at the second level). The marginal (or compound) sampling density of a single group, with $\varphi$ held fixed, is

$$m_{\varphi }(y_j)=\int p(y_j\mid {\theta }_j)p({\theta }_j\mid \varphi )d{\theta }_j,$$

the data-generating density after integrating out the group parameter; it is the object 45.03.02 computes in closed form for conjugate pairs and the object empirical Bayes will estimate $\varphi$ from.

The Gaussian hierarchical normal-means model. Take the conjugate normal hierarchy with known sampling variance ${\sigma }^2$ and hyperparameter $\varphi =(\mu ,{\tau }^2)$:

$$y_j\mid {\theta }_j\sim N({\theta }_j,{\sigma }^2),{\theta }_j\mid \mu ,{\tau }^2\sim N(\mu ,{\tau }^2),j=1,\dots ,J.$$

With $\varphi$ fixed, the conjugate normal-normal update of 45.03.02 gives the group posterior

$${\theta }_j\mid y,\mu ,{\tau }^2\sim N(\mu +B(y_j-\mu ),(1-B){\sigma }^2),B=\frac{{\tau }^2}{{\sigma }^2+{\tau }^2},$$

where $B\in (0,1)$ is the shrinkage factor and $1-B$ measures how far each estimate is pulled from its own datum toward the population center $\mu$. The posterior mean ${\hat{\theta }}_j=\mu +B(y_j-\mu )=(1-B)\mu +B{y}_j$ is a convex combination of the population center and the group's datum: partial pooling interpolates between complete pooling ($B\to 0$, all ${\hat{\theta }}_j=\mu$) and no pooling ($B\to 1$, all ${\hat{\theta }}_j=y_j$).

Empirical Bayes. Parametric empirical Bayes (EB) replaces the integration over $\varphi$ by a point estimate $\hat{\varphi }$ obtained from the data through the marginal: the marginal MLE maximises the marginal likelihood

$$\hat{\varphi }=\arg \underset{\varphi }{\max }\prod _{j=1}^J{m}_{\varphi }(y_j),$$

and the EB estimator of each ${\theta }_j$ is the conjugate posterior mean evaluated at $\hat{\varphi }$. For the Gaussian hierarchy the marginal of $y_j$ given $\varphi$ is $y_j\sim N(\mu ,{\sigma }^2+{\tau }^2)$, so $\hat{\mu }=\bar{y}$ and $\widehat{{\sigma }^2+{\tau }^2}$ is read off the marginal sample variance, yielding $\hat{B}={\tau }^2/({\sigma }^2+{\tau }^2)$ estimated from the spread of the $y_j$ themselves.

Nonparametric (Robbins) empirical Bayes does not assume a parametric form for the population distribution $G(\theta )$ at all. It estimates the posterior functional directly from the marginal $m(y)=\int p(y\mid \theta )dG(\theta )$ — the $f$-modeling strategy — using identities that express the posterior mean in terms of $m$ alone. For the Poisson sampling model $y_j\mid {\lambda }_j\sim \text{Poisson}({\lambda }_j)$ with an arbitrary mixing distribution $G$ on the rates, Robbins' formula

$$\mathbb{E}[{\lambda }_j\mid y_j=y]=(y+1)\frac{m(y+1)}{m(y)}$$

recovers the posterior mean of the rate from ratios of the marginal frequencies, with no parametric assumption on $G$.

Counterexamples to common slips Intermediate+

• Empirical Bayes is not full hierarchical Bayes. EB plugs a point estimate $\hat{\varphi }$ into the posterior and treats it as known; hierarchical Bayes integrates $\varphi$ against its hyperprior. The EB posterior $p({\theta }_j\mid y,\hat{\varphi })$ therefore understates uncertainty, because it ignores the variability in $\hat{\varphi }$. The two agree to first order when $J$ is large (the hyperprior is overwhelmed and $\hat{\varphi }$ is sharp), and diverge when $J$ is small.

• A hyperprior is not optional decoration. Omitting the hyperprior — fixing $\varphi$ at a guessed value — collapses the hierarchy to an ordinary conjugate model with no borrowing of strength across groups in the choice of $\varphi$. The borrowing comes precisely from letting all groups jointly inform the population level.

• Improper hyperpriors can produce improper posteriors. The flat hyperprior $p({\tau }^2)\propto 1$ or $p(\tau )\propto 1$ may seem harmless, but for the Gaussian hierarchy the joint posterior can fail to be integrable, especially the prior $p({\tau }^2)\propto 1/{\tau }^2$, which puts infinite mass near $\tau =0$ and yields an improper posterior. Propriety must be checked, not assumed.

Key theorem with proof Intermediate+

The organising result is that the hierarchical posterior mean is exactly the conjugate posterior mean evaluated at the hyperparameter, and that empirical Bayes estimates that hyperparameter from the marginal, so the EB shrinkage factor is a data-driven version of the hierarchical one.

Theorem (the Gaussian hierarchical posterior and its empirical-Bayes estimate). Let $y_j\mid {\theta }_j\sim N({\theta }_j,{\sigma }^2)$ independently and ${\theta }_j\mid \mu ,{\tau }^2\sim N(\mu ,{\tau }^2)$ independently, $j=1,\dots ,J$, with ${\sigma }^2$ known.

(i) (Group posterior given the hyperparameter.) For fixed $(\mu ,{\tau }^2)$ the posterior of ${\theta }_j$ is $N(\mu +B(y_j-\mu ),(1-B){\sigma }^2)$ with $B={\tau }^2/({\sigma }^2+{\tau }^2)$, and the posterior mean is the partial-pooling estimate ${\hat{\theta }}_j=(1-B)\mu +B{y}_j$.

(ii) (Marginal.) Integrating out ${\theta }_j$, the marginal of $y_j$ given $(\mu ,{\tau }^2)$ is $N(\mu ,{\sigma }^2+{\tau }^2)$, with the $y_j$ independent across $j$.

(iii) (Empirical-Bayes shrinkage.) The marginal MLE is $\hat{\mu }=\bar{y}$ and $\widehat{{\sigma }^2+{\tau }^2}=\frac{1}{J}{\sum }_j(y_j-\bar{y}{)}^2$ (taking the positive part of $\widehat{{\tau }^2}$), so the EB shrinkage factor $\hat{B}=\widehat{{\tau }^2}/({\sigma }^2+\widehat{{\tau }^2})$ is estimated from the empirical spread of the $y_j$; the EB estimator is ${\hat{\theta }}_j^{EB}=\bar{y}+\hat{B}(y_j-\bar{y})$.

Proof of (i). With $\varphi =(\mu ,{\tau }^2)$ fixed, the joint of $({\theta }_j,y_j)$ is the conjugate normal-normal model of 45.03.02 with prior ${\theta }_j\sim N(\mu ,{\tau }^2)$ and likelihood $y_j\mid {\theta }_j\sim N({\theta }_j,{\sigma }^2)$. Working in precisions, the posterior precision is the sum of prior and data precisions, ${\tau }_{\text{post}}^{-2}={\tau }^{-2}+{\sigma }^{-2}$, and the posterior mean is the precision-weighted average

$${\hat{\theta }}_j=\frac{{\tau }^{-2}\mu +{\sigma }^{-2}{y}_j}{{\tau }^{-2}+{\sigma }^{-2}}=\frac{{\sigma }^2\mu +{\tau }^2{y}_j}{{\sigma }^2+{\tau }^2}=(1-B)\mu +B{y}_j,$$

with $B={\tau }^2/({\sigma }^2+{\tau }^2)$. The posterior variance is ${\tau }_{\text{post}}^2=({\tau }^{-2}+{\sigma }^{-2}{)}^{-1}={\sigma }^2{\tau }^2/({\sigma }^2+{\tau }^2)=(1-B){\sigma }^2$. $\square$

Proof of (ii). Write $y_j={\theta }_j+{\epsilon }_j$ with ${\theta }_j\sim N(\mu ,{\tau }^2)$ and ${\epsilon }_j\sim N(0,{\sigma }^2)$ independent. A sum of independent normals is normal with summed mean and variance, so $y_j\sim N(\mu ,{\tau }^2+{\sigma }^2)$. Independence across $j$ follows from the conditional-independence factorisation: $m_{\varphi }(y)={\prod }_j\int p(y_j\mid {\theta }_j)p({\theta }_j\mid \varphi )d{\theta }_j={\prod }_{j}N(y_j;\mu ,{\sigma }^2+{\tau }^2)$. $\square$

Proof of (iii). By part (ii) the $y_j$ are i.i.d. $N(\mu ,\nu )$ with $\nu :={\sigma }^2+{\tau }^2$. The likelihood is the ordinary normal likelihood in $(\mu ,\nu )$, whose MLEs are the sample mean $\hat{\mu }=\bar{y}$ and sample variance $\hat{\nu }=\frac{1}{J}{\sum }_j(y_j-\bar{y}{)}^2$. Inverting $\nu ={\sigma }^2+{\tau }^2$ gives $\widehat{{\tau }^2}=(\hat{\nu }-{\sigma }^2{)}^{+}$, the positive part guarding against a negative variance estimate when the groups cluster more tightly than sampling noise alone would predict. Substituting into $B$ yields $\hat{B}=\widehat{{\tau }^2}/\hat{\nu }$ and ${\hat{\theta }}_j^{EB}=\bar{y}+\hat{B}(y_j-\bar{y})$, the conjugate posterior mean of part (i) with $(\mu ,{\tau }^2)$ replaced by their marginal MLEs. $\square$

Bridge. This theorem builds toward the parametric and nonparametric empirical-Bayes constructions of the Advanced section and the variance correction that hierarchical Bayes supplies over the plug-in EB estimate, and it appears again in the James-Stein estimator 45.06.04, whose shrinkage factor $1-(J-2){\sigma }^2/\Vert y-\bar{y}{\Vert }^2$ is precisely an unbiased estimate of the EB factor $1-B$ derived here. This is exactly the device that turns a fixed conjugate prior into a data-driven one: the marginal $m_{\varphi }$ of 45.03.02 is the bridge between the parameter level and the population level, and maximising it over $\varphi$ lets all $J$ groups jointly determine the prior. The central insight is that hierarchical and empirical Bayes are the same model read at two levels of completeness — the hierarchy integrates the hyperparameter, EB point-estimates it — so partial pooling, shrinkage, and borrowing strength are one phenomenon: a group with little information is pulled toward a population center the other groups reveal. Putting these together, the partial-pooling estimate $(1-B)\mu +B{y}_j$ exhibits the shrinkage of 45.03.01 as a convex combination whose weight $B$ generalises the single-group conjugate weight by tying it to the cross-group spread, and the prior-precision dial ${\tau }^{-2}$ is dual to the ridge regularisation strength of 45.06.03, with EB choosing that dial automatically from the data.

Exercises Intermediate+

Advanced results Master

The hierarchy as a compound model and the marginal as the EB objective

The hierarchical model $\varphi \sim p(\varphi )$, ${\theta }_j\mid \varphi \sim p({\theta }_j\mid \varphi )$, $y_j\mid {\theta }_j\sim p(y_j\mid {\theta }_j)$ separates inference into two coupled problems: estimating the population level $\varphi$ and estimating the group level ${\theta }_j$. The coupling runs through the marginal $m_{\varphi }(y_j)=\int p(y_j\mid {\theta }_j)p({\theta }_j\mid \varphi )d{\theta }_j$, which is the only place the data touch $\varphi$ once $\theta$ is integrated out. Full hierarchical Bayes computes the posterior $p(\varphi \mid y)\propto p(\varphi ){\prod }_j{m}_{\varphi }(y_j)$ and then the group posteriors $p({\theta }_j\mid y)=\int p({\theta }_j\mid y,\varphi )p(\varphi \mid y)d\varphi$, averaging the conjugate group posterior of 45.03.02 over the hyperparameter posterior. Parametric empirical Bayes truncates this to the single value $\hat{\varphi }=\arg {\max }_{\varphi }{\prod }_j{m}_{\varphi }(y_j)$, the marginal MLE, and reports $p({\theta }_j\mid y,\hat{\varphi })$. The two coincide in the limit $J\to \infty$, where $p(\varphi \mid y)$ concentrates at $\hat{\varphi }$ and the average over $\varphi$ degenerates to evaluation at the mode. The marginal-MLE optimisation is itself a latent-variable problem — the ${\theta }_j$ are the latents — solved cleanly by the EM algorithm 45.08.07, whose E-step computes the conjugate group posteriors and whose M-step updates $\hat{\varphi }$ from the expected complete-data sufficient statistics; for the Gaussian hierarchy EM reproduces the closed-form $\hat{\mu }=\bar{y}$, $\widehat{{\tau }^2}=(\hat{\nu }-{\sigma }^2{)}^{+}$ of the Key theorem.

Parametric EB, the variance correction, and the cost of point-estimating $\varphi$

The defect of plug-in empirical Bayes is that it treats the estimated hyperparameter as known, so the reported group-posterior variance $(1-\hat{B}){\sigma }^2$ omits the contribution of the uncertainty in $\hat{\varphi }$. The honest variance, by the law of total variance, adds a term of order $1/J$ that records the randomness of $\hat{\mu }$ and ${\hat{\tau }}^2$ ^{[Morris 1983]}. Carl Morris' parametric-EB programme supplies the correction explicitly for the Gaussian hierarchy: the shrinkage factor is estimated by $\hat{B}=1-(J-3){\sigma }^2/S$ with $S={\sum }_j(y_j-\bar{y}{)}^2$ (the $(J-3)$ rather than $(J-2)$ accounting for the estimation of $\mu$), and the corrected posterior variance inflates $(1-\hat{B}){\sigma }^2$ by a data-based estimate of the $\varphi$-uncertainty term. This is the precise sense in which empirical Bayes is approximate hierarchical Bayes: the point estimate is asymptotically efficient but throws away the hyperparameter's posterior spread, and the correction restores the part that matters at finite $J$. The connection to shrinkage is exact — $\hat{B}$ is the James-Stein factor 45.06.04, and the EB derivation is the cleanest route to understanding why estimating many means jointly dominates estimating each alone.

Nonparametric (Robbins) empirical Bayes and $f$-modeling

Robbins' 1956 insight was that for many sampling models the posterior functional of interest depends on the mixing distribution $G$ only through the marginal $m$, so one may estimate $m$ from the data and never estimate $G$ at all ^{[Robbins 1956]}. For the Poisson model the posterior mean is $(y+1)m(y+1)/m(y)$, recoverable from marginal frequencies; for the normal model with $y\mid \theta \sim N(\theta ,{\sigma }^2)$ and arbitrary $G$, Tweedie's formula gives the posterior mean as $\mathbb{E}[\theta \mid y]=y+{\sigma }^2\frac{d}{dy}\log m(y)$, again a functional of the marginal alone. This is the $f$-modeling strategy — model the marginal density $f=m$ and differentiate, as opposed to $g$-modeling, which posits a parametric or flexible form for $G$ and deconvolves ^{[Efron 2010]}. The $f$-modeling route is robust because the marginal is directly observed, but it is delicate near the tails where $m$ is poorly estimated, and the derivative or ratio amplifies estimation error; modern large-scale-inference practice smooths the marginal (e.g. by a low-order exponential-family fit) before applying the formula. The deep point is that with thousands of parallel problems the marginal $m$ is estimated very accurately, so the empirical-Bayes correction becomes nearly as good as knowing the true prior — the regime Efron calls learning the prior from the ensemble.

Hierarchical generalisations and computation

Beyond the conjugate Gaussian and gamma-Poisson hierarchies, multilevel models stack arbitrarily many levels (groups within regions within countries), allow non-conjugate population distributions and link functions (hierarchical logistic and Poisson regression), and place structured hyperpriors (the half-Cauchy on $\tau$ is a standard weakly-informative default that avoids the improper-posterior hazards of the flat or $1/{\tau }^2$ choices). Outside the conjugate case the marginal $m_{\varphi }$ has no closed form and neither full hierarchical Bayes nor the marginal MLE is analytic; the resolution is Markov chain Monte Carlo 45.03.05, where Gibbs sampling exploits the conditional conjugacy exhibited in Exercise 8 — each full conditional of a conjugate-within-level hierarchy is a closed-form update — and Hamiltonian Monte Carlo handles the non-conjugate and highly correlated cases. The hierarchy reorganises a high-dimensional posterior into a sequence of low-dimensional conditional updates, which is precisely what makes Gibbs sampling natural for multilevel models. Empirical Bayes survives in this computational landscape as the fast approximation: where MCMC integrates $\varphi$, EB maximises over it, trading a known understatement of uncertainty for a substantial reduction in cost.

Synthesis. The central insight is that hierarchical and empirical Bayes are one model read at two levels of completeness: the hierarchy of 45.03.02's conjugate update repeated across groups, with a hyperprior on the population, integrated in full Bayes and point-estimated in EB, so partial pooling, shrinkage, and borrowing strength are a single phenomenon. The foundational reason the marginal $m_{\varphi }$ organises everything is that, after integrating out $\theta$, it is the only channel through which data inform the population level — full Bayes computes $p(\varphi \mid y)\propto p(\varphi ){\prod }_j{m}_{\varphi }(y_j)$, EB maximises ${\prod }_j{m}_{\varphi }(y_j)$, and Robbins estimates a posterior functional from $m$ alone — so the closed-form conjugate marginal of 45.03.02 is exactly the EB objective and the Robbins/Tweedie $f$-modeling input. Putting these together, the Gaussian partial-pooling estimate $(1-B)\mu +B{y}_j$ generalises the single-group conjugate posterior mean by tying its weight to the cross-group spread, its EB factor $\hat{B}$ is dual to the James-Stein factor 45.06.04 and to the ridge penalty of 45.06.03 under prior-precision-equals-regularisation, and the understatement of variance EB incurs appears again in the Morris correction and is resolved by integrating $\varphi$. The bridge is computation: the marginal MLE is an EM problem 45.08.07, the conditional conjugacy of the hierarchy is what a Gibbs sampler cycles, and the boundary where closed forms end — non-conjugate or deeply nested hierarchies — is where MCMC 45.03.05 takes over.

Full proof set Master

Proposition 1 (Gaussian hierarchical posterior and partial pooling). For $y_j\mid {\theta }_j\sim N({\theta }_j,{\sigma }^2)$ and ${\theta }_j\mid \mu ,{\tau }^2\sim N(\mu ,{\tau }^2)$ with $(\mu ,{\tau }^2)$ fixed and ${\sigma }^2$ known, ${\theta }_j\mid y,\mu ,{\tau }^2\sim N(\mu +B(y_j-\mu ),(1-B){\sigma }^2)$ with $B={\tau }^2/({\sigma }^2+{\tau }^2)$.

Proof. By the conditional-independence factorisation, $p({\theta }_j\mid y,\mu ,{\tau }^2)\propto p(y_j\mid {\theta }_j)p({\theta }_j\mid \mu ,{\tau }^2)\propto \exp (-\frac{(y_j-{\theta }_j{)}^2}{2{\sigma }^2}-\frac{({\theta }_j-\mu {)}^2}{2{\tau }^2})$. Collecting the quadratic in ${\theta }_j$, the coefficient of ${\theta }_j^2$ is $\frac{1}{2}({\sigma }^{-2}+{\tau }^{-2})$ and the coefficient of ${\theta }_j$ is ${\sigma }^{-2}{y}_j+{\tau }^{-2}\mu$, so the posterior is normal with precision ${\sigma }^{-2}+{\tau }^{-2}$ and mean $({\sigma }^{-2}{y}_j+{\tau }^{-2}\mu )/({\sigma }^{-2}+{\tau }^{-2})$. Multiplying numerator and denominator by ${\sigma }^2{\tau }^2$ gives mean $({\tau }^2{y}_j+{\sigma }^2\mu )/({\sigma }^2+{\tau }^2)=\mu +B(y_j-\mu )$ and variance $({\sigma }^{-2}+{\tau }^{-2}{)}^{-1}={\sigma }^2{\tau }^2/({\sigma }^2+{\tau }^2)=(1-B){\sigma }^2$. $\square$

Proposition 2 (closed-form marginals: normal and negative-binomial). In the Gaussian hierarchy $y_j\mid \mu ,{\tau }^2\sim N(\mu ,{\sigma }^2+{\tau }^2)$; in the gamma-Poisson hierarchy $y_j\mid \alpha ,\beta$ is negative binomial with size $\alpha$ and success probability $\beta /(\beta +1)$.

Proof. For the Gaussian case $y_j={\theta }_j+{\epsilon }_j$ is a sum of independent $N(\mu ,{\tau }^2)$ and $N(0,{\sigma }^2)$ variables, hence $N(\mu ,{\sigma }^2+{\tau }^2)$ by convolution of normals. For the gamma-Poisson case,

$$m(y)={\int }_0^{\infty }\frac{e^{-\lambda }{\lambda }^y}{y!}\frac{{\beta }^{\alpha }}{\Gamma (\alpha )}{\lambda }^{\alpha -1}{e}^{-\beta \lambda }d\lambda =\frac{{\beta }^{\alpha }}{y!\Gamma (\alpha )}\frac{\Gamma (y+\alpha )}{(\beta +1{)}^{y+\alpha }}=\left(\genfrac{}{}{0ex}{}{y+\alpha -1}{y}\right)\left(\frac{\beta }{\beta +1}{)}^{\alpha }\right(\frac{1}{\beta +1}{)}^y,$$

using ${\int }_0^{\infty }{\lambda }^{y+\alpha -1}{e}^{-(\beta +1)\lambda }d\lambda =\Gamma (y+\alpha )/(\beta +1{)}^{y+\alpha }$ and $\Gamma (y+\alpha )/(y!\Gamma (\alpha ))=\left(\genfrac{}{}{0ex}{}{y+\alpha -1}{y}\right)$ for the generalised binomial coefficient. This is the negative binomial. $\square$

Proposition 3 (Robbins' formula and Tweedie's formula). For $y\mid \lambda \sim \text{Poisson}(\lambda )$, $\lambda \sim G$, the posterior mean is $\mathbb{E}[\lambda \mid y=k]=(k+1)m(k+1)/m(k)$. For $y\mid \theta \sim N(\theta ,{\sigma }^2)$, $\theta \sim G$, the posterior mean is $\mathbb{E}[\theta \mid y]=y+{\sigma }^2(\log m{)}^{\prime }(y)$.

Proof. Poisson: with $m(k)=\int \frac{e^{-\lambda }{\lambda }^k}{k!}dG(\lambda )$, the posterior mean is $\int \lambda \frac{e^{-\lambda }{\lambda }^k}{k!}dG/m(k)$. The numerator equals $(k+1)\int \frac{e^{-\lambda }{\lambda }^{k+1}}{(k+1)!}dG=(k+1)m(k+1)$, giving the ratio. Normal: write $m(y)=\int {\phi }_{\sigma }(y-\theta )dG(\theta )$ with ${\phi }_{\sigma }$ the $N(0,{\sigma }^2)$ density. Differentiating under the integral, $m^{\prime }(y)=\int {\phi }_{\sigma }^{\prime }(y-\theta )dG(\theta )=\int -\frac{y-\theta }{{\sigma }^2}{\phi }_{\sigma }(y-\theta )dG(\theta )$. Hence ${\sigma }^2{m}^{\prime }(y)=\int (\theta -y){\phi }_{\sigma }(y-\theta )dG(\theta )=(\mathbb{E}[\theta \mid y]-y)m(y)$, where the last step uses that the posterior density of $\theta$ given $y$ is ${\phi }_{\sigma }(y-\theta )dG(\theta )/m(y)$. Rearranging, $\mathbb{E}[\theta \mid y]=y+{\sigma }^2{m}^{\prime }(y)/m(y)=y+{\sigma }^2(\log m{)}^{\prime }(y)$. $\square$

Proposition 4 (asymptotic agreement of EB and hierarchical Bayes). In the Gaussian hierarchy with a fixed proper hyperprior $p(\varphi )$ and $\hat{\varphi }$ the marginal MLE, the hierarchical-Bayes group posterior mean $\mathbb{E}[{\theta }_j\mid y]$ and the empirical-Bayes mean $\mathbb{E}[{\theta }_j\mid y,\hat{\varphi }]$ differ by $O_p(1/J)$ as $J\to \infty$.

Proof (sketch). The hierarchical posterior mean is $\mathbb{E}[{\theta }_j\mid y]=\int \mathbb{E}[{\theta }_j\mid y,\varphi ]p(\varphi \mid y)d\varphi$. By the Bernstein-von Mises theorem applied to the marginal likelihood ${\prod }_j{m}_{\varphi }(y_j)$ — a regular i.i.d. model in $\varphi$ with $J$ observations — the posterior $p(\varphi \mid y)$ is asymptotically $N(\hat{\varphi },J^{-1}I(\hat{\varphi }{)}^{-1})$, concentrating at $\hat{\varphi }$ at rate $J^{-1/2}$. Expanding the smooth map $\varphi \mapsto \mathbb{E}[{\theta }_j\mid y,\varphi ]$ around $\hat{\varphi }$ in a second-order Taylor series and integrating against this posterior, the first-order term vanishes (it multiplies the posterior-mean deviation $\mathbb{E}[\varphi \mid y]-\hat{\varphi }=O_p(1/J)$) and the second-order term contributes the curvature against the posterior covariance $O(1/J)$. Hence $\mathbb{E}[{\theta }_j\mid y]=\mathbb{E}[{\theta }_j\mid y,\hat{\varphi }]+O_p(1/J)$. The discrepancy is exactly the hyperparameter-uncertainty term the plug-in EB variance omits, of the same order. $\square$

Proposition 5 (improper hyperprior can give improper posterior). For the Gaussian hierarchy with $J$ groups and the hyperprior $p(\mu ,{\tau }^2)\propto 1/{\tau }^2$ on $\mu \in \mathbb{R}$, ${\tau }^2>0$, the joint posterior is improper.

Proof. Integrate the joint posterior over $(\theta ,\mu )$ to obtain the marginal posterior of ${\tau }^2$, proportional to $p({\tau }^2)({\sigma }^2+{\tau }^2{)}^{-(J-1)/2}\exp (-S/[2({\sigma }^2+{\tau }^2)])$ up to the $\mu$-integration constant, where $S={\sum }_j(y_j-\bar{y}{)}^2$; this follows from Proposition 2 (the marginal of the $y_j$ is $N(\mu ,{\sigma }^2+{\tau }^2)$) and integrating the flat $\mu$. With $p({\tau }^2)\propto 1/{\tau }^2$ the integrand behaves like ${\tau }^{-2}\cdot ({\sigma }^2{)}^{-(J-1)/2}\exp (-S/2{\sigma }^2)$ as ${\tau }^2\to 0^{+}$, and ${\int }_0{\tau }^{-2}d{\tau }^2$ diverges at the origin. The marginal posterior of ${\tau }^2$ is therefore not integrable near ${\tau }^2=0$, so the joint posterior is improper. The flat prior $p({\tau }^2)\propto 1$ avoids this particular divergence at the origin but must still be checked at infinity; the half-Cauchy on $\tau$ is proper at both ends. $\square$

Connections Master

• Conjugate priors and exponential-family Bayesian inference 45.03.02 is the engine each level of the hierarchy runs: the group posterior $N(\mu +B(y_j-\mu ),(1-B){\sigma }^2)$ is the conjugate normal-normal update with $(\mu ,{\tau }^2)$ as the prior, the closed-form marginal $m_{\varphi }$ that couples the levels is the conjugate marginal likelihood proved there, and the conditional conjugacy a Gibbs sampler exploits is that same update applied block by block. Hierarchical Bayes is conjugate Bayes repeated across groups with a hyperprior on the shared parameters.

• Bayes estimation under loss 45.03.01 supplies the decision-theoretic reading: the partial-pooling estimate is the squared-error Bayes estimator of ${\theta }_j$ under the hierarchical prior, and the EB estimator is the feasible approximation that estimates that prior from the ensemble, so shrinkage toward the population center is the optimal compromise between the group's datum and the borrowed crowd information.

• The James-Stein estimator and inadmissibility of the MLE 45.06.04 is the same shrinkage seen from frequentist decision theory: the James-Stein factor $1-(J-2){\sigma }^2/\Vert y{\Vert }^2$ is an unbiased estimate of the hierarchical shrinkage $1-B$, so James-Stein is the empirical-Bayes plug-in of the Gaussian hierarchy, and the dimension threshold $J\ge 3$ is the condition under which the cross-group pooling strictly improves on estimating each mean alone.

• The EM algorithm 45.08.07 solves the marginal-likelihood maximisation that defines parametric empirical Bayes: with the group parameters ${\theta }_j$ as latent variables, the E-step computes the conjugate group posteriors and the M-step updates the hyperparameter $\hat{\varphi }$ from the expected sufficient statistics, reproducing the closed-form $\hat{\mu }=\bar{y}$, $\widehat{{\tau }^2}=(\hat{\nu }-{\sigma }^2{)}^{+}$ in the Gaussian case.

• Markov chain Monte Carlo 45.03.05 is the resolution when the hierarchy is non-conjugate: full hierarchical Bayes integrates the hyperparameter $\varphi$ rather than point-estimating it, and outside the conjugate case that integration is performed by Gibbs sampling — cycling the conditional-conjugate full conditionals derived here — or by Hamiltonian Monte Carlo for the correlated and non-conjugate levels.

Historical & philosophical context Master

The empirical-Bayes idea was introduced by Herbert Robbins at the Third Berkeley Symposium in 1956 ^{[Robbins 1956]}, who observed that in a sequence of parallel estimation problems the prior could be estimated from the accumulated data without ever being assumed, deriving the formula $\mathbb{E}[\lambda \mid y]=(y+1)m(y+1)/m(y)$ for the Poisson mean from marginal frequencies alone. The connection to the Stein phenomenon was made by Bradley Efron and Carl Morris in a sequence of papers in the early 1970s, who showed the James-Stein estimator is the empirical-Bayes rule for the Gaussian hierarchy with the prior variance estimated from the ensemble; Carl Morris' 1983 Journal of the American Statistical Association paper ^{[Morris 1983]} gave the systematic parametric-EB theory, including the variance correction that accounts for the uncertainty in the estimated hyperparameter and so closes the honesty gap between empirical and full Bayes.

The full hierarchical (multilevel) treatment, in which the hyperparameter carries its own prior and is integrated rather than point-estimated, was developed in the Bayesian decision-theory tradition of I. J. Good, who coined the term hierarchical Bayes, and Dennis Lindley and Adrian Smith's 1972 Journal of the Royal Statistical Society B paper on Bayes estimates for the linear model, which derived the exchangeable normal hierarchy and its shrinkage estimators. The modern computational and applied synthesis is in Andrew Gelman and colleagues' Bayesian Data Analysis ^{[Gelman 2013]}, whose eight-schools example became the canonical illustration of partial pooling, and Bradley Efron's Large-Scale Inference ^{[Efron 2010]} reframed empirical Bayes for the thousands-of-parallel-problems regime of genomics and high-throughput testing, distinguishing the $f$-modeling and $g$-modeling routes and establishing that with enough parallel problems the estimated prior is nearly as good as the true one.

Bibliography Master

@inproceedings{robbins1956empirical,
  author    = {Robbins, Herbert},
  title     = {An Empirical Bayes Approach to Statistics},
  booktitle = {Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability},
  volume    = {1},
  pages     = {157--163},
  publisher = {University of California Press},
  year      = {1956}
}

@article{lindleysmith1972bayes,
  author  = {Lindley, Dennis V. and Smith, Adrian F. M.},
  title   = {Bayes Estimates for the Linear Model},
  journal = {Journal of the Royal Statistical Society. Series B},
  volume  = {34},
  number  = {1},
  pages   = {1--41},
  year    = {1972}
}

@article{efronmorris1973stein,
  author  = {Efron, Bradley and Morris, Carl},
  title   = {Stein's Estimation Rule and Its Competitors --- An Empirical Bayes Approach},
  journal = {Journal of the American Statistical Association},
  volume  = {68},
  number  = {341},
  pages   = {117--130},
  year    = {1973}
}

@article{morris1983parametric,
  author  = {Morris, Carl N.},
  title   = {Parametric Empirical Bayes Inference: Theory and Applications},
  journal = {Journal of the American Statistical Association},
  volume  = {78},
  number  = {381},
  pages   = {47--65},
  year    = {1983}
}

@book{efron2010large,
  author    = {Efron, Bradley},
  title     = {Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction},
  series    = {Institute of Mathematical Statistics Monographs},
  publisher = {Cambridge University Press},
  year      = {2010}
}

@book{gelman2013bayesian,
  author    = {Gelman, Andrew and Carlin, John B. and Stern, Hal S. and Dunson, David B. and Vehtari, Aki and Rubin, Donald B.},
  title     = {Bayesian Data Analysis},
  edition   = {3rd},
  publisher = {CRC Press},
  year      = {2013}
}

@book{casellaberger2002statistical,
  author    = {Casella, George and Berger, Roger L.},
  title     = {Statistical Inference},
  edition   = {2nd},
  publisher = {Duxbury},
  year      = {2002}
}