Gelman-Pardoe (2004) Bayesian Measures of Explained Variance and Pooling in Multilevel (Hierarchical) Models

Summary

Gelman and Pardoe (2004) propose two scalar summary measures for evaluating and communicating multilevel model fit at each level of the hierarchy. The first, R² (eq. 7), measures the proportion of level-specific variance explained by the linear predictor μ; the second, λ (eq. 8), measures the fraction of residual variance that is reduced by hierarchical pooling relative to independent estimation. Both quantities are computed as free byproducts of Markov chain Monte Carlo (MCMC) posterior simulation — no null model or refitting is required. The paper also derives a Bayesian adjusted R² for linear regression with a $(n-3)/(n-p-2)$ correction factor (strictly lower than the classical Wherry 1931 formula using $(n-1)/(n-p)$ ). Applied to the Minnesota radon data (919 houses, 85 counties, 4 nested models), the measures reveal that the county-level uranium predictor explains 83% of cross-county slope variation, while hierarchical pooling nearly completely determines the slope estimates (λ_β = 0.97).

Key Claims

Level-m R² (eq. 7): $R^2 = 1 - E[V_k(\varepsilon_k)] / E[V_k(\theta_k)]$ , where $V_k$ is variance across level-m units $k$ , $\theta_k$ is the quantity of interest, and $\varepsilon_k = \theta_k - \mu_k$ is the residual after the predictor. Measures the fraction of level-m variance explained by $\mu$ ; reduces to classical R² in ordinary least squares (OLS) (no pooling), consistent with Xu (2003).
Pooling factor λ (eq. 8): $\lambda = 1 - V_k(E(\varepsilon_k)) / E[V_k(\varepsilon_k)]$ . Measures the fraction of residual variance that is reduced by the hierarchical model for the residuals $\varepsilon_k$ . λ ≈ 0 means each group is estimated independently; λ ≈ 1 means all residual structure is captured by the group-level model.
Parallel structure: R² and λ share a common denominator structure (denominator of R² = $E[V_k(\theta_k)]$ ; denominator of λ = $E[V_k(\varepsilon_k)]$ ); R² measures predictor fit to the raw quantity, λ measures hierarchical shrinkage of the residuals.
Bayesian adjusted R²: $R^2_{\text{adj}} = 1 - \frac{n-3}{n-p-2} \cdot \frac{y^T(I-H)y}{y^T I_c y}$ , strictly lower than the Wherry (1931) classical adjusted R² (factor $(n-1)/(n-p)$ ). The Bayesian formula arises from putting a uniform prior on the $p$ -dimensional parameter and integrating out the intercept separately.
Bühlmann credibility: For the balanced one-way hierarchical model with $J$ groups of size $n_j = n$ , λ specializes to $\omega = 1 - \hat\sigma^2_\alpha / (\hat\sigma^2_\alpha + \hat\sigma^2_y/n)$ , the classical credibility weight (Efron-Morris 1975; Morris 1983).
MCMC computation: Both R² and λ are computed at each MCMC draw and averaged over the posterior; no null model or separate fitting is required; the quantities come free from the standard posterior simulation.
Minnesota radon example (J=85 counties, n=919 houses; 4 nested models). Model 4 (uranium predictor at county level):

Level	$R^2$	$\lambda$	Interpretation
Household ( $y$ )	0.21	0.03	Predictor log(radon)₀ explains 21% of house-level variation; very little pooling at household level
County intercept ( $\alpha$ )	0.53	0.81	Uranium explains 53% of county-intercept variation; strong pooling
County slope ( $\beta$ )	0.83	0.97	Uranium explains 83% of county-slope variation; slopes nearly entirely from county-level model

Negative values: Both R² and λ can be negative, signaling misspecification (anti-informative predictors or variance-expanding hierarchical terms).
Cross-level interaction caveat: When individual-level predictors interact with group-level predictors, the parallel structure of R² and λ is disrupted and interpretation requires care.
Multiple imputation connection: λ is related to Rubin's (1987) fraction of missing information $\gamma$ — both measure how much information is "supplied" by the hierarchical structure relative to total residual variation (Hodges 1998; Louis 1984).

Concepts Introduced or Extended

Bayesian Hierarchical Model

Entities Mentioned

Quotes

"We propose two scalar summaries, $R^2$ and $\lambda$ ... The first, $R^2$ , is intended to measure the proportion of variance explained by the predictors in the model; the second, $\lambda$ , is intended to measure the degree of pooling from the hierarchical model for the unexplained variation." (§1)

"As with $R^2$ , the quantity $\lambda$ can be negative, and a negative value of $\lambda$ would be a sign of model misspecification." (§3)

"For a balanced one-way hierarchical model... $\lambda$ reduces to $\omega = 1 - \sigma^2_\alpha/(\sigma^2_\alpha + \sigma^2_y/n)$ , the classical pooling factor." (§3)

My Take

The two measures are elegant and practically valuable: they give multilevel modelers a simple language for communicating fit at each level ("uranium explains 83% of slope variation"), analogous to R² in OLS. The MCMC-free-byproduct property is especially useful — unlike classical model comparison, neither R² nor λ requires refitting. The Bayesian adjusted R² correction $(n-3)/(n-p-2)$ is tighter than Wherry's but the difference is small in large samples. The main limitation is the independence-of-errors assumption embedded in $V_k(\cdot)$ : in models with correlated residuals (spatial, temporal hierarchies), the measures can misrepresent explained variance. The explicit connection to Bühlmann credibility (ω = λ for balanced one-way models) helpfully bridges the actuarial and Bayesian multilevel literatures. This paper predates the ubiquity of brms/Stan; its ideas are partly embedded in posterior R² measures in those toolchains.