Bayesian Hierarchical Model

Definition

A Bayesian hierarchical model organizes uncertainty into two or more levels: observable data governed by parameters, and those parameters governed by hyperparameters with their own prior distribution. Each level borrows information from the others through the joint posterior, enabling partial pooling — estimates for sparse groups are shrunk toward a common mean while rich groups remain data-driven.

Key Ideas

Partial pooling (Bühlmann credibility): the posterior estimate for group $i$ is a weighted average of the group's own data and the grand mean — the weight determined by the signal-to-noise ratio at each level.
Conjugate structures: Normal/Normal and Poisson/Gamma hierarchies yield closed-form Gibbs conditionals; non-conjugate cases require Griddy-Gibbs or Metropolis steps.
Proper priors at every level guarantee a proper joint posterior and valid Gibbs convergence; improper priors at the top level may still yield proper posteriors when data are informative.
Predictive distributions: $f(\tilde{Y} \mid Y) = \int f(\tilde{Y} \mid \theta)\, f(\theta \mid Y)\, d\theta$ — automatically propagates parameter uncertainty into future predictions; estimated by averaging $f(\tilde{Y} \mid \theta^{(i)})$ over posterior draws.

How It Works

Three-Level Structure

Let $Y_{ij}$ be the $j$ th observation in group $i$ :

$\text{Level 1 (data model):} \quad Y_{ij} \mid \theta_i \sim f(y \mid \theta_i)$ $\text{Level 2 (latent parameters):} \quad \theta_i \mid \mu, \tau \sim g(\theta \mid \mu, \tau)$ $\text{Level 3 (hyperprior):} \quad \mu, \tau \sim h(\mu, \tau)$

Bayes' theorem ties the levels together:

$p(\theta_1,\ldots,\theta_n, \mu, \tau \mid Y) \propto h(\mu,\tau) \prod_i g(\theta_i|\mu,\tau) \prod_{i,j} f(y_{ij}|\theta_i)$

Three-Component Marketing Framework and Heterogeneity Distributions (Rossi-Allenby 2003)

Rossi and Allenby articulate that every Bayesian marketing model has three coherent components: (1) a within-unit likelihood $p(y_i|\theta_i)$ modelling individual consumer behavior, (2) a cross-unit heterogeneity distribution $p(\theta_i|\tau)$ , and (3) an action rule derived from a loss function applied to the posterior.

Normal hierarchical model. The canonical three-level specification (where $W$ denotes the Wishart distribution):

$\theta_i \sim N(\bar\theta, V_\theta), \quad \bar\theta \sim N(\bar{\bar\theta}, A^{-1}), \quad V_\theta^{-1} \sim W(\nu, V)$

Two-block Gibbs: draw $\{\theta_i|\tau, y_i\}$ (conjugate given the likelihood) and $\tau|\{\theta_i\}$ (conjugate from the Wishart-Normal hyperprior). Shrinkage toward $\bar\theta$ is strong when within-unit data are sparse and between-unit variation is small.

Observable heterogeneity. When unit-level covariates $z_i$ are available:

$\theta_i = Bz_i + u_i, \quad u_i \sim N(0, V_\theta)$

$B$ and $V_\theta$ are estimated jointly; residual heterogeneity $u_i$ is the unexplained part after covariate adjustment.

Mixture-of-normals heterogeneity. To accommodate multimodal or heavy-tailed heterogeneity:

$p(\theta_i|\tau) = \sum_{k=1}^K r_k\, \phi(\theta_i|\bar\theta_k, V_k)$

A normal mixture avoids the convex-hull limitation of finite-mixture models (in which predictions are constrained to the convex hull of the component means). Component assignments are auxiliary variables in the Gibbs sampler; the number of components $K$ is selected by marginal likelihood, not BIC (Bayesian Information Criterion, which is "extremely inaccurate" for correlated-parameter models).

Structural heterogeneity. When individuals switch between behavioral regimes within the panel:

$p(y_{it}|\{\theta_{ik}\}) = \sum_{k=1}^K r_k\, p_k(y_{it}|\theta_{ik})$

The mixture operates at the observation level, not the parameter level — each $\theta_{ik}$ is a complete behavioral type, and the assignment probability $r_k$ governs regime prevalence.

Diagnostic check (Allenby-Rossi 1999). Compare unit posterior means $\hat\theta_i$ to the marginal predictive:

$p(\theta_{i'}|\text{data}) = \iint \phi(\theta|\bar\theta, V_\theta)\, p(\bar\theta, V_\theta|\text{data})\, d\bar\theta\, dV_\theta$

Unit posterior means that are outliers relative to this predictive distribution signal that the heterogeneity model is too restrictive. Properly accounting for hyperparameter uncertainty (integrating over $\bar\theta$ and $V_\theta$ ) widens the predictive relative to a plug-in check.

Bayesian decision theory. The optimal action is:

$a^* = \arg\min_a \int l(a, \theta)\, p(\theta|\text{data})\, d\theta$

Evaluating profit only at the posterior mean (plug-in) overstates expected profit because it ignores parameter uncertainty — the loss function is typically nonlinear in $\theta$ . Full posterior integration is required for valid targeting decisions. The value of disaggregate (unit-level) information is quantified by $\Pi_{\text{disagg}} - \Pi_{\text{agg}}$ where $\Pi_{\text{disagg}}$ uses individual posterior $\theta_i$ and $\Pi_{\text{agg}}$ uses a pooled estimate.

Workers Compensation Credibility (Scollnik ~1996, Example 4)

133 occupation classes, 7 years of claim frequency data; Klugman (1992) Data Set 2. Let $Y_{ij}$ = claim count, $P_{ij}$ = scaled exposure, $x_{ij} = Y_{ij}/P_{ij}$ = relative frequency.

HNM (Hierarchical Normal Model):

$x_{ij} \mid \theta_i, \tau_1 \sim \mathcal{N}(\theta_i,\, P_{ij}/\tau_1), \quad \theta_i \mid \mu, \tau_2 \sim \mathcal{N}(\mu,\, 1/\tau_2)$

Diffuse proper priors: $\mu \sim \mathcal{N}(0, 0.001)$ ; $\tau_1, \tau_2 \sim \text{Gamma}(0.001, 0.001)$ . Full conditionals are Normal/Gamma — closed-form Gibbs. Drawback: allows negative frequencies.

HFLPM (Hierarchical First-Level Poisson Model):

$Y_{ij} \mid \lambda_i \sim \text{Poisson}(\lambda_i P_{ij}), \quad \lambda_i \mid \mu, \tau \sim \mathcal{N}(\mu, 1/\tau)$

Non-negative counts by construction. Best predictive performance on the hold-out year: overall measure of prediction success (OMPS) posterior mean = 12.96 (sd = 0.64).

VSHNM (Variance-Stabilized HNM):

Apply $z_{ij} = 2\sqrt{x_{ij}}$ before fitting a Normal hierarchy. Stabilizes the variance of Poisson-like data. OMPS = 13.72 — better than HNM but worse than HFLPM.

Gibbs Sampling for Hierarchical Models

The block structure of the posterior corresponds naturally to Gibbs blocks. For the HNM:

$\theta_i \mid \mu, \tau_1, \tau_2, Y \sim \mathcal{N}$ — weighted average of group data and prior mean
$\mu \mid \theta_{1:n}, \tau_2 \sim \mathcal{N}$ — grand mean update
$\tau_1 \mid \theta_{1:n}, Y \sim \text{Gamma}$ — within-group precision
$\tau_2 \mid \theta_{1:n}, \mu \sim \text{Gamma}$ — between-group precision

All four blocks have closed-form conjugate conditionals. BUGS (Bayesian Inference Using Gibbs Sampling) implements this automatically from the model declaration.

Gamma-Poisson Model for Geographic Mortality Data (Tsutakawa 1988)

For count data $y_{ij}$ (deaths in region $i$ , demographic group $j$ ) with exposure $n_{ij}$ , the model adds two random-effect layers above the Poisson observation:

$y_{ij} \mid p_{ij} \sim \text{Poisson}(n_{ij} p_{ij})$ $p_{ij} \mid \theta_j, \sigma_j^2, z_i \sim \text{Gamma}(\sigma_j^{-2},\; \sigma_j^2 \theta_j z_i)$ $z_i \sim \text{Inv-Gamma-1}(\gamma), \quad E(z_i) = 1$

Marginalising over $p_{ij}$ gives a negative binomial for $y_{ij}$ with variance $= n_{ij}\theta_j z_i + \sigma_j^2 (n_{ij}\theta_j z_i)^2$ — Poisson plus extra-Poisson. The regional effect $z_i$ is shared across all demographic groups $j$ in region $i$ , inducing positive cross-age correlation within a county.

Posterior mean (partial pooling): $E(p_{ij} \mid y_{ij}, \xi, z_i)$ is a weighted average of the observed rate $y_{ij}/n_{ij}$ (weight $\propto n_{ij}\bar{p}_{ij}$ ) and the expected rate $\theta_j z_i$ (weight $\propto \sigma_j^{-2}$ ). Small counties shrink heavily toward expected; large counties stay data-driven.

Estimation: Empirical Bayes via profile likelihood — estimate hyperparameters $\xi = (\theta, \sigma^2, \gamma)$ by ML (profile over $\gamma$ , integrate over $z_i$ ), then estimate rates by posterior means given $\hat{\xi}$ . Missouri lung cancer: 115 counties × 4 age groups, 1972–1981; likelihood-ratio (LR) test vs. fixed-effects Poisson = 51.4; Empirical Bayes (EB) standard deviations up to 15% smaller than raw rates.

Hierarchical Pooling of Air Pollution–Mortality Evidence (Dominici-Samet-Zeger 2000)

The NMMAPS (National Morbidity, Mortality, and Air Pollution Study) exemplifies the two-stage approach where Stage 1 delivers sufficient statistics and Stage 2 applies the hierarchical model to those statistics rather than to raw data.

Stage 1 — City-specific Poisson generalized additive model (GAM). For each of 20 US cities, city $c$ , age group $a$ , and day $t$ :

$\log(\mu^{ac}_t) = x^{ac\prime}_t \beta_c + \gamma_c \text{DOW}_t + \sum_{k=1}^{5} S^c_k(\cdot) + \delta^c_a + \sum_a S^{ac}(\text{time}_t,8/\text{yr})$

where $\text{DOW}_t$ is a day-of-week indicator, $S^c_k$ are natural splines for long-run time trend (7 df/yr), same-day temperature (6 df), 1–3 day lagged temperature (6 df), same-day dew point (3 df), lagged dew point (3 df), and $x^{ac}_t$ = (PM₁₀, O₃) with lags 0, 1, 2. Output per city: $\hat{\beta}_c \in \mathbb{R}^p$ and sample covariance $V_c$ . The large sample ( $T = 550$ – $2550$ days) justifies the normal approximation $\hat{\beta}_c \mid \beta_c \sim N_p(\beta_c, V_c)$ .

Stage 2 — Bayesian hierarchical pooling. Base-line independence model:

$\hat{\beta}_c \mid \beta_c \sim N_p(\beta_c, V_c)$ $\beta_c = z^\prime_c \alpha + \varepsilon_c, \quad \varepsilon_c \mid \Sigma \sim N_p(0, \Sigma)$ $\alpha \sim N(0, 100I), \quad \Sigma \sim \mathrm{IW}_p(df=p+1=3,\; D=3I) \quad \text{(IW = Inverse-Wishart)}$

where $z_c = (1, P^c_\text{poverty}, P^c_{>65}, \bar{X}^c_{\text{PM}_{10}})$ . All full conditionals are conjugate (normal for $\alpha$ and $\beta_c$ , IW for $\Sigma$ ): a pure Gibbs sampler runs 3000 iterations with burn-in 100 and thinning by 5. Convergence via Raftery-Lewis (1992) CODA ( $N_\text{min} \approx 2000$ ).

Spatial extension. City correlations follow an exponential decay model:

$\text{corr}(\beta_c, \beta_{c'}) = \exp\{-\theta\, d(c,c')\}$

The decay parameter $\theta \sim \text{LogNormal}(0.2, 0.5^2)$ lacks a conjugate conditional; a Metropolis-Hastings (M-H) step with gamma proposal (acceptance 0.3–0.5) is required. A regional step-function alternative recovers closed-form conditionals: $\tau^2 \sim IG(5, 8.5)$ (IG = Inverse-Gamma).

Key results. Univariate base-line, lag 0: +10 μg/m³ PM₁₀ → +0.48% daily mortality (95% credible interval [CI]: 0.05–0.92), Pr(α₀ > 0) = 0.98. Bivariate (PM₁₀ + O₃): +0.52% (0.16–0.85). Between-city standard deviation $\hat{\sigma} \approx 0.76$ (0.41–1.37), approximately twice the pooled effect — substantial heterogeneity. City-level covariates do not explain it. Robust across 27 sensitivity scenarios (3 models × 3 lags × 3 priors).

Limitation. The IW prior with only $p+1=3$ degrees of freedom is influential relative to 20 cities. The profile likelihood for $\sigma$ has 95% mass approximately on $(0, 0.4)$ , narrower than the IW posterior's $(0.41, 1.37)$ — a standard small- $n$ problem for IW hyperpriors.

Level-Specific R² and Pooling Factor λ (Gelman-Pardoe 2004)

For a multilevel model with level-m units $k = 1,\ldots,K_m$ , let $\theta_k$ be the quantity of interest and $\varepsilon_k = \theta_k - \mu_k$ the residual after the linear predictor $\mu_k$ . Two Markov Chain Monte Carlo (MCMC)-computed summary measures (posterior expectations of variance ratios) are:

Level-m $R^2$ (eq. 7): $R^2 = 1 - \frac{E[V_k(\varepsilon_k)]}{E[V_k(\theta_k)]}$

Measures the fraction of level-m variance explained by the predictor $\mu$ . $V_k$ denotes variance across units $k$ . Reduces to classical ordinary least squares (OLS) R² with no pooling (Xu 2003).

Pooling factor λ (eq. 8): $\lambda = 1 - \frac{V_k(E(\varepsilon_k))}{E[V_k(\varepsilon_k)]}$

Measures the fraction of residual variance reduced by the hierarchical model for the residuals $\varepsilon_k$ . $\lambda \approx 0$ : each group is estimated independently (no pooling). $\lambda \approx 1$ : all residual structure is captured by the group-level model.

Parallel structure: R² and λ have the same denominator structure. R² measures predictor fit to $\theta_k$ ; λ measures hierarchical shrinkage of $\varepsilon_k$ .

Bayesian adjusted $R^2$ for linear regression ( $n$ obs., $p$ predictors): $R^2_{\text{adj,Bayes}} = 1 - \frac{n-3}{n-p-2} \cdot \frac{y^T(I-H)y}{y^T I_c y}$

Strictly lower than the classical Wherry (1931) adjusted $R^2$ (factor $(n-1)/(n-p)$ ). For the balanced one-way hierarchical model, λ reduces to the Bühlmann credibility weight $\omega = 1 - \hat\sigma^2_\alpha / (\hat\sigma^2_\alpha + \hat\sigma^2_y/n)$ .

Minnesota radon example (919 houses, J=85 counties; Model 4 with uranium predictor):

Level	$R^2$	$\lambda$
Household ( $y$ )	0.21	0.03
County intercept ( $\alpha$ )	0.53	0.81
County slope ( $\beta$ )	0.83	0.97

Uranium explains 83% of cross-county slope variation; hierarchical pooling accounts for 97% of slope residual variance (slopes are nearly entirely determined by the county-level uranium model).

Both measures are free MCMC byproducts (no null model needed) and can be negative (misspecification signal).

Varying-Slope Multilevel Logistic Regression: Rich-State/Poor-State Voting Paradox (Gelman-Shor-Bafumi-Park 2007)

The apparent contradiction — rich states vote Democratic, rich individuals vote Republican — is resolved by allowing the income slope to vary across states and to depend on average state income.

Varying-intercept model (eq. 1): Constant slope $\beta$ , state-specific intercepts:

$\Pr(y_i = 1) = \text{logit}^{-1}(\alpha_{s[i]} + \beta x_i)$

Varying-intercept, varying-slope model (eq. 2): Both intercept and slope vary by state, each modeled as a linear function of average state income $u_s$ :

$\Pr(y_i = 1) = \text{logit}^{-1}(\alpha_{s[i]} + \beta_{s[i]} x_i)$ $\alpha_s = a^\alpha + b^\alpha u_s + \varepsilon^\alpha_s, \quad \beta_s = a^\beta + b^\beta u_s + \varepsilon^\beta_s$

Key finding: $\hat b^\beta < 0$ — the income–Republican slope is steepest in poor states (Mississippi) and near-zero in rich states (Connecticut). The aggregate relationship runs in the opposite direction because richer states have higher $\alpha_s$ (baseline Democratic lean), which dominates the cross-state comparison even as $\beta_s > 0$ within each state.

County-level model (§2.4): $y_c = \alpha_{s[c]} + \beta_{s[c]} x_c + \varepsilon_c$ fits county average income and vote shares with the same hierarchical structure — confirms individual-level inference from aggregated data.

Pattern emergence: Essentially absent 1968–1988; emerged post-1992 (Figures 14–16). The geographic income–vote correlation is a post-Cold-War phenomenon. Estimation via WinBUGS; model checking via binned residual plots.

Random-Effects Logit for Conjoint Analysis (Allenby-Rossi-McCulloch 2005, eqs. 16–19)

The standard hierarchical Bayes (HB) conjoint model applies the three-level hierarchy to a logit observation layer:

$\Pr_h(i) = \frac{\exp(x_i'\beta_h)}{\sum_j \exp(x_j'\beta_h)}, \quad \beta_h = \Gamma z_h + \xi_h, \quad \xi_h \sim \mathcal{N}(0, V_\beta)$

Priors: $\Gamma \sim \mathcal{N}(a, A)$ ; $V_\beta^{-1} \sim \text{Wishart}(w, W)$ . Gibbs cycle — (1) each $\beta_h$ one respondent at a time (Metropolis for non-conjugate logit likelihood); (2) $\Gamma|\{\beta_h\}, V_\beta \sim \mathcal{N}$ ; (3) $V_\beta|\{\beta_h\}, \Gamma \sim \text{IW}$ ; (4) repeat.

Credit card application: N=946 bank customers, 14,799 paired comparisons, 7 attributes; 20,000 iterations, convergence by ~6,000. Key finding: comparing low annual fee vs. low fixed interest as incentives to overcome the −3.614 out-of-state penalty, despite low interest having higher mean utility, the low-fee incentive produces 7.5% of respondents above a net-utility profitability threshold vs. 4.5% for low interest — nearly 2×. The explanation is the off-diagonal of $\hat V_\beta$ : $\text{Cov}(\text{low fee, out-of-state}) = 8.5$ (ρ=0.80) vs. $\text{Cov}(\text{low interest, out-of-state}) = 1.9$ (ρ=0.23). Respondents who are fee-sensitive also tend to be less penalized by the out-of-state attribute.

Blocked MCMC for Longitudinal Models (Chib-Carlin 1999)

The standard Gaussian linear mixed model $y_i = X_i\beta + W_i b_i + \varepsilon_i$ , $b_i\sim\mathcal{N}_q(0,D)$ , $\varepsilon_i\sim\mathcal{N}(0,\sigma^2 I)$ with Wishart prior on $D^{-1}$ and Gamma prior on $\sigma^{-2}$ , admits progressively blocked Gibbs samplers that dramatically reduce autocorrelation relative to the standard 4-block algorithm (Gelfand-Smith 1990).

Algorithm 2 (partial blocking): Marginalize $y_i$ over $b_i$ : $y_i|\beta,\sigma^2,D\sim\mathcal{N}_{n_i}(X_i\beta, V_i)$ where $V_i=\sigma^2I+W_iDW_i'$ . This gives a closed-form $\beta|y,\sigma^2,D\sim\mathcal{N}_p(\hat\beta,B)$ that does not condition on the random effects. After sampling β, draw $b_i|\beta,\sigma^2,D$ from its standard conditional. Autocorrelation time $\kappa(\beta)$ drops to 1.0 (essentially iid) in the HIV CD4 example.

Algorithm 3 (single block): The joint marginal $\pi(\sigma^2,D^{-1}|y)$ has no closed form, but the integrated likelihood $f(y|\sigma^2,D)$ can be evaluated via the Chib (1995) identity as a ratio of three tractable Normal densities: $f(y|\sigma^2,D) = \frac{\phi_p(\beta^*|\beta_0,B_0)\prod_i\phi_{n_i}(y_i|X_i\beta^*,V_i)}{\phi_p(\beta^*|\hat\beta,B)}$ evaluated at the Algorithm 2 posterior mean $\beta^*$ . Update $\theta=(\sigma^2,\psi)$ jointly by M-H with a multivariate-t proposal calibrated from a short Algorithm 2 warmup run. Reduces $\kappa(D_{21})$ from 60 → 11.

Binary longitudinal extension: For binary responses with probit link, data augmentation (Albert-Chib 1993) introduces latent $z_{it}|b_i\sim\mathcal{N}(x_{it}'\beta+w_{it}'b_i,1)$ . Marginalizing $Z_i$ over $b_i$ yields $Z_i\sim\mathcal{N}_{n_i}(X_i\beta, I_{n_i}+W_iDW_i')$ — a special case of the Chib-Greenberg (1998) multivariate probit. Algorithm 7 marginalizes over both $b_i$ and $Z$ by computing $f(y_i|\beta,D)$ via the GHK (Geweke–Hajivassiliou–Keane) method (Cholesky + sequential truncated-Normal draws + product of Normal CDFs); lag-10 β autocorrelation drops from 0.934 to 0.009 (PSID (Panel Study of Income Dynamics), two random effects (REs)) but at high computational cost. Algorithm 5 (partial blocking only) is the recommended practical choice.

Universal rule: fixed effects should always be sampled marginalized over random effects. See Gibbs Sampler and Chib-Carlin (1999).

Conjugate Likelihood Distributions and Hyperparameter Gibbs Sampling (George-Makov-Smith 1993)

In a conjugate exponential-family hierarchical model — data layer $x_i|\theta_i\sim f(x|\theta_i)$ and prior layer $\theta_i|\alpha,\beta\sim\pi(\theta|\alpha,\beta)$ — the Gibbs conditionals for $\theta_i$ are standard conjugate draws. The obstacle is the hyperparameter block: its full conditional is proportional to the conjugate likelihood distribution,

$\pi(\alpha,\beta|\theta_1,\ldots,\theta_p) \propto L(\alpha,\beta|\theta_1,\ldots,\theta_p)\cdot\pi(\alpha,\beta)$

where $L(x_0,n_0|\theta_{1:p}) \propto \exp\!\left[x_0\sum_i\theta_i - n_0\sum_i\psi(\theta_i) - p\varphi(x_0,n_0)\right]$ .

Theorem 3.3 (George-Makov-Smith 1993): $L$ is log-concave in $(\alpha,\beta)$ — enabling Gilks-Wild (1992) adaptive rejection sampling (ARS) — and jointly proper iff $p\geq 2$ . This resolves the hyperparameter bottleneck: the shape hyperparameter samples by ARS; the scale hyperparameter often draws conjugately; every step is exact.

Gamma-Poisson pump failure application ( $p=10$ , §4): With $x_i|\theta_i\sim\text{Poisson}(\theta_i t_i)$ and $\theta_i|\alpha,\beta\sim\text{Gamma}(\alpha,\beta)$ , the explicit conjugate likelihood is $L^G(\alpha,\beta|\theta_{1:p})\propto(\beta^\alpha/\Gamma(\alpha))^p(\prod_i\theta_i)^{\alpha-1}\exp(-\beta\sum_i\theta_i)$ . The complete Gibbs cycle: (a) $\theta_i|\cdot\sim\text{Gamma}(\alpha+x_i,\beta+t_i)$ ; (b) $\alpha|\cdot$ by ARS; (c) $\beta|\cdot\sim\text{Gamma}$ . 100,000 draws across four hyperprior specifications confirm shrinkage of pump-specific rates toward a common value — shrinkage calibrated by data through $\alpha$ and $\beta$ .

Why It Matters

Credibility theory (actuarial): the hierarchical posterior exactly formalizes Bühlmann's classical credibility formula — prior mean weighted by credibility factor $Z = n/(n + k)$ where $k$ = within/between variance ratio.
Small-area estimation: the Tsutakawa gamma-Poisson model shows how hierarchical partial pooling prevents both over-reliance on sparse county data and the smoothing-away of true geographic signal; a direct precursor to the mortality-modelling literature (Lee-Carter, Vector Error Correction Model [VECM] mortality).
Shrinkage in econometrics: the Minnesota Prior in Bayesian Vector Autoregressions (VARs) is conceptually analogous — each equation's coefficients are shrunk toward a common prior mean (zero), with the tightness parameter playing the role of $k$ .
Data-sparse subgroups: groups with few observations borrow heavily from the global distribution, avoiding the overfitting that OLS would produce on each group separately.
Software: BUGS makes hierarchical models almost trivially programmable; the model declaration maps directly to the graphical model structure.

Open Questions

Choice of hyperprior: improper priors at the hyperparameter level sometimes produce proper posteriors but make convergence analysis difficult.
Choice of model (HNM vs. HFLPM vs. VSHNM): model selection via marginal data density or posterior predictive loss.
Extension to crossed or nested hierarchies (multiple grouping factors).