Lewis-Raftery (1997) Estimating Bayes Factors via Posterior Simulation With the Laplace-Metropolis Estimator

Summary

Lewis and Raftery (1997) describe the Laplace-Metropolis (LM) estimator, a way to compute integrated (marginal) likelihoods — and hence Bayes factors — from existing Markov chain Monte Carlo (MCMC) output without additional simulation runs. The method applies the classical Laplace approximation to the marginal likelihood integral, estimating the required posterior mode and posterior variance matrix directly from the MCMC sample using robust statistics. For hierarchical models with many random-effect nuisance parameters, a compound variant integrates out the random effects analytically (inner Laplace) before applying the basic LM estimator to the remaining parameters (outer Laplace). Validated against a 50-million-draw Monte Carlo gold standard on an Iran Fertility Survey example, the method recovers the correct log-integrated likelihood to within one unit using only 5,500 MCMC iterations.

Key Claims

Basic Laplace-Metropolis Estimator (fixed-effects models)

The Laplace approximation to $f(Y) = \int f(Y|\theta)f(\theta)\,d\theta$ at the posterior mode $\theta^*$ gives (eq. 4):

$$\log\hat{f}(Y) = \log f(Y|\theta^*) + \log f(\theta^*) + \frac{P}{2}\log(2\pi) + \frac{1}{2}\log|H^*|$$

where $P$ is the parameter dimension, $H^*$ is the posterior covariance matrix, and $\theta^*$ is the posterior mode. Both quantities are estimated from the MCMC output already in hand:

• Mode $\theta^*$: estimated as the multivariate median ($L_1$ center, $\arg\min_a \sum_l \|\theta^{(l)} - a\|$) rather than the posterior mean or direct likelihood maximization. Preferred because MCMC chains take occasional distant excursions, and the median is more robust to these outliers than the mean, and is generally closer to the mode for asymmetric distributions.
• Variance $H^*$: estimated using the minimum volume ellipsoid (Rousseeuw-van Zomeren 1990) weighted variance matrix — again robust to outliers.

No additional MCMC runs are required beyond the original posterior simulation. This distinguishes LM from Chib (1995), which requires $B-1$ additional reduced Gibbs runs, one per parameter block.

Compound Laplace-Metropolis Estimator (hierarchical/random-effects models)

For models with many nuisance random effects $a = (a_1,\ldots,a_N)$, direct application of the basic LM is infeasible (too many parameters). The compound LM marginalizes over the random effects analytically using a nested Laplace approximation (Section 4).

Partition parameters into: $\eta$ (fixed effects), $\xi$ (variance hyperparameter), $a$ (random effects). The compound integrated likelihood (eq. 6–7) is:

$$\log\hat{f}(Y) = \log\hat{f}^c(Y|\theta^*,\xi^*) + \log f(\theta^*,\xi^*) + \frac{P_1}{2}\log(2\pi) + \frac{1}{2}\log|H_1^*|$$

where $\hat{f}^c(Y|\theta,\xi) = \prod_i \hat{f}(Y_i|\theta,\xi)$ and each factor $\hat{f}(Y_i|\theta,\xi)$ is a Laplace approximation over $a_i$ at its conditional mode $a_i^*$.

For the logistic hierarchical model (eq. 8–11): $\text{logit}(\pi_{it}) = X_{it}\eta + a_i$, $a_i \sim \mathcal{N}(0,\xi)$:

• Conditional mode $a_i^*$ found by Newton's method
• Hessian factor: $|H_{2i}^*|^{1/2} = \prod_t \left(\frac{e^{X_{it}\eta+a_i^*}}{(1+e^{X_{it}\eta+a_i^*})^2}\right)^{-1/2}$
• The compound log-conditional likelihood for subject $i$ is closed-form (eq. in Section 5)

Variance Assessment via Batch Means

Divide the MCMC output of length $G$ into $B$ batches of $G/B$ iterations each. Compute a compound LM estimate $\widehat{LM}_b$ within each batch. These $B$ estimates are approximately independent; their mean and variance give the overall estimate and its uncertainty (Section 7).

Empirical Validation — Iran Fertility Survey

Sample of 29 women from the Iran Fertility Survey (IFS) dataset (part of the World Fertility Survey); response = birth in each exposure-year; mixed logistic model with 4 fixed effects (centered age, first-year indicator, parity, woman's education) + random effects per woman + variance hyperparameter.

Model	log-integrated likelihood	Method
4 fixed effects	−220.5	Compound LM (5,500 MCMC iters)
4 fixed effects	−221.8	Gold standard (50M prior draws)
4 fixed effects + husband's education	−223.5	Compound LM
4 fixed effects, no random effects	−222.15	Analytic Laplace (Raftery 1996b)

• BF for dropping husband's education ≈ $e^3 \approx 20$ — strong evidence for smaller model
• BF for including random effects ≈ $e^{1.65} \approx 1.35$ — weak evidence for random effects
• 95% batch-means CI for compound LM: $(−222.0, −219.0)$ — includes the gold-standard value

Position Relative to Other Methods

Method	Requires extra runs?	Works with Metropolis-Hastings (M-H)?	Random effects?	Notes
Chib (1995)	Yes ($B-1$ reduced runs)	No (Gibbs only)	Via compound version	Exact; requires tractable conditionals
Chib-Jeliazkov (2001)	Yes (short reduced run)	Yes	Via compound version	Extends Chib to any MCMC
Lewis-Raftery (1997) LM	No	Yes	Yes (compound LM)	Unimodal posterior required
Gelfand-Dey (1994)	No	Yes	Harder	Single pass; can have infinite variance
Harmonic mean	No	Yes	Harder	Inconsistent variance; unreliable
Prior MC (eq. 2)	No	—	—	Exact but requires ~50M draws

Lewis-Raftery (1997) occupy a useful middle ground: no additional runs (unlike Chib/CJ), more stable than harmonic mean, and extendable to hierarchical models. The main limitation is the Laplace approximation's reliance on approximate posterior normality (unimodal, not heavy-tailed).

Concepts Introduced or Extended

• Marginal Data Density

Entities Mentioned

• Adrian E. Raftery
• Raftery-Lewis (1992)

Quotes

"The key quantity needed for Bayesian hypothesis testing and model selection is the integrated, or marginal, likelihood of a model."

"Using this approximation, a 95% confidence interval for the compound Laplace–Metropolis estimate is (–222.0, –219.0). The actual log-integrated likelihood was –221.8. So the compound Laplace–Metropolis estimator worked well in this example."

My Take

LM is the practical, low-overhead option when (a) an MCMC run has already been done, (b) the posterior is approximately unimodal, and (c) additional reduced runs (as in Chib/CJ) would be expensive. The compound variant for hierarchical models is genuinely useful and underappreciated — most practitioners revert to the Bayesian information criterion (BIC) or harmonic mean in this setting. The batch-means uncertainty assessment is a well-motivated complement. The limitation is real: for mixtures, separated posteriors, or highly non-normal posteriors the Laplace approximation can fail silently, making Chib (1995) or Chib-Jeliazkov (2001) preferable when tractable. The "gold standard" comparison (50M-draw prior MC) is a convincing validation, making this a trustworthy reference for approximate marginal likelihood computation.