Reconstructing past climate variability from lake sediments can be hindered by gaps in multiproxy records. We present a modular Bayesian two-stage approach to impute missing values in a target proxy by leveraging information from all coeval proxies while explicitly propagating uncertainty. As a case study, we reconstruct \(\delta^{18}\) O and \(\delta^{13}\) C over the last 624 ka (kiloannum, meaning thousand years) of the Lake Ohrid DEEP sequence, where endogenic-calcite scarcity primarily during glacials leaves substantial isotope gaps. The workflow comprises (i) a latent Gaussian state-space model that maps irregularly sampled proxies onto a common, evenly spaced depth grid (0.25 m in this study) and accommodates reported measurement uncertainties, with optional propagation of any reported uncertainty; and (ii) a multivariate Bayesian linear regression with stochastic search variable selection (SSVS), based on the assumption that each proxy serves as an indirect indicator of multiple environmental processes that identifies a sparse, predictive subset of proxies to reconstruct the isotope series. The reconstructed stable isotope reconstructions exhibit coherent glacial–interglacial variability and gap-bridging skill across data-free intervals, with agreement in adjacent observed segments. By integrating diverse datasets within a coherent probabilistic framework and providing credible intervals for imputed values, the method improves the temporal continuity of paleoenvironmental datasets and offers a promising framework for future studies.
