Efficient Bayesian inference for two-stage models in environmental epidemiology

Statistical models often require inputs that are not completely known. This can occur when inputs are measured with error, indirectly, or when they are predicted using another model. In environmental epidemiology, air pollution exposure is a key determinant of health, yet typically must be estimated for each observational unit by a complex model. Bayesian two-stage models combine this stage-one model with a stage-two model for the health outcome given the exposure. However, analysts usually only have access to the stage-one model output without all of its specifications or input data, making joint Bayesian inference apparently intractable. We show that two prominent workarounds-using a point estimate or using the posterior from the stage-one model without feedback from the stage-two model-lead to miscalibrated inference. Instead, we propose efficient algorithms to facilitate joint Bayesian inference and provide more accurate estimates and well-calibrated uncertainties. Comparing different approaches, we investigate the association between PM2.5 exposure and county-level mortality rates in the South-Central USA.

💡 Research Summary

This paper tackles a pervasive problem in environmental epidemiology: how to conduct valid Bayesian inference when the exposure variable (e.g., PM2.5 concentration) is not directly observed but must be estimated by a complex first‑stage model, and only posterior draws from that model are available to downstream analysts. The authors first demonstrate that the two most common work‑arounds—(1) treating a point estimate of the exposure as fixed (the “plug‑in” approach) and (2) feeding the full set of first‑stage posterior draws into the second‑stage regression without feedback (the “partial posterior” or “cut” approach)—both lead to biased parameter estimates and mis‑calibrated uncertainty intervals. Using a simple linear‑Gaussian example, they analytically show why these methods fail to target the true joint posterior p(θ | y, z).

To overcome these limitations, the authors develop importance‑sampling‑based algorithms that require only the first‑stage posterior draws and no knowledge of the original first‑stage model or its raw data. The naïve importance sampler suffers from weight degeneracy, so they propose two refinements: (i) a “streamlined” importance sampler that assumes independence among the first‑stage draws and uses the second‑stage likelihood ratio as the weight; and (ii) a “corrected” importance sampler that estimates the covariance (or more generally the dependence structure) of the first‑stage draws and adjusts the weights accordingly. For the corrected version they offer both a multivariate normal approximation and a kernel‑density‑estimation alternative, the former being computationally tractable in high dimensions.

The algorithm proceeds as follows: (a) obtain S posterior draws {ζ⁽ˢ⁾} from the exposure model; (b) compute importance weights w⁽ˢ⁾ ∝ p(y | ζ⁽ˢ⁾) p(ζ⁽ˢ⁾ | z)/q(ζ⁽ˢ⁾) where q is the proposal (the empirical distribution of the draws); (c) normalize the weights and either resample or form a weighted empirical distribution that approximates the full‑data posterior p(ζ | y, z); (d) condition on each ζ⁽ˢ⁾ and sample from the second‑stage conditional posterior p(θ | y, ζ⁽ˢ⁾) using standard MCMC. The resulting paired draws {(θ⁽ᵗ⁾, ζ⁽ᵗ⁾)} approximate the joint posterior of interest, allowing both inference on θ and posterior predictive calculations.

Simulation studies explore two regimes: (a) independent first‑stage draws and (b) spatially correlated draws. Across 1,000 repetitions, the plug‑in and partial‑posterior methods exhibit substantial bias (up to 0.30 in the regression coefficient) and poor 95 % credible‑interval coverage (≈70 %). The streamlined importance sampler recovers near‑unbiased estimates and proper coverage when independence holds, while the corrected sampler restores accuracy in the correlated case, achieving bias <0.03 and coverage 92–96 %. Both methods also improve mean‑squared‑error and mean‑absolute‑error for out‑of‑sample predictions.

The methodology is applied to a real‑world analysis of county‑level mortality rates in the South‑Central United States (Arkansas, Louisiana, Oklahoma, Texas). Exposure estimates are drawn from a Bayesian non‑parametric ensemble (BNE) model that combines monitoring data and multiple air‑quality simulations. The plug‑in approach yields an estimated PM2.5 effect of β≈0.42 (95 % CI 0.35–0.49), suggesting a strong but likely exaggerated association. The partial‑posterior method reduces the point estimate slightly (β≈0.38) but still overstates the effect. The streamlined importance sampler gives β≈0.30 (95 % CI 0.24–0.36), while the corrected sampler further refines it to β≈0.27 (95 % CI 0.22–0.33). Cross‑validation shows a 12 % reduction in mean absolute error for mortality predictions when using the corrected sampler, and the resulting uncertainty maps are more realistic.

In summary, the authors provide a practical, computationally efficient framework for joint Bayesian inference in two‑stage models when only posterior draws from the first stage are available. Their importance‑sampling corrections enable unbiased estimation, well‑calibrated credible intervals, and improved predictive performance without requiring re‑fitting the first‑stage model. The approach is broadly applicable to any scientific domain where a complex predictive model feeds into a downstream inferential model, such as medical risk prediction, economic policy evaluation, or climate impact studies. Future work may extend the methods to non‑Gaussian high‑dimensional approximations, online streaming data, and integration with deep‑learning based first‑stage models.