Significance of uncertain phasing between the onsets of stadial–interstadial transitions in different Greenland ice core proxies

Different paleoclimate proxy records evidence repeated abrupt climate transitions during previous glacial intervals. These transitions are thought to comprise abrupt warming and increase in local precipitation over Greenland, sudden reorganization of the Northern Hemisphere atmospheric circulation, and retreat of sea ice in the North Atlantic. The physical mechanism underlying these so-called Dansgaard–Oeschger (DO) events remains debated. A recent analysis of Greenland ice core proxy records found that transitions in Na⁺ concentrations and δ¹⁸O values are delayed by about 1 decade with respect to corresponding transitions in Ca²⁺ concentrations and in the annual layer thickness during DO events. These delays are interpreted as a temporal lag of sea-ice retreat and Greenland warming with respect to a synoptic- and hemispheric-scale atmospheric reorganization at the onset of DO events and may thereby help constrain possible triggering mechanisms for the DO events. However, the explanatory power of these results is limited by the uncertainty of the transition onset detection in noisy proxy records. Here, we extend previous work by testing the significance of the reported lags with respect to the null hypothesis that the proposed transition order is in fact not systematically favored. If the detection uncertainties are averaged out, the temporal delays in the δ¹⁸O and Na⁺ transitions with respect to their counterparts in Ca²⁺ and the annual layer thickness are indeed pairwise statistically significant. In contrast, under rigorous propagation of uncertainty, three statistical tests cannot provide evidence against the null hypothesis. We thus confirm the previously reported tendency of delayed transitions in the δ¹⁸O and Na⁺ concentration records. Yet, given the uncertainties in the determination of the transition onsets, it cannot be decided whether these tendencies are truly the imprint of a prescribed transition order or whether they are due to chance. The analyzed set of DO transitions can therefore not serve as evidence for systematic lead–lag relationships between the transitions in the different proxies, which in turn limits the power of the observed tendencies to constrain possible physical causes of the DO events.

1 Introduction

In view of anthropogenic global warming, concerns have been raised that several subsystems of the earth's climate system may undergo abrupt and fundamental state transitions if temperatures exceed corresponding critical thresholds (Lenton and Schellnhuber, 2007; Lenton et al., 2008, 2019). Under sustained warming, the Atlantic Meridional Overturning Circulation (AMOC), the Amazon rainforest, or the Greenland ice sheet are, among others, possible candidates to abruptly transition to new equilibrium states that may differ strongly from their current states (Lenton et al., 2008). Understanding the physical mechanisms behind abrupt shifts in climatic subsystems is crucial for assessing the associated risks and for defining safe operating spaces in terms of cumulative greenhouse gas emissions. To date, empirical evidence for abrupt climate transitions only comes from paleoclimate proxy records encoding climate variability in the long-term past. First discovered in the δ¹⁸O records from Greenland ice cores, the so-called Dansgaard–Oeschger (DO) events are considered the archetype of past abrupt climate changes (see Fig. 1) (Johnsen et al., 1992; Dansgaard et al., 1993; Bond et al., 1993; North Greenland Ice Core Project members, 2004). These events constitute a series of abrupt regional warming transitions that punctuated the last and previous glacial intervals at millennial recurrence periods. Amplitudes of these decadal-scale temperature increases reach from 5 to 16.5 ^∘C over Greenland (Kindler et al., 2014; Huber et al., 2006; Landais et al., 2005). The abrupt warming is followed by gradual cooling over centuries to millennia before the climate abruptly transitions back to cold conditions. The relatively cold (warm) intervals within the glacial episodes have been termed Greenland stadials (GSs) (Greenland interstadials (GIs)). GSs typically show millennial-scale persistence before another abrupt warming starts a new cycle (Rasmussen et al., 2014; Ditlevsen et al., 2007). Despite being less pronounced, a global impact of DO events on climate and ecosystems is evident in many proxy records (e.g. Moseley et al., 2020; Buizert et al., 2015; Lynch-Stieglitz, 2017; Kim et al., 2012; Fleitmann et al., 2009; Voelker, 2002; Cheng et al., 2013).

Figure 1 Time series of δ¹⁸O (blue), annual layer thickness λ (cyan), Ca²⁺ (orange), and Na⁺ (green) from the NGRIP ice core, together with time series of Ca²⁺ (red) and Na⁺ (light green) from the NEEM ice core on the GICC05 timescale in ky b2k, at 10-year resolution. Light blue vertical lines mark the timings of major DO events. All time series are retrieved from Erhardt et al. (2019), and for the DO event timings and Greenland interstadial (GI) notation we followed Rasmussen et al. (2014). Originally, the δ¹⁸O record was published by NGRIP project members (2004) as 50-year mean values and later at higher resolution (5 cm) as a Supplement to Gkinis et al. (2014). The GICC05 age scale for the NGRIP ice core was compiled by Vinther et al. (2006), Rasmussen et al. (2006), Andersen et al. (2006), and Svensson et al. (2008). For the NEEM ice core, the GICC05 presented by Rasmussen et al. (2013) is used here.

Apart from δ¹⁸O, other Greenland ice core proxy records, such as Ca²⁺ and Na⁺ concentrations as well as the annual layer thickness λ, also bear the signature of DO cycles, as can be seen in Fig. 1 (e.g., Erhardt et al., 2019; Fuhrer et al., 1999; Ruth et al., 2007). While δ¹⁸O is interpreted as a qualitative proxy for ice core site temperatures (e.g. Gkinis et al., 2014; Jouzel et al., 1997; Johnsen et al., 2001), changes in Ca²⁺ concentrations – or equivalently dust – are believed to reflect changes in the atmospheric circulation (Ruth et al., 2007; Erhardt et al., 2019). Na⁺ concentration records indicate past sea-salt aerosol concentrations and are thought to negatively correlate with North Atlantic sea-ice cover (Erhardt et al., 2019; Schüpbach et al., 2018). The annual layer thickness depends on past accumulation rates at the drilling site and hence indicates local precipitation driven by synoptic circulation patterns (Erhardt et al., 2019). According to this proxy record interpretation, DO events comprise not only sudden warming but also a sudden increase in local precipitation amounts, retreat of the North Atlantic sea-ice cover, and changes in hemispheric circulation patterns.

In the search for the mechanism(s) causing or triggering DO events, several attempts have been made to deduce the relative temporal order of these abrupt changes by analyzing the phasing of corresponding abrupt shifts detected in multi-proxy time series from Greenland ice cores (Erhardt et al., 2019; Thomas et al., 2009; Steffensen et al., 2008; Ruth et al., 2007). While Thomas et al. (2009) and Steffensen et al. (2008) report delayed Greenland warming with respect to atmospheric changes for the onsets of GI-8 and GI-1 and the Holocene, Ruth et al. (2007) find no systematic lead or lag across the onsets of GI-1 to GI-24. However, the comprehensive study conducted by Erhardt et al. (2019) concludes that, on average, initial changes in both terrestrial dust aerosol concentrations (Ca²⁺) and local precipitation (λ) have preceded the changes in local temperatures (δ¹⁸O) and sea-salt aerosol concentrations (Na⁺) by roughly 1 decade at the onset of DO events during the last glacial cycle.

These observation-based studies are complemented by numerous conceptual theories and modeling studies that explore a variety of mechanisms to explain the DO events. Many studies emphasize the role of the AMOC in the emergence of DO events (Broecker et al., 1985; Clark et al., 2002; Ganopolski and Rahmstorf, 2001; Henry et al., 2016). In this context, Vettoretti and Peltier (2018) identified a self-sustained sea-salt oscillation mechanism to initiate transitions between stadials and interstadials in a comprehensive general circulation model (GCM) run, while Boers et al. (2018) proposed a coupling between sea-ice growth, subsurface warming, and AMOC changes to explain the DO cycles. Moreover, Li and Born (2019) draw attention to the subpolar gyre, a sensitive region that features strong interactions between atmosphere, ocean, and sea ice. In line with the empirical studies that suggest a delayed Greenland warming with respect to atmospheric changes, Kleppin et al. (2015) and Zhang et al. (2014) find DO-like transitions in GCM studies triggered by an abrupt reorganization of atmospheric circulation patterns.

Here, we refine the investigation of a potential pairwise lead–lag relationship between the four climate proxies Ca²⁺, Na⁺, δ¹⁸O, and the annual layer thickness λ at DO transition onsets, as previously presented by Erhardt et al. (2019), by rigorously taking into account the uncertainties of the DO onset detection in the different proxy records. We use the same data and the same probabilistic transition onset detection method as provided by Erhardt et al. (2019). The data comprise piecewise high-resolution (7 years or higher) multi-proxy time series around 23 major DO events from the later half of the last glacial cycle, from the NEEM and the NGRIP ice cores (Erhardt et al., 2019). The fact that high-frequency internal climate variability blurs abrupt transitions limits the ability to precisely detect their onset in the proxy data and thereby constitutes the main obstacle for the statistical analysis of the succession of events. The method designed by Erhardt et al. (2019) very conveniently takes this into account and instead of returning scalar estimators it quantifies the transition onsets in terms of Bayesian posterior probability densities that indicate the plausibility of a transition onset at a certain time in view of the data. This gives rise to a set of uncertain DO transition onset lags for each pair of proxies under study, whose statistical interpretation is the goal of this study.

While Erhardt et al. (2019) report transition onsets, midpoints, and endpoints, we restrict our investigation to the transition onset points, since we consider the leads and lags between the initial changes in the different proxy records to be the relevant quantity for a potential identification of the physical trigger of the DO events. We extend the previous work by interpreting the sets of uncertain lags as samples generated in random experiments from corresponding unknown populations – each proxy pair is associated with its own population of lags. This allows for the investigation of whether the reported average lags (Erhardt et al., 2019) are a systematic feature or whether they might have emerged by chance. In order to review the statistical evidence for potential systematic lags, we formalize the notion of a “systematic lag”: we call a lag systematic if it is enshrined in the random experiment in form of a population mean different from 0. Samples generated from such a population with a non-zero mean would systematically (and not by chance) exhibit sample means different from 0. Accordingly, we formulate the null hypothesis that the proposed transition sequence is in fact not physically favored. In mathematical terms this corresponds to an underlying population of lags with a mean equal to 0 or with reversed sign with respect to the observed lags. A rejection of this null hypothesis would statistically corroborate the interpretation that transitions in δ¹⁸O and Na⁺ systematically lag their counterparts in λ and Ca²⁺. On the other hand, acceptance of the hypothesis would prevent us from ruling out that the observed lag tendencies are a coincidence and not a systematic feature. We have identified three different statistical tests suitable for this task, which all rely on slightly different assumptions. Therefore, in combination they yield a robust assessment of the observations. Most importantly, we propagate the uncertainties that arise from the transition onset detection to the level of p values of the different tests.

We will show that, if the uncertainties are averaged out at the level of the individual transition onset lags – thus ignoring the uncertainties in the onset detection – all tests indicate statistical significance (at 5 % confidence level) of the observed tendencies toward delayed δ¹⁸O and Na⁺ transition onsets with respect to the corresponding onsets in λ and Ca²⁺. Rigorous uncertainty propagation, however, yields substantial probabilities for the observed transition onset lags to be non-significant with respect to the null hypothesis. We thus argue that the uncertainties in the transition onset detection are too large to infer a population mean different from 0 in the direction of the observed tendencies. In turn, this prevents the attribution of the observed lead–lag relations to a fundamental mechanism underlying the DO events. We discuss the difference between our approach and the one followed by Erhardt et al. (2019) in detail below.

In addition to the quantitative uncertainty discussed here, there is always qualitative uncertainty about the interpretation of climate proxies. Clearly, there is no one-to-one mapping between proxy variables and the climate variables they are assumed to represent. To give an example, changes in the atmospheric circulation will simultaneously impact the transport efficiency of sea-salt aerosols to Greenland. Schüpbach et al. (2018) discuss in detail the entanglement of transport efficiency changes and source emission changes for aerosol proxies measured in Greenland ice cores. We restrict our analysis to those proxy pairs that have been found to show decadal-scale time lags by Erhardt et al. (2019) and leave aside those pairs which show almost simultaneous transition onsets according to Erhardt et al. (2019).

This article is structured as follows: first, the data used for the study are described. Second, we introduce our methodology in general terms, in order to facilitate potential adaptation to structurally similar problems. Within this section, we pay special attention to clarifying the differences between the approaches chosen in this study and by Erhardt et al. (2019). This is followed by the presentation of our results including a comparison to previous results. In the subsequent discussion, we give a statistical interpretation and explain how the two lines of inference lead to different conclusions. The last section summarizes the key conclusions that can be drawn from our analysis.

2 Data

In conjunction with their study, Erhardt et al. (2019) published 23 highly resolved time series for Ca²⁺ and Na⁺ concentrations from the NGRIP and NEEM ice cores for time intervals of 250 to 500 years centered around DO events from the later half of the last glacial. The data set covers all major interstadial onsets from GI-17.2 to the Holocene, as determined by Rasmussen et al. (2014). The time resolution decreases from 2 to 4 years with increasing depth in the ice cores due to the thinning of the core. In addition, Erhardt et al. (2019) derived the annual layer thickness from the NGRIP aerosol data and published these records likewise for the time intervals described above. Furthermore, continuous 10-year resolution versions of the proxy records were published, which cover the period 60–10 kyr BP, shown in Fig. 1 (Erhardt et al., 2019). Finally, the NGRIP δ¹⁸O record at 5 cm resolution (corresponding to 4–7 years for the respective time windows) (North Greenland Ice Core Project members, 2004) completes the data set used in the study by Erhardt et al. (2019) and correspondingly in our study.

While Ca²⁺ and Na⁺ mass concentrations are interpreted as indicators of the past state of the atmospheric large-scale circulation and the past North Atlantic sea-ice extent, respectively, the annual layer thickness and δ¹⁸O records give qualitative measures of the local precipitation and temperature, respectively (Erhardt et al., 2019, and references therein). The high resolution and the shared origin of the time series make them ideally suited to study the succession of events at the beginning of DO transitions. On top of that, the aerosol data have been co-registered in a continuous flow analysis allowing for the highest possible comparability (Erhardt et al., 2019).

Figure 2 DO events (Greenland interstadial onsets) for which Erhardt et al. (2019) provide high-resolution proxy data (${\mathrm{Ca}}^{\mathrm{2}+},{\mathrm{Na}}^{+}$, and λ) for windows centered around the transitions. δ¹⁸O data for the corresponding windows were retrieved from continuous δ¹⁸O time series measured in 5 cm steps in the NGRIP ice core (see Fig. 1). The posterior probability densities for the transition onsets with respect to the timing of the DO events according to Rasmussen et al. (2014) are shown in arbitrary units for all proxies. They were recalculated using the data and the method provided by Erhardt et al. (2019). The uncertain transition onsets are only shown for those transitions investigated in this study – the selection is adopted from Erhardt et al. (2019) to guarantee comparability.

For their analysis, Erhardt et al. (2019) only considered time series around DO events that do not suffer from substantial data gaps. For the sake of comparability, we adopt their selection. From Fig. 2 it can be inferred which proxy records around which DO events have been included in this study. For details on the data and the proxy interpretations we refer to Erhardt et al. (2019) and the manifold references therein.

3 Methods

We first briefly review the probabilistic method that we adopted from Erhardt et al. (2019) in order to estimate the transition onset time t₀ of each proxy variable for each DO event comprised in the data (see Fig. 3). The Bayesian method accounts for the uncertainty inherent to the determination of t₀ by returning probability densities ${\mathit{\rho }}_{T_{\mathrm{0}}}\left(t_{\mathrm{0}}\right)$ instead of scalar estimators. From these distributions, corresponding probability distributions for the pairwise time lags between two proxies can be derived for all DO events. Second, a statistical perspective on the series of DO events is established. For a given proxy pair, the set of transition onset lags from the different DO events is treated as a sample of observations from an unknown underlying population. In this very common setup, naturally one would use hypothesis tests to constrain the population. In particular, the question of whether any lag tendencies observed in the data are a systematic feature or whether they have instead occurred by chance can be assessed by testing a convenient null hypothesis. However, the particularity that the individual observations that comprise the sample are themselves subject to uncertainty requires a generalization of the hypothesis tests. We propagate the uncertainty of the transition onset timings to the p values of the tests and hence obtain uncertain p values in terms of probability densities (see Figs. 4 and 7). While in common hypothesis tests the scalar p value is compared to a predefined significance level, here we propose two criteria to project the p-value distribution onto the binary decision between acceptance and rejection of the null hypothesis. After this general characterization of the statistical problem, we introduce the tests which we employ for the analysis. Finally, we compare our approach to the one followed by Erhardt et al. (2019).

Figure 3 (a) Posterior probability distribution ${\mathit{\rho }}_{T_{\mathrm{0}}}\left(t_{\mathrm{0}}\right)$ for the onset of NGRIP Ca²⁺ and Na⁺ transitions associated with the onset of GI-12c, derived from Ca²⁺ (orange) and Na⁺ (green) values around the GI-12c onset at 2-year resolution, using the probabilistic ramp-fitting shown in (b). The black lines in (b) indicate the expected ramp, i.e., the average over all ramps determined by the posterior distributions of the ramp parameters. The grey shaded area indicates the 5th–95th percentiles of these ramps. (c) Histogram sampled from the posterior distribution for the transition onset lag Δt between the two proxies (violet), together with the corresponding Gaussian kernel density estimate (KDE, blue).

Figure 4 (a) Schematic representation of an uncertain observation of a sample (purple) generated from a population (blue) in a random experiment. The blue line indicates the probability density of the generating population 𝒫_X. Green lines indicate the true value of a sample $\mathit{x}=(x_{\mathrm{1}},\mathrm{\dots },x_{\mathrm{6}})$ realized from 𝒫_X. If the observational process involves uncertainty, a second level of randomness is introduced and the values can at best be approximated by probability density functions depicted in purple. These uncertainty distributions indicate the informed estimate of the observer about how likely a certain value ${\widehat{y}}_i$ for the estimator ${\widehat{Y}}_i$ is to coincide with the true value x_i. Depending on the measurement process, the uncertainty distributions of the sample members may all exhibit individual shapes or they may share a common one. (b) Distribution of the uncertain test statistic $\widehat{\mathrm{\Phi }}=\mathit{\varphi }\left(\widehat{\mathit{Y}}\right)$ derived from the uncertain sample (purple) together with the corresponding value derived from the true sample (green). In olive, the distribution of Φ under the null hypothesis is shown. The dotted red line separates the rejection region (left) from the acceptance region in a one-sided test setup. (c) Distribution of the uncertain p value corresponding to the uncertain sample. In green, the p value of the certain sample is marked. The red line indicates the significance level α.

3.1 Transition onset detection

Consider a fluctuating time series $\mathcal{D}=\mathit{\left\{}x\right(t_i){\mathit{\}}}_{i=\mathrm{1},\mathrm{\dots },n}$ with n data points, which includes one abrupt transition from one level of values to another, as shown in Fig. 3b. For this setting, Erhardt et al. (2019) have designed a method to estimate the transition onset time t₀ in a probabilistic, Bayesian sense. The application of the method to NGRIP Ca²⁺ and Na⁺ concentration data around the onset of GI-12c is illustrated in Fig. 3. Instead of a point estimate, their method returns a so-called posterior probability density that indicates the plausibility of the respective onset time in view of the data (see Fig. 3a). For technical reasons, this probability density cannot be derived in form of a continuous function but only in form of a representative set of values generated from it by means of a Markov chain Monte Carlo (MCMC) algorithm (Goodman and Weare, 2010).

The key idea is to model the transition as a linear ramp ℒ(t_i) perturbed by red noise ϵ(t_i), which is an autoregressive process of first order:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }y\left(t_i\right)=\underset{\text{liner ramp}\mathcal{L}\left(t_i\right)}{\underbrace{\left\{\begin{array}{ll}{y}_{\mathrm{0}}& t_i\le t_{\mathrm{0}}\\ y_{\mathrm{0}}+\mathrm{\Delta }y\frac{t_i-t_{\mathrm{0}}}{\mathit{\tau }}& t_{\mathrm{0}}<t_i<t_{\mathrm{0}}+\mathit{\tau }\\ y_{\mathrm{0}}+\mathrm{\Delta }y& t\ge t_i+\mathit{\tau }\end{array}\right.}}+\underset{\text{AR}(\mathrm{1}{)}_{\mathit{\sigma },\mathit{\alpha }}}{\underbrace{\mathit{ϵ}\left(t_i\right)}}.\end{array}$$

This model is fully determined by the four ramp parameters $\mathit{\{}{t}_{\mathrm{0}},y_{\mathrm{0}},\mathit{\tau },\mathrm{\Delta }y\mathit{\}}$, the amplitude σ, and the autoregressive coefficient α of the AR(1) process. For a given configuration θ of these six parameters, the probability of this stochastic model to exactly reproduce the data 𝒟 reads

$$\begin{array}{ll}{\displaystyle }\mathit{\pi }\left(\mathcal{D}\right|\mathit{\theta })& {\displaystyle }:=\mathit{\pi }\left(y\right(t_i)=x(t_i)\forall i\in \mathit{\{}\mathrm{1},\mathrm{\dots },n\mathit{\left\}}\right|\mathit{\theta })\\ \text{(2)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{{\left(\sqrt{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}^{\mathrm{2}}}\right)}^n}}\prod _{i=\mathrm{1}}^n\exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{({\mathit{\delta }}_i-\mathit{\alpha }{\mathit{\delta }}_{i-\mathrm{1}}{)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}\right),\end{array}$$

where ${\mathit{\delta }}_i=x\left(t_i\right)-\mathcal{L}\left(t_i\right)$ denotes the residuals between the linear ramp and the observations and δ₀=0. Bayes' theorem immediately yields the posterior probability density for the model parameters π(θ|𝒟) upon introduction of convenient priors π(θ):

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }\mathit{\pi }\left(\mathit{\theta }\right|\mathcal{D})={\displaystyle \frac{\mathit{\pi }\left(\mathcal{D}\right|\mathit{\theta }\left)\mathit{\pi }\right(\mathit{\theta })}{\mathit{\pi }\left(\mathcal{D}\right)}},\end{array}$$

where the normalization constant $\mathit{\pi }\left(\mathcal{D}\right)=\int \mathit{\pi }\left(\mathcal{D}\right|\mathit{\theta }\left)\mathit{\pi }\right(\mathit{\theta })\mathrm{d}\mathit{\theta }$ is the a priori probability of the observations. Since the parameter space is six-dimensional, Eq. (3) cannot be evaluated explicitly on a grid with reasonably fine spacing. Instead, an MCMC algorithm is used to sample a representative set $\mathit{\{}{\mathit{\theta }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\theta }}_m\mathit{\}}$ of parameter configurations from the posterior distribution that approximates the continuous distribution in the sense that for smooth functions f

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }\int f\left(\mathit{\theta }\right){\mathit{\rho }}_{\mathbf{\Theta }}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\simeq \int f\left(\mathit{\theta }\right){\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathbf{\Theta }}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }={\displaystyle \frac{\mathrm{1}}{m}}\sum _{j=\mathrm{1}}^{m}f\left({\mathit{\theta }}_j\right),\end{array}$$

where the notion of a so-called empirical distribution ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathbf{\Theta }}\left(\mathit{\theta }\right)=\frac{\mathrm{1}}{m}{\sum }_{j=\mathrm{1}}^m\mathit{\delta }(\mathit{\theta }-{\mathit{\theta }}_j)$ has been used. The use of the MCMC algorithm further allows us to omit the normalization constant π(𝒟). The number m of individuals comprised in the MCMC sample must be chosen large enough to ensure a good approximation in Eq. (4). The marginal distribution for the parameter t₀ relevant for our study can be obtained by integration over the remaining parameters θ^*:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }{\mathit{\rho }}_{T_{\mathrm{0}}|\mathcal{D}}\left(t_{\mathrm{0}}\right)=\int \mathit{\pi }\left(\mathit{\theta }\right|\mathcal{D})\mathrm{d}{\mathit{\theta }}^{*},\end{array}$$

which reads

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{T_{\mathrm{0}}}\left(t_{\mathrm{0}}\right)={\displaystyle \frac{\mathrm{1}}{m}}\sum _{j=\mathrm{1}}^m\mathit{\delta }(t_{\mathrm{0}}-t_{\mathrm{0},j})\end{array}$$

in terms of the empirical density induced by the MCMC sample.

Given the probability densities for the transition onsets of two proxy variables p and q at a chosen DO event i, the probability density for the lag $\mathrm{\Delta }{t}_i^{p,q}=t_{\mathrm{0}}^{p,i}-t_{\mathrm{0}}^{q,i}$ between them reads

$$\begin{array}{ll}{\displaystyle }{\mathit{\rho }}_{\mathrm{\Delta }{T}_i^{p,q}}\left(\mathrm{\Delta }{t}_i^{p,q}\right)=& {\displaystyle }\iint \mathit{\delta }\left(t_{\mathrm{0}}^{p,i}-t_{\mathrm{0}}^{q,i}-\mathrm{\Delta }{t}_i^{p,q}\right){\mathit{\rho }}_{T_{\mathrm{0}}}^{p,i}\left(t_{\mathrm{0}}^{p,i}\right)\\ \text{(7)}& {\displaystyle }& {\displaystyle }\times {\mathit{\rho }}_{T_{\mathrm{0}}}^{q,i}\left(t_{\mathrm{0}}^{q,i}\right)\mathrm{d}{t}_{\mathrm{0}}^{p,i}\mathrm{d}{t}_{\mathrm{0}}^{q,i}.\end{array}$$

$\mathrm{\Delta }{T}_i^{p,q}$ was chosen to denote the time lag which inherits the uncertainty from the transition onset detection and must thus mathematically be treated as a random variable. $\mathrm{\Delta }{t}_i^{p,q}$ denotes a potential value that $\mathrm{\Delta }{T}_i^{p,q}$ may assume. The set of probability densities $\mathit{\left\{}{\mathit{\rho }}_{\mathrm{\Delta }{T}_i^{p,q}}\right(\mathrm{\Delta }{t}_i^{p,q}){\mathit{\}}}_i$ derived from the different DO events conveniently describes the random vector of uncertain DO onset lag observations $\mathbf{\Delta }{\mathit{T}}^{p,q}=(\mathrm{\Delta }{T}_{\mathrm{1}}^{p,q},\mathrm{\dots },\mathrm{\Delta }{T}_n^{p,q})$ for the (p,q) proxy pair in the sense that

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }{\mathit{\rho }}_{\mathbf{\Delta }{\mathit{T}}^{p,q}}\left(\mathbf{\Delta }{\mathit{t}}^{p,q}\right)=\prod _{i=\mathrm{1}}^n{\mathit{\rho }}_{\mathrm{\Delta }{T}_i^{p,q}}\left(\mathrm{\Delta }{t}_i^{p,q}\right).\end{array}$$

Note that the entries $\mathrm{\Delta }{T}_i^{p,q}$ of the random vector ΔT^p,q are independent from each other and follow their individual distributions ${\mathit{\rho }}_{\mathrm{\Delta }{T}_i^{p,q}}\left(\mathrm{\Delta }{t}_i^{p,q}\right)$, such that the joint distribution is given by the product of the individual distributions. A cross-core comparison is not possible because the relative dating uncertainties between the cores exceed the magnitude of the potential time lags.

For the sake of simplicity, we omit the difference between the posterior density distribution and the empirical posterior density distribution induced by an MCMC sample. It is shown in Appendix A that all methods can be equivalently formulated in terms of the empirical posterior density distribution. The numerical computations themselves have of course been carried out with the empirical densities obtained from the MCMC sampler. Appendix B discusses the construction of numerically manageable empirical densities ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathbf{\Delta }{\mathit{T}}^{p,q}}\left(\mathbf{\Delta }{\mathit{t}}^{p,q}\right)$. Since a substantial reduction in the available MCMC sampled data is required, a control group of alternative realizations of ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathbf{\Delta }{\mathit{T}}^{p,q}}\left(\mathbf{\Delta }{\mathit{t}}^{p,q}\right)$ is introduced. The high agreement of the results obtained from the control group with the results discussed in the main text confirms the validity of the initial ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathbf{\Delta }{\mathit{T}}^{p,q}}\left(\mathbf{\Delta }{\mathit{t}}^{p,q}\right)$ construction.

In the following, all probability densities that represent uncertainties with origin in the transition onset observation will be referred to as uncertainty distributions or uncertainty densities. This helps to distinguish them from probability distributions that generically characterize random experiments. The random variables described by uncertainty distributions will be termed uncertain variables and will be marked with a hat. Generally, we denote all random (uncertain) variables by capital letters X ($\widehat{X}$), while realizations will be denoted with lower-case letters x ($\widehat{x}$). Furthermore, distributions will always be subscripted with the random variables that they characterize, e.g., ρ_X(x) (${\mathit{\rho }}_{\widehat{X}}\left(\widehat{x}\right)$). For the sake of readability, sometimes we omit the index p,q when it is clear that a quantity refers to a pair of proxies (p,q).

3.2 Statistical setting

Despite their diversity in terms of temperature amplitude, duration, and frequency across the last glacial, the reoccurring patterns and their common manifestation in different proxies suggest that the DO events follow a common physical mechanism. If this assumption holds true, this mechanism prescribes a fixed pattern of causes and effects for all DO events – at least on the scale of interactions between climatic subsystems represented by the proxies under study. However, natural variability randomly delays or advances the individual parts of the event chain of the DO mechanism in each single realization, without violating the mechanistic causality. The observed pairwise transition onset lags can thus be regarded as realizations of independent and identically distributed (i.i.d.) random variables generated in a random experiment $(\mathrm{\Omega },\mathcal{F},{\mathcal{P}}_{\mathrm{\Delta }T}^{p,q})$ on the sample space Ω=ℝ. Here, ℱ is a σ algebra defined on Ω and may be taken as the Borel algebra. ${\mathcal{P}}_{\mathrm{\Delta }T}^{p,q}$ – the so-called population – denotes a probability measure with respect to ℱ and fully characterizes the random lag ΔT^p,q between the proxies p and q. Importantly, if any of the proxy variables investigated here was to represent a climate variable associated with the DO event trigger, we would expect an advanced initial change in the record of this proxy with respect to other proxies at DO events. In turn, a pronounced delay of a proxy record's transition onset contradicts the assumption that the proxy represents a climate variable associated with the trigger. Therefore, the identification of leads and lags between the transition onsets in the individual proxy time series may help in the search for the trigger of the DO events. Here, we formalize the investigation of systematic lead–lag relationships between the proxy transitions. The random experiment framework allows us to relate a suspected transition sequence to a mean of the generating population ${\mathcal{P}}_{\mathrm{\Delta }T}^{p,q}$ differing from 0 in the according direction. Evidence for the suspected sequence can then be obtained by testing the null hypothesis of a population mean equal to 0 or with a sign opposed to the suspected lag direction. If this null hypothesis can be rejected based on the observations, this would constitute a strong indication of a systematic, physical lag and would hence potentially yield valuable information on the search for the mechanism(s) and trigger(s) of the DO transitions.

According to the data selection by Erhardt et al. (2019) as explained in Sect. 2, for all studied pairs of proxies we compute either 16 or 20 transition lags from the different DO events, which we interpret as samples $\mathbf{\Delta }{\mathit{t}}^{p,q}=(\mathrm{\Delta }{t}_{\mathrm{1}}^{p,q},\mathrm{\dots },\mathrm{\Delta }{t}_n^{p,q})$ from their respective populations ${\mathcal{P}}_{\mathrm{\Delta }T}^{p,q}$. Studying these samples, Erhardt et al. (2019) deduced a decadal-scale delay in the transition onsets in Na⁺ and δ¹⁸O records with respect to their counterparts in Ca²⁺ and λ. In order to test if the data support evidence for this lag to be systematic in a statistical sense, the notion of a “systematic lag” first needs to be formalized mathematically. We argue that a physical process that systematically delays one of the proxy variable transitions with respect to another must in the random experiment framework be associated with a population that exhibits a mean different from 0: ${\mathit{\mu }}^{p,q}=E\left(\mathrm{\Delta }{T}^{p,q}\right)\ne \mathrm{0}$. The outcomes of such a random experiment will systematically exhibit sample means different from 0 in accordance with the population mean. Samples generated from a population with a mean equal to 0 may as well yield sample means that strongly differ from 0. However, their occurrence is not systematic but rather a coincidence. Given a limited number of observations, hypothesis tests provide a consistent yet not unambiguous way to distinguish systematic from random features. If the mean of the observed sample $u^{p,q}\left(\mathbf{\Delta }{\mathit{t}}^{p,q}\right)=\frac{\mathrm{1}}{n}\sum \mathrm{\Delta }{t}_i^{p,q}$ indicates an apparent lag between the proxies p and q, testing whether the sample statistically contradicts a population that favors no (${\mathit{\mu }}^{p,q}=\mathrm{0}$) or the opposite lag (or $\text{sign}\left({\mathit{\mu }}^{p,q}\right)\ne \text{sign}\left(u^{p,q}\right)$) can provide evidence at the significance level α for the observed mean lag to be systematic in the sense that $\text{sign}\left({\mathit{\mu }}^{p,q}\right)=\text{sign}\left(u^{p,q}\right)$. However, as long as the null hypothesis cannot be rejected, the observed average sample lag cannot be regarded as statistical evidence for a systematic lag.

Before we introduce the tests deployed for this study, we discuss the particularity that the individual observations of the i.i.d. variables that comprise our samples are themselves subject to uncertainty and hence are represented by probability densities instead of scalar values. The common literature on hypothesis tests assumes that an observation of a random variable yields a scalar value. Given a sample of n scalar observations

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }\mathit{x}=\left(x_{\mathrm{1}},x_{\mathrm{2}},\mathrm{\dots },x_n\right),\end{array}$$

the application of hypothesis tests to the sample is in general straight forward and has been abundantly discussed (e.g. Lehmann and Romano, 2006). In short, a test statistic ϕ_x=ϕ(x) is computed from the observed sample, where

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }\mathit{\varphi }:{\mathbb{R}}^n\to \mathbb{R},\mathit{x}\mapsto \mathit{\varphi }\left(\mathit{x}\right)\end{array}$$

denotes the mapping from the space of n-dimensional samples to the space of the test statistic and ϕ_x denotes the explicit value of the function when applied to the observed sample x. Subsequently, integration of the so-called null distribution over all values ϕ^′, which under the null hypothesis H₀ are more extreme than the observed ϕ_x, yields the test's p value. In this study, a hypothesis on the lower limit of a parameter will be tested. In this one-sided left-tailed application of hypothesis testing, the p value explicitly reads

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }{p}_{\mathit{x}}=\underset{-\mathrm{\infty }}{\overset{{\mathit{\varphi }}_{\mathit{x}}}{\int }}{\mathit{\rho }}_{\mathrm{\Phi }}^{\mathrm{0}}\left({\mathit{\varphi }}^{\prime }\right)\mathrm{d}{\mathit{\varphi }}^{\prime },\end{array}$$

which defines the mapping

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }p:\mathbb{R}\to [\mathrm{0},\mathrm{1}],{\mathit{\varphi }}_{\mathit{x}}\mapsto p\left({\mathit{\varphi }}_{\mathit{x}}\right)=p_{\mathit{x}}.\end{array}$$

Analogous expressions may be given for one-sided right-tailed and two-sided tests. The null distribution ${\mathit{\rho }}_{\mathrm{\Phi }}^{\mathrm{0}}\left({\mathit{\varphi }}^{\prime }\right)$ is the theoretical distribution of the random test statistic Φ=ϕ(X) under the assumption that the null hypothesis on the population 𝒫_X holds true. If p_x is less than a predefined significance level α, the observed sample x is said to contradict the null hypothesis at the significance level α, and the null hypothesis should be rejected.

In contrast to this setting, the DO transition onset lags $\mathrm{\Delta }{t}_i^{p,q}$ between the proxies p and q, which are thought to have been generated from the population ${\mathcal{P}}_{\mathrm{\Delta }T}^{p,q}$, are observed with uncertainty. In our case, the entries in the vector of observations are uncertain variables themselves, which are characterized by the previously introduced uncertainty distributions ${\mathit{\rho }}_{\mathrm{\Delta }{\widehat{T}}_i^{p,q}}\left(\mathrm{\Delta }{\widehat{t}}_i^{p,q}\right)$. Figure 4a illustrates this situation: from an underlying population 𝒫_X a sample $\mathit{x}=(x_{\mathrm{1}},\mathrm{\dots },x_{\mathrm{6}})$ is realized, with x_i denoting the true values of the individual realizations. However, the exact value of x_i cannot be measured precisely due to measurement uncertainties. Instead, an estimator ${\widehat{Y}}_i$ is introduced together with the uncertainty distribution ${\mathit{\rho }}_{{\widehat{Y}}_i}\left({\widehat{y}}_i\right)$ that expresses the observer's belief about how likely a specific value ${\widehat{y}}_i$ for the estimator ${\widehat{Y}}_i$ is to agree with the true value x_i. The ${\widehat{Y}}_i$ correspond to the $\mathrm{\Delta }{\widehat{T}}_i^{p,q}$. For the x_i there is no direct correspondence in the problem at hand because this quantity cannot be accessed in practice and is thus not denoted explicitly. We call the vector of estimators $\widehat{\mathit{Y}}=({\widehat{Y}}_{\mathrm{1}},\mathrm{\dots },{\widehat{Y}}_n)$ an uncertain sample in the following. Omitting the (p,q) notation, we denote an uncertain sample of time lags as

$$\begin{array}{}\text{(13a)}& \mathbf{\Delta }\widehat{\mathit{T}}=\left(\mathrm{\Delta }{\widehat{T}}_{\mathrm{1}},\mathrm{\Delta }{\widehat{T}}_{\mathrm{2}},\mathrm{\dots },\mathrm{\Delta }{\widehat{T}}_n\right),\end{array}$$

with

$$\begin{array}{}\text{(13b)}& {\mathit{\rho }}_{\mathbf{\Delta }\widehat{\mathit{T}}}\left(\mathbf{\Delta }\widehat{\mathit{t}}\right)=\prod _{i=\mathrm{1}}^n{\mathit{\rho }}_{\mathrm{\Delta }{\widehat{T}}_i}\left(\mathrm{\Delta }{\widehat{t}}_i\right).\end{array}$$

Note that the uncertainty represented by the uncertain sample originates from the observation process – the sample no longer carries the generic randomness of the population 𝒫_ΔT it was generated from. The $\mathrm{\Delta }{\widehat{T}}_i$ are no longer identically but yet independently distributed.

A simplistic approach to test hypotheses on an uncertain sample would be to average over the uncertainty distribution and subsequently apply the test to the resulting expected sample

$$\begin{array}{}\text{(14)}& \begin{array}{rl}\mathrm{E}\left(\mathbf{\Delta }\widehat{\mathit{T}}\right)& =\left(\mathrm{E}\left(\mathrm{\Delta }{\widehat{T}}_{\mathrm{1}}\right),\mathrm{\dots },\mathrm{E}\left(\mathrm{\Delta }{\widehat{T}}_n\right)\right)\\ & =\left(\int \mathrm{\Delta }{\widehat{t}}_{\mathrm{1}}{\mathit{\rho }}_{\mathrm{\Delta }{\widehat{T}}_{\mathrm{1}}}\left(\mathrm{\Delta }{\widehat{t}}_{\mathrm{1}}\right)\mathrm{d}\mathrm{\Delta }{\widehat{t}}_{\mathrm{1}},\mathrm{\dots },\right.\\ & \left.\int \mathrm{\Delta }{\widehat{t}}_n{\mathit{\rho }}_{\mathrm{\Delta }{\widehat{T}}_n}\left(\mathrm{\Delta }{\widehat{t}}_n\right)\mathrm{d}\mathrm{\Delta }{\widehat{t}}_n\right).\end{array}\end{array}$$

Averaging out uncertainties, however, essentially implies that the uncertainties are ignored and is thus always associated with a loss of information. The need for a more thorough treatment, with proper propagation of the uncertainties, may be illustrated by a simple consideration. Assume that a random variable X can be observed at a given precision σ_obs, where σ_obs may be interpreted as the typical width of the corresponding uncertainty distribution. For any finite number of observations of X, the measurement uncertainty limits the ability to infer properties of the population. For example, one cannot distinguish between potential candidates μ₀ and μ₁ for the population mean, whose difference $\mathrm{\Delta }\mathit{\mu }=|{\mathit{\mu }}_{\mathrm{0}}-{\mathit{\mu }}_{\mathrm{1}}|$ is a lot smaller than the observational precision, unless the number of observations tends to infinity. If uncertainties are averaged out, testing $H_{\mathrm{0}}:\mathit{\mu }={\mathit{\mu }}_{\mathrm{0}}$ against the alternative $H_{\mathrm{1}}:\mathit{\mu }={\mathit{\mu }}_{\mathrm{1}}$ can eventually still yield significance, even in the case where $|{\mathit{\mu }}_{\mathrm{0}}-{\mathit{\mu }}_{\mathrm{1}}|\ll {\mathit{\sigma }}_{\text{obs}}$. For relatively small sample sizes, such an attested significance would be statistically meaningless. The uncertainties in the individual measurements considered here are on the order of a few decades, while the proposed size of the investigated time lag is roughly 1 decade. In combination with the relatively small sample sizes of either 16 or 20 events, the scales involved in the analysis require a suitable treatment of the measurement uncertainties.

The uncertainty propagation relies on the fact that applying a function $f:\mathbb{R}\to \mathbb{R}$ to a real valued random (uncertain) variable X yields a new random (uncertain) variable G=f(X), which is distributed according to

$$\begin{array}{}\text{(15)}& {\displaystyle }{\displaystyle }{\mathit{\rho }}_G\left(g\right)=\int \mathit{\delta }\left(f\left(x\right)-g\right){\mathit{\rho }}_X\left(x\right)\mathrm{d}x.\end{array}$$

Analogously, the uncertain test statistic $\widehat{\mathrm{\Phi }}=\mathit{\varphi }\left(\mathbf{\Delta }\widehat{\mathit{T}}\right)$ follows the distribution

$$\begin{array}{}\text{(16)}& {\displaystyle }{\displaystyle }{\mathit{\rho }}_{\widehat{\mathrm{\Phi }}}\left(\widehat{\mathit{\varphi }}\right)=\int \mathit{\delta }\left(\mathit{\varphi }\right(\mathbf{\Delta }\widehat{\mathit{t}})-\widehat{\mathit{\varphi }}){\mathit{\rho }}_{\mathbf{\Delta }\widehat{\mathit{T}}}\left(\mathbf{\Delta }\widehat{\mathit{t}}\right)\mathrm{d}\mathbf{\Delta }\widehat{\mathit{t}}.\end{array}$$

Repeated application of Eq. (15) yields the uncertainty distribution of a given test's p value $\widehat{P}=p\left(\mathit{\varphi }\right(\mathbf{\Delta }\widehat{\mathit{T}}\left)\right)$:

$$\begin{array}{ll}{\displaystyle }{\mathit{\rho }}_{\widehat{P}}\left(\widehat{p}\right)& {\displaystyle }=\int \mathit{\delta }\left(p\left(\widehat{\mathit{\varphi }}\right)-\widehat{p}\right){\mathit{\rho }}_{\widehat{\mathrm{\Phi }}}\left(\widehat{\mathit{\varphi }}\right)\mathrm{d}\widehat{\mathit{\varphi }}\\ {\displaystyle }& {\displaystyle }=\int \int \mathit{\delta }\left(p\left(\widehat{\mathit{\varphi }}\right)-\widehat{p}\right)\mathit{\delta }\left(\mathit{\varphi }\left(\mathbf{\Delta }\widehat{\mathit{t}}\right)-\widehat{\mathit{\varphi }}\right)\\ {\displaystyle }& {\displaystyle }\times {\mathit{\rho }}_{\mathbf{\Delta }\widehat{\mathit{T}}}\left(\mathbf{\Delta }\widehat{\mathit{t}}\right)\mathrm{d}\mathbf{\Delta }\widehat{\mathit{t}}\mathrm{d}\widehat{\mathit{\varphi }}\\ \text{(17)}& {\displaystyle }& {\displaystyle }=\int \mathit{\delta }\left(p\left(\mathit{\varphi }\right(\mathbf{\Delta }\widehat{\mathit{t}}\left)\right)-\widehat{p}\right)\mathrm{d}\mathbf{\Delta }\widehat{\mathit{t}}.\end{array}$$

In the example shown in Fig. 4 the initial uncertainties in the observations translate into an uncertain p value that features both probability for significance and probability for non-significance. This illustrates the need for a criterion to project the uncertain p value onto a binary decision space comprised of rejection and acceptance of the null hypothesis. We propose to consider the following criteria to facilitate an informed decision:

• The hypothesis shall be rejected at the significance level α if and only if the expected p value is less than α, that is

  $$\begin{array}{}\text{(18)}& {\displaystyle }{\displaystyle }\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\widehat{p}{\mathit{\rho }}_{\widehat{P}}\left(\widehat{p}\right)\mathrm{d}\widehat{p}<\mathit{\alpha }.\end{array}$$

• The hypothesis shall be rejected at the significance level α if and only if the probability of p to be less than α is greater than a predefined threshold η (we propose η=90 %), that is

  $$\begin{array}{}\text{(19)}& {\displaystyle }{\displaystyle }\mathit{\pi }(\widehat{P}<\mathit{\alpha })=\underset{\mathrm{0}}{\overset{\mathit{\alpha }}{\int }}{\mathit{\rho }}_{\widehat{P}}\left(\widehat{p}\right)\mathrm{d}\widehat{p}>\mathit{\eta }.\end{array}$$

While the p value of a certain sample indicates its extremeness with respect to the null distribution, the expected p value may be regarded as a measure of the uncertain sample's extremeness. Given the measurement uncertainty, the quantity $\mathit{\pi }(\widehat{P}<\mathit{\alpha })$ constitutes an informed assessment of how likely or plausible the true value of the measured sample is to be statistically significant with respect to the null hypothesis. Thus, the first criterion assesses how “strongly” the uncertain sample contradicts the null hypothesis, while the second criterion evaluates the likelihood of the uncertain sample to contradict the null hypothesis. The choice of η is arbitrary and may be tailored to the specific circumstances of the test. In some situations, one might want to reject the null hypothesis only in the case of a high probability of significance and therefore choose a large value for η – e.g., when mistakenly attested significance is associated with high costs. In other situations, even small probabilities for significance may be important, e.g., if a significant test result would be associated with high costs or with high risks. We propose to assess the two criteria in combination. In the case of them not yielding the same result, the weighing between the criteria must be adapted to the statistical problem at hand.

3.3 Hypothesis tests

We have introduced the notion of uncertain samples and its consequences for the application of hypothesis tests. Here, we briefly introduce the tests used to test our null hypothesis that the observed tendency for delayed transition onsets in Na⁺ and δ¹⁸O with respect to Ca²⁺ and λ has occurred by chance and that the corresponding populations ${\mathcal{P}}_{\mathrm{\Delta }T}^{p,q}$ that characterize the pairwise random lags ΔT^p,q in fact do not favor the tentative transition orders apparent from the observations. Mathematically, this can be formulated as follows:

• Let ${\mathit{\rho }}_{\mathrm{\Delta }T}^{p,q}\left(\mathrm{\Delta }t\right)$ be the probability density associated with the population of DO transition onset lags ${\mathcal{P}}_{\mathrm{\Delta }{t}^{p,q}}$ between the proxy variables p and q and let the observations $\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}$ suggest a delayed transition of the proxy q – that is, the corresponding uncertainty distributions ${\mathit{\rho }}_{\mathrm{\Delta }{\widehat{T}}_i^{p,q}}\left(\mathrm{\Delta }{\widehat{t}}_i^{p,q}\right)$ indicate high probabilities for negative $\mathrm{\Delta }{\widehat{T}}_i^{p,q}$ across the sample according to Eq. (7). We then test the hypothesis H₀: The mean value ${\mathit{\mu }}^{p,q}=\int {\mathit{\rho }}_{\mathrm{\Delta }T}^{p,q}\left(\mathrm{\Delta }t\right)\mathrm{d}\mathrm{\Delta }t$ of the population ${\mathcal{P}}_{\mathrm{\Delta }T}^{p,q}$ is greater than or equal to zero.

We identified three tests that are suited for this task, namely the t test, the Wilcoxon signed-rank (WSR) test, and a bootstrap test. The WSR and the t test are typically formulated in terms of paired observation {x_i,y_i} that give rise to a sample of differences $\mathit{\{}{d}_i=x_i-y_i\mathit{\}}$ which correspond to the time lags $\mathit{\left\{}\mathrm{\Delta }{t}_i^{p,q}\mathit{\right\}}$ of different DO events (Rice, 2007; Lehmann and Romano, 2006, e.g.). The null distributions of the tests rely on slightly different assumptions regarding the populations. Since we cannot guarantee the compliance of these assumptions, we apply the tests in combination to obtain a robust assessment.

3.3.1 t test

The t test (Student, 1908) relies on the assumption that the population of differences 𝒫_D is normally distributed with mean μ and standard deviation σ. For a random sample $\mathit{D}=(D_{\mathrm{1}},\mathrm{\dots },D_n)$ the test statistic

$$\begin{array}{}\text{(20)}& {\displaystyle }{\displaystyle }Z\left(\mathit{D}\right)={\displaystyle \frac{U\left(\mathit{D}\right)-\mathit{\mu }}{S\left(\mathit{D}\right)/\sqrt{n}}}\end{array}$$

follows a t distribution t_n−1(z) with n−1 degrees of freedom. Here, $U=\frac{\mathrm{1}}{n}\sum D_i$ is the sample mean and $S=\frac{\mathrm{1}}{n-\mathrm{1}}\sum (U-D_i{)}^{\mathrm{2}}$ is the samples' standard deviation. This allows us to test whether an observed sample $\mathit{d}=(d_{\mathrm{1}},\mathrm{\dots },d_n)$ contradicts a hypothesis on the mean μ. To compute the p value for the hypothesis $H_{\mathrm{0}}:\mathit{\mu }\ge \mathrm{0}$ (left-handed application) the null distribution is integrated from −∞ to the observed value z(d):

$$\begin{array}{}\text{(21)}& {\displaystyle }{\displaystyle }{p}_z\left(z\right(\mathit{d}\left)\right)=\underset{-\mathrm{\infty }}{\overset{z\left(\mathit{d}\right)}{\int }}{t}_{n-\mathrm{1}}\left(z^{\prime }\right)\mathrm{d}{z}^{\prime }.\end{array}$$

The resulting p value must then be compared to the predefined significance level α.

The t test can be generalized for application to an uncertain sample of the form $\mathbf{\Delta }\widehat{\mathit{T}}=(\mathrm{\Delta }{\widehat{T}}_{\mathrm{1}},\mathrm{\dots },\mathrm{\Delta }{\widehat{T}}_n)$ as follows: let ${\mathit{\rho }}_{\mathbf{\Delta }\widehat{\mathit{T}}}\left(\mathbf{\Delta }\widehat{\mathit{t}}\right)$ denote the uncertainty distribution of $\mathbf{\Delta }\widehat{\mathit{T}}$. Then according to Eq. (15) the distribution of the uncertain statistic $\widehat{Z}\left(\mathbf{\Delta }\widehat{\mathit{T}}\right)$ reads

$$\begin{array}{}\text{(22)}& {\displaystyle }{\displaystyle }{\mathit{\rho }}_{\widehat{Z}}\left(\widehat{z}\right)=\int \mathit{\delta }\left({\displaystyle \frac{u\left(\mathbf{\Delta }\widehat{\mathit{t}}\right)}{s\left(\mathbf{\Delta }\widehat{\mathit{t}}\right)/\sqrt{n}}}-\widehat{z}\right){\mathit{\rho }}_{\mathbf{\Delta }\widehat{\mathit{T}}}\left(\mathbf{\Delta }\widehat{\mathit{t}}\right)\mathrm{d}\mathbf{\Delta }\widehat{\mathit{t}}.\end{array}$$

Finally, the distribution of the uncertain p value may again be computed according to Eq. (15)

$$\begin{array}{ll}{\displaystyle }{\mathit{\rho }}_{{\widehat{P}}_z}\left({\widehat{p}}_z\right)& {\displaystyle }=\int \mathit{\delta }\left(p_z\left(\widehat{z}\right)-{\widehat{p}}_z\right){\mathit{\rho }}_{\widehat{Z}}\left(\widehat{z}\right)\mathrm{d}\widehat{z}\\ \text{(23)}& {\displaystyle }& {\displaystyle }=\int \mathit{\delta }\left(\underset{-\mathrm{\infty }}{\overset{\widehat{z}}{\int }}{t}_{n-\mathrm{1}}\left(z\right)\mathrm{d}z-{\widehat{p}}_z\right){\mathit{\rho }}_{\widehat{Z}}\left(\widehat{z}\right)\mathrm{d}\widehat{z}\end{array}$$

and then be evaluated according to the two criteria formulated above.

3.3.2 Wilcoxon signed rank

Compared to the t test, the WSR test (Wilcoxon, 1945) allows us to relax the assumption of normality imposed on the generating population 𝒫_D and replaces it by the weaker assumption of symmetry with respect to its mean μ in order to test the null hypothesis $H_{\mathrm{0}}:\mathit{\mu }\ge \mathrm{0}$. The test statistic W for this test is defined as

$$\begin{array}{}\text{(24)}& {\displaystyle }{\displaystyle }W\left(\mathit{D}\right)=\sum _{i=\mathrm{1}}^{n}R\left(\right|D_i\left|\right)\mathbf{\Theta }\left(D_i\right),\end{array}$$

where $R\left(\right|D_i\left|\right)$ denotes the rank of $\left|D_i\right|$ within the sorted set of the absolute values of differences $\mathit{\left\{}\right|D_i\left|\mathit{\right\}}$. The Heaviside function Θ(D_i) guarantees that exclusively D_i>0 values are summed. The derivation of the null distribution is a purely combinatorial problem and its explicit form can be found in lookup tables. Because $W\in {\mathbb{N}}_{[\mathrm{0},n(n+\mathrm{1})/\mathrm{2}]}$ we denote the null distribution by ${\mathcal{P}}_W^{\mathrm{0}}\left(w\right)$ to signal that this is not a continuous density. Explicitly, the null distribution can be derived as follows: first, the assumption of symmetry around 0 (for the hypothesis $H_{\mathrm{0}}:\mathit{\mu }\ge \mathrm{0}$, the relevant null distribution builds on μ=0) guarantees that the chance for D_i to be positive is equal to $\frac{\mathrm{1}}{\mathrm{2}}$. Hence, the number of positive outcomes m follows a symmetric binomial distribution $\mathit{\pi }\left(m\right)=\left(\genfrac{}{}{0ex}{}{n}{m}\right)(\frac{\mathrm{1}}{\mathrm{2}}{)}^n$. For m positive observations, there are $\left(\genfrac{}{}{0ex}{}{n}{m}\right)$ different sets of ranks $\mathit{\{}{r}_{\mathrm{1}},\mathrm{\dots },r_m\mathit{\}}$ that they may assume and which are again, due to the symmetry of 𝒫_D, equally likely. Hence, for a given number of positive outcomes m the probability to obtain a test statistic w is given by the share of those $\left(\genfrac{}{}{0ex}{}{n}{m}\right)$ configurations that yield a rank sum equal to w. Summing these probabilities over all possible values of m yields the null distribution for the test statistic w.

For a given sample d we test the hypothesis $H_{\mathrm{0}}:\mathit{\mu }\ge \mathrm{0}$ by computing the corresponding one-sided p value p_w, which is given by the cumulative probability that the null distribution assigns to w^′ values smaller than the observed w(d):

$$\begin{array}{}\text{(25)}& {\displaystyle }{\displaystyle }{p}_w\left(w\right(\mathit{d}\left)\right)=\sum _{i=\mathrm{1}}^n{\mathcal{P}}_W^{\mathrm{0}}\left(w_i^{\prime }\right)\mathbf{\Theta }\left(w\right(\mathit{d})-w_i^{\prime }).\end{array}$$

Since $W\in {\mathbb{N}}_{[\mathrm{0},n(n+\mathrm{1})/\mathrm{2}]}$ it follows that P_w assumes only discrete values in [0,1] with the null distribution determining the mapping between these two sets.

The generalization of the WSR test to the uncertain sample $\mathbf{\Delta }\widehat{\mathit{T}}$ can be carried out almost analogously to the t test. However, the fact that $W\in {\mathbb{N}}_{[\mathrm{0},n(n+\mathrm{1})/\mathrm{2}]}$ makes it inconvenient to use a continuous probability density distribution. We denote the distribution for the uncertain $\widehat{W}\left(\mathbf{\Delta }\widehat{\mathit{T}}\right)$ by

$$\begin{array}{}\text{(26)}& \begin{array}{rl}{\mathcal{P}}_{\widehat{W}}\left(\widehat{w}\right)& =\int \mathit{\delta }\left(\sum _{i=\mathrm{1}}^{n}R\left(\right|\mathrm{\Delta }{\widehat{t}}_i\left|\right)\mathbf{\Theta }\left(\mathrm{\Delta }{\widehat{t}}_i\right)-\widehat{w}\right)\\ & {\mathit{\rho }}_{\mathbf{\Delta }\widehat{\mathit{T}}}\left(\mathbf{\Delta }\widehat{\mathit{t}}\right)\mathrm{d}\mathbf{\Delta }\widehat{\mathit{t}}.\end{array}\end{array}$$

Given the one-to-one map from all $w\in {\mathbb{N}}_{[\mathrm{0},n(n+\mathrm{1})/\mathrm{2}]}$ to the set of discrete potential values p_w for P_w in [0,1] determined by Eq. (25), the probability to obtain ${\widehat{p}}_w$ is already given by the probability to obtain the corresponding $\widehat{w}$. Hence, we find

$$\begin{array}{}\text{(27)}& {\displaystyle }{\displaystyle }{\mathcal{P}}_{{\widehat{P}}_w}\left(p_w\right(\widehat{W})={\widehat{p}}_w)={\mathcal{P}}_{\widehat{W}}\left(\widehat{w}\right).\end{array}$$

3.3.3 Bootstrap test

Given an observed sample of differences $\mathit{d}=(d_{\mathrm{1}},\mathrm{\dots },d_n)$, a bootstrap test constitutes a third option to test the compatibility of the sample with the hypothesis that the population of differences features a mean equal to or greater than 0: $H_{\mathrm{0}}:={\mathit{\mu }}_{\mathrm{0}}\ge \mathrm{0}$. Guidance for the construction of a bootstrap hypothesis test can be found in Lehmann and Romano (2006) and Hall and Wilson (1991). The advantage of the bootstrap test lies in its independence from assumptions regarding the distributions' shape. Lehmann and Romano (2006) propose the test statistic

$$\begin{array}{}\text{(28)}& {\displaystyle }{\displaystyle }v=\sqrt{n}u,\end{array}$$

with $u\left(\mathit{d}\right)=\frac{\mathrm{1}}{n}{\sum }_{i=\mathrm{1}}^n{d}_i$ denoting the sample mean. In contrast to the above two tests, the bootstrap test constructs the null distribution directly from the observed data. In the absence of assumptions, the best available approximation of the population 𝒫_D is given by the empirical density

$$\begin{array}{}\text{(29)}& {\displaystyle }{\displaystyle }{\mathcal{P}}_D\left(d\right)\sim {\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n\mathit{\delta }(d-d_i).\end{array}$$

However, the empirical density does not necessarily comply with the null hypothesis and it thus has to be shifted accordingly:

$$\begin{array}{}\text{(30)}& {\displaystyle }{\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_D\left(d\right)=\sum _{i=\mathrm{1}}^n\mathit{\delta }(d-d_i+u).\end{array}$$

${\stackrel{\mathrm{̃}}{\mathit{\rho }}}_D\left(d\right)$ corresponds to the borderline case of the null hypothesis μ=0. The null distribution for v is then derived by resampling m synthetic samples ${\stackrel{\mathrm{̃}}{\mathit{d}}}_j=({\stackrel{\mathrm{̃}}{d}}_{\mathrm{1}},\mathrm{\dots },{\stackrel{\mathrm{̃}}{d}}_n{)}_j$ of size n from ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}_D\left(d\right)$ and computing ${\stackrel{\mathrm{̃}}{v}}_j=v\left({\stackrel{\mathrm{̃}}{\mathit{d}}}_j\right)$ for each of them. This corresponds to randomly drawing n values from the set d−u with replacement and computing v for the resampled vectors m times, where the index j labels the iteration of this process. The resulting set $\mathit{\{}{\stackrel{\mathrm{̃}}{v}}_j{\mathit{\}}}_j$ induces the data-driven null distribution for the test statistic

$$\begin{array}{}\text{(31)}& {\displaystyle }{\displaystyle }{\mathit{\rho }}_V^{\mathrm{0}}\left(v\right)={\displaystyle \frac{\mathrm{1}}{m}}\sum _{j=\mathrm{1}}^m\mathit{\delta }(v-{\stackrel{\mathrm{̃}}{v}}_j).\end{array}$$

Setting m=10 000 we obtain robust null distributions for the two cases relevant here (n=16 and n=20). The p value of this bootstrap test is then computed as before in a one-sided manner:

$$\begin{array}{}\text{(32)}& {\displaystyle }{\displaystyle }{p}_v\left(v\right(\mathit{d}\left)\right)=\underset{-\mathrm{\infty }}{\overset{v\left(\mathit{d}\right)}{\int }}{\mathit{\rho }}_V^{\mathrm{0}}\left(v\right)\mathrm{d}v={\displaystyle \frac{\mathrm{1}}{m}}\sum _{j=\mathrm{1}}^m\mathbf{\Theta }\left(v\left(\mathit{d}\right)-{\stackrel{\mathrm{̃}}{v}}_j\right),\end{array}$$

where the right-hand side equals the fraction of resampled ${\stackrel{\mathrm{̃}}{v}}_j$ that are smaller than v(d) of the original sample.

In the case where the sample of differences is uncertain, as for $\mathbf{\Delta }\widehat{\mathit{T}}=(\mathrm{\Delta }{\widehat{T}}_{\mathrm{1}},\mathrm{\dots },\mathrm{\Delta }{\widehat{T}}_n)$, the construction scheme for ${\mathit{\rho }}_V^{\mathrm{0}}$ needs to be adjusted to reflect these uncertainties. In principle, each possible value $\mathbf{\Delta }\widehat{\mathit{t}}$ for the uncertain $\mathbf{\Delta }\widehat{\mathit{T}}$ is associated with its own null distribution ${\mathit{\rho }}_V^{\mathrm{0}}(v,\mathbf{\Delta }\widehat{\mathit{t}})$. In this sense, the value for the test statistic $v\left(\mathbf{\Delta }\widehat{\mathit{t}}\right)$ should be compared to the corresponding ${\mathit{\rho }}_V^{\mathrm{0}}(v,\mathbf{\Delta }\widehat{\mathit{t}})$ to derive a p value for this $\mathbf{\Delta }\widehat{\mathit{t}}$. Equations (31) and (32) define a mapping from $\mathbf{\Delta }\widehat{\mathit{t}}$ to its corresponding p value. To compute the uncertainty distribution for the p value, this map has to be evaluated for all potential $\mathbf{\Delta }\widehat{\mathit{t}}$, weighted by the uncertainty distribution ${\mathit{\rho }}_{\mathbf{\Delta }\widehat{\mathit{T}}}\left(\mathbf{\Delta }\widehat{\mathit{t}}\right)$:

$$\begin{array}{}\text{(33)}& {\displaystyle }{\displaystyle }{\mathit{\rho }}_{{\widehat{P}}_v}\left({\widehat{p}}_v\right)=\int \mathit{\delta }({\widehat{p}}_v-p_v(\mathbf{\Delta }\widehat{\mathit{t}}\left)\right){\mathit{\rho }}_{\mathbf{\Delta }\widehat{\mathit{T}}}\left(\mathbf{\Delta }\widehat{\mathit{t}}\right)\mathrm{d}\mathbf{\Delta }\widehat{\mathit{t}}.\end{array}$$

The three tests are applied in combination in order to compensate for their individual deficits. If the population 𝒫_ΔT was truly Gaussian, the t test would be the most powerful test; i.e., its rejection region would be the largest across all tests on the population mean (Lehmann and Romano, 2006). Since normality of 𝒫_ΔT cannot be guaranteed, the less powerful Wilcoxon signed-rank test constitutes a meaningful supplement to the t test, relying on the somewhat weaker assumption that 𝒫_ΔT is symmetric around 0. Finally, the bootstrap test is non-parametric and in view of its independence from any assumptions adds a valuable contribution.

3.4 Comparison to previous analysis

For the derivation of the transition lag uncertainty distributions ${\mathit{\rho }}_{\mathrm{\Delta }{\widehat{T}}_i^{p,q}}\left(\mathrm{\Delta }{\widehat{t}}_i^{p,q}\right)$ of the ith DO event between the proxies p and q, we have directly adopted the methodology designed by Erhardt et al. (2019). However, our statistical interpretation of the resulting sets of uncertainty distributions $\mathit{\left\{}{\mathit{\rho }}_{\mathrm{\Delta }{\widehat{T}}_{\mathrm{1}}^{p,q}}\right(\mathrm{\Delta }{\widehat{t}}_{\mathrm{1}}^{p,q}),\mathrm{\dots },{\mathit{\rho }}_{\mathrm{\Delta }{\widehat{T}}_n^{p,q}}(\mathrm{\Delta }{\widehat{t}}_n^{p,q}\left)\mathit{\right\}}$ derived from the set of DO events differs from the one proposed by Erhardt et al. (2019). In this section we explain the subtle yet important differences between the two statistical perspectives.

Given a pair of variables (p,q), Erhardt et al. (2019) define what they call “combined estimate” ${\mathit{\rho }}_{\mathrm{\Delta }{T}^{*}}\left(\mathrm{\Delta }{t}^{*}\right)$ as the product over all corresponding lag uncertainty distributions:

$$\begin{array}{}\text{(34)}& {\displaystyle }{\displaystyle }{\mathit{\rho }}_{\mathrm{\Delta }{T}^{*}}\left(\mathrm{\Delta }{t}^{*}\right)\propto \prod _{i=\mathrm{1}}^n{\mathit{\rho }}_{\mathrm{\Delta }{\widehat{T}}_i}\left(\mathrm{\Delta }{t}^{*}\right).\end{array}$$

This implicitly assumes that all DO events share the exact same time lag Δt^* between the variables p and q. This is realized by inserting a single argument Δt^* into the different distributions ${\mathit{\rho }}_{\mathrm{\Delta }{\widehat{T}}_i}(\cdot )$. Hence, the product on the right-hand side of Eq. (34) in fact indicates the probability that all DO events assume the time lag Δt^*, provided that they all assume the same lag:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\rho }}_{\mathrm{\Delta }{T}^{*}}\left(\mathrm{\Delta }{t}^{*}\right)={\mathit{\rho }}_{\mathrm{\Delta }{T}^{*}}\left(\mathrm{\Delta }{t}^{*}\right|\mathrm{\Delta }{\widehat{t}}_{\mathrm{1}}=\mathrm{\dots }=\mathrm{\Delta }{\widehat{t}}_n=\mathrm{\Delta }{t}^{*})\\ \text{(35)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\prod {\mathit{\rho }}_{\mathrm{\Delta }{\widehat{T}}_i}\left(\mathrm{\Delta }{t}^{*}\right)}{{\int }_{\mathrm{\Omega }}\prod {\mathit{\rho }}_{\mathrm{\Delta }{\widehat{T}}_i}\left(\mathrm{\Delta }{\widehat{t}}_i\right)\mathrm{d}\mathrm{\Delta }{\widehat{t}}_i}},\mathrm{\Omega }=\mathit{\{}\mathbf{\Delta }\widehat{\mathit{t}}:\mathrm{\Delta }{\widehat{t}}_i=\mathrm{\Delta }{\widehat{t}}_j\forall i,j\mathit{\}}.\end{array}$$

The denominator on the right-hand side equals the probability that all DO events share a common time lag. Equation (34) strongly emphasizes those regions where all uncertainty distributions ${\mathit{\rho }}_{\mathrm{\Delta }{\widehat{T}}_i}\left(\mathrm{\Delta }{\widehat{t}}_i\right)$ are simultaneously substantially larger than 0. The combined estimate answers the question: provided that all DO events exhibit the same lag between the transition onsets of p and q, then how likely is it that this lag is given by Δt^*. Drawing on this quantity, Erhardt et al. (2019) conclude that δ¹⁸O and Na⁺ “on average” lag Ca²⁺ and λ by about 1 decade.

Thinking of the DO transition onset lags as i.i.d. random variables of a repeatedly executed random experiment takes into account the natural variability between different DO events, and hence it removes the restricting a priori assumption $\mathrm{\Delta }{\widehat{t}}_{\mathrm{1}}=\mathrm{\dots }=\mathrm{\Delta }{\widehat{t}}_n$. In our approach we have related the potentially systematic character of lags to the population mean. Since the sample mean is the best point estimate of a population mean, we consider it to reasonably indicate potential leads and lags, whose significance should be tested in a second step. Thus, we ascribe to the sample mean a similar role as Erhardt et al. (2019) ascribe to the combined estimate, and therefore we present a comparison of these two quantities in Sect. 4.1.

The mean of an uncertain sample $\widehat{U}=u\left(\mathbf{\Delta }\widehat{\mathit{T}}\right)$ is again an uncertain quantity and its distribution reads

$$\begin{array}{}\text{(36)}& {\displaystyle }{\displaystyle }{\mathit{\rho }}_{\widehat{U}}\left(\widehat{u}\right)=\int \mathit{\delta }(\widehat{u}-u(\mathbf{\Delta }\widehat{\mathit{t}}\left)\right){\mathit{\rho }}_{\mathbf{\Delta }\widehat{\mathit{T}}}\left(\mathbf{\Delta }\widehat{\mathit{t}}\right)\mathrm{d}\mathbf{\Delta }\widehat{\mathit{t}}.\end{array}$$

While the combined estimate multiplies the distributions ${\mathit{\rho }}_{\mathrm{\Delta }{\widehat{T}}_i}\left(\mathrm{\Delta }{t}^{*}\right)$, the uncertain sample mean convolutes them pairwise (see Appendix C). We thus expect the distributions for uncertain sample means to be broader than the corresponding distributions for the combined estimate. This can be motivated by considering the simple example of two Gaussian variables X and Y. According to the convolution their sample mean $U=\frac{X+Y}{\mathrm{2}}$ is normally distributed with variance ${\mathit{\sigma }}_{x*y}^{\mathrm{2}}=\frac{{\mathit{\sigma }}_x^{\mathrm{2}}+{\mathit{\sigma }}_y^{\mathrm{2}}}{\mathrm{4}}$ . In contrast, a combined estimate would yield a normal distribution with variance ${\mathit{\sigma }}_{xy}^{\mathrm{2}}=\frac{{\mathit{\sigma }}_x^{\mathrm{2}}+{\mathit{\sigma }}_y^{\mathrm{2}}}{{\mathit{\sigma }}_x^{\mathrm{2}}{\mathit{\sigma }}_y^{\mathrm{2}}}$. Thus, the convolution will appear broader for all ${\mathit{\sigma }}_x^{\mathrm{2}}{\mathit{\sigma }}_y^{\mathrm{2}}>\mathrm{4}$, which is the case for the distributions considered in this study.

4 Results

In the following we apply the above methodology to the different pairs of proxies that Erhardt et al. (2019) found to exhibit a decadal-scale time lag, based on an assessment of the combined estimate, namely $({\mathrm{Ca}}^{\mathrm{2}+},{\mathrm{Na}}^{+})$, $(\mathit{\lambda },{\mathrm{Na}}^{+})$, $({\mathrm{Ca}}^{\mathrm{2}+},{\mathit{\delta }}^{\mathrm{18}}\mathrm{O})$, and (λ,δ¹⁸O) from the NGRIP ice core and $({\mathrm{Ca}}^{\mathrm{2}+},{\mathrm{Na}}^{+})$ from the NEEM ice core. For each individual proxy we estimate the uncertain transition onsets relative to the timing of the DO events as given by Rasmussen et al. (2014) (see Fig. 2). From these uncertain transition onsets, the uncertainty distributions for the sets of uncertain lags $\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}$ between the proxies p and q are derived according to Eq. (7). As mentioned previously, we study the same selection of transitions evidenced in the multi-proxy records as Erhardt et al. (2019). This selection yields sample sizes of either 16 or 20 lags per pair of proxies but not 23, which is the total number of DO events present in the data.

We first study the uncertain sample means. As already mentioned, the sample mean is the best available point estimate for the population mean. Hence, sample means different from 0 may be regarded as first indications for potential systematic lead–lag relationships and thus motivate the application of hypothesis tests. We compare the results obtained for the uncertain sample means with corresponding results for the combined estimate. Both quantities indicate a tendency towards a delayed transition in Na⁺ and δ¹⁸O. Accordingly, in the subsequent section we apply the generalized hypothesis tests introduced above to the uncertain samples of transition lags to test the null hypothesis that pairwise, the apparent transition sequence is not systematically favored, that is, that the populations have mean equal to or greater than 0.

4.1 Uncertain sample mean and combined estimate

Based on their assessment of the combined estimate, Erhardt et al. (2019) concluded that on average, transitions in Ca²⁺ and λ started approximately 1 decade earlier than their counterparts in Na⁺ and δ¹⁸O. Figure 5 shows a reproduction of their results together with the uncertainty distributions of the sample means for all proxy pairs under study ($({\mathrm{Ca}}^{\mathrm{2}+},{\mathit{\delta }}^{\mathrm{18}}\mathrm{O})$ and (λ,δ¹⁸O) are not shown in Erhardt et al., 2019). For an uncertain sample of lags $\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}$ between the proxies p and q, the combined estimate and the uncertain sample mean are computed according to Eqs. (35) and (36), respectively. The reproduction of the combined estimate deviates from the original publication by no more than 1 year with respect to the mean and the 5th and 95th percentiles across all pairs. These deviations might originate from the stochastic MCMC-sampling process used for the analysis.

Figure 5 Comparison between the uncertain sample means (this study) and “combined estimates” according to Erhardt et al. (2019). The probability densities for the combined estimate are derived from the samples of uncertain time lags according to Eq. (34). Correspondingly, the uncertain sample means are computed according to Eq. (36). The numbers in the plots indicate the mean, the 5th, and the 95th percentile of the respective quantity. Both computations use Gaussian kernel density estimates of the MCMC-sampled transition onsets lags. Panels (a–d) refer to proxy pairs from the NGRIP ice core, and panel (e) shows results from the NEEM ice core. The distributions for both the combined estimate and the uncertain sample mean point towards a delayed transition onset in δ¹⁸O and Na⁺ with respect to λ and Ca²⁺.

With the sample mean being the best point estimator of the population mean, it serves as a suitable indicator for a potential population mean different from 0. The expectations

$$\begin{array}{}\text{(37)}& {\displaystyle }{\displaystyle }\mathrm{E}\left(\widehat{U}\right)=\int \widehat{u}{\mathit{\rho }}_{\widehat{U}}\left(\widehat{u}\right)\mathrm{d}\widehat{u}\end{array}$$

for the sample means of all proxy pairs do in fact suggest a tendency towards negative values in all distributions, i.e., a delay of the Na⁺ and δ¹⁸O transition onsets with respect to Ca²⁺ and λ. This indication is weakest for (${\mathrm{Ca}}^{\mathrm{2}+},{\mathrm{Na}}^{+}$) and (${\mathrm{Ca}}^{\mathrm{2}+},{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$) from NGRIP, since for these pairs we find a non-zero probability of a positive sample mean. For the other pairs the indication is comparably strong, with the 95th percentiles of the uncertainty distributions for the sample mean still being less than 0. Overall, the results for the uncertain mean confirm the previously reported tendencies, and in very rough terms, the distributions qualitatively agree with those for the combined estimate. In agreement with the heuristic example from Sect. 3.4, we find the sample mean distributions to be broader than the combined estimate distributions in all cases. The expected sample means indicate less pronounced lags for $({\text{Ca}}^{\mathrm{2}},{\mathrm{Na}}^{+})$ (Fig. 5a) and $({\mathrm{Ca}}^{\mathrm{2}+},{\mathit{\delta }}^{\mathrm{18}}\mathrm{O})$ (Fig. 5c) from the NGRIP ice core compared to the expectations of the corresponding combined estimate. In combination with the broadening of the distribution, this yields considerable probabilities for U>0 of 12 % and 14 %, respectively, indicating a delayed transition of Ca²⁺ in the sample mean with respect to Na⁺ or δ¹⁸O. Contrarily, for $(\mathit{\lambda },{\mathrm{Na}}^{+})$ (NGRIP, Fig. 5b) and $({\mathrm{Ca}}^{\mathrm{2}+},{\mathrm{Na}}^{+})$ (NEEM, Fig. 5e) the expected sample means point towards more distinct lags than reported by Erhardt et al. (2019) based on the combined estimate. For (λ,δ¹⁸O) (NGRIP, Fig. 5d) the sample mean and the combined estimate are very close. Note that the analysis of the uncertain sample values yields a more inconsistent picture with regard to the $({\mathrm{Ca}}^{\mathrm{2}+},{\mathrm{Na}}^{+})$ lag in the two different cores. While the distribution is shifted to less negative (less pronounced lag) for the NGRIP data, it tends to more negative values in the case of NEEM (stronger lag), suggesting a slight discrepancy between the cores.

Both quantities, the uncertain sample mean and the combined estimate point towards delayed transition onsets in Na⁺ and δ¹⁸O with respect to Ca²⁺ and λ, with major fractions of their uncertainty densities being allocated to negative values. This provides motivation to test whether the observations significantly contradict the hypothesis of a population mean equal to or greater than 0. Accordingly, the subsequent section presents the results obtained from the application of three different hypothesis tests that target the population mean. As discussed in Sect. 3, the tests have been modified to allow for a rigorous uncertainty propagation and return an uncertainty distribution for their corresponding p values rather than scalars.

4.2 Statistical significance of the proposed lead–lag relations

Above, we identified three tests for testing the hypothesis that the samples $\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}$ were actually generated from populations that on average feature no or even reversed time lags compared to what the sign of the corresponding uncertain sample mean suggests. Mathematically, this is equivalent to testing the hypothesis that the mean μ^p,q of the population ${\mathcal{P}}_{\mathrm{\Delta }T}^{p,q}$ is greater than or equal to 0: $H_{\mathrm{0}}:{\mathit{\mu }}^{p,q}\ge \mathrm{0}$. A rejection of this hypothesis would confirm that the assessed sample is very unlikely to stem from a population with ${\mathit{\mu }}^{p,q}\ge \mathrm{0}$ and would thereby provide evidence for a systematic lag. Under the constraints indicated above this would in turn yield evidence for an actual lead of the corresponding climatic process. We have chosen a significance level of α=0.05, which is a typical choice. Figure 7 summarizes the final uncertainty distributions of the three tests for all proxy pairs under study. Corresponding values are given in Table 1.

Table 1 Results from the application of the t test, the WSR test, and a bootstrap test to uncertain samples of DO transition onset lags $\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}$. $\mathrm{E}\left(\widehat{P}\right)$ denotes the expected p value, derived from the uncertainty-propagated p-value distribution. The probability of significant test results associated with the same distribution is indicated by $\mathit{\pi }(\widehat{P}<\mathrm{0.05})$. For comparison, the p values from the application of the tests to the expected sample $\mathrm{E}\left(\mathbf{\Delta }\widehat{\mathit{T}}\right)=\int {\mathit{\rho }}_{\mathbf{\Delta }\widehat{\mathit{T}}}\left(\mathbf{\Delta }\widehat{\mathit{t}}\right)\mathbf{\Delta }\widehat{\mathit{t}}\mathrm{d}\mathbf{\Delta }\widehat{\mathit{t}}$ are given in the bottom row.

Figure 6 Exemplary application of the analysis to the proxy pair $({\mathrm{Ca}}^{\mathrm{2}+},{\mathrm{Na}}^{+})$ from the NGRIP ice core. Panel (a) shows 16 uncertain time lags $\mathrm{\Delta }{\widehat{T}}_i$ derived from the proxy data around DO events. The continuous densities have been obtained via a Gaussian kernel density estimate from the corresponding MCMC samples (see Sect. 3.1). In panel (b) the uncertain test statistics induced by the uncertain sample are shown for the t test (blue), the WSR test (orange), and a bootstrap test (green). The values that comprise the histograms are immediately derived from the MCMC samples. Panel (c) shows the empirical uncertainty distribution for the p values of the three tests, following from the uncertain test statistics in panel (b). Dotted red lines separate rejection from acceptance regions in panels (b) and (c). In the case of the bootstrap test, the rejection regions cannot be defined consistently on the level of the test statistic, since each possible value $\mathbf{\Delta }\widehat{\mathit{t}}$ for the uncertain $\mathbf{\Delta }\widehat{\mathit{T}}$ generates its individual null distribution. The null distribution shown here is in fact the pooled distribution of resampled ${\stackrel{\mathrm{̃}}{v}}_j$ obtained from all MCMC- sampled values for $\mathbf{\Delta }\widehat{\mathit{T}}$. For the other proxy pairs investigated in this study, corresponding plots would appear structurally similar.

Figure 6 exemplarily illustrates the application of the three tests to the empirical densities obtained for $\mathbf{\Delta }\widehat{\mathit{T}}({\mathrm{Ca}}^{\mathrm{2}+},{\mathrm{Na}}^{+})$ (NGRIP). In Fig. 6 the initial uncertainty in the observations – i.e., the uncertainty encoded by the distributions of transition onset lags – is propagated to an uncertain test statistic according to Eq. (16). In turn, the uncertain test statistic yields an uncertain p value (see Eq. 17). Since the numerical computation is based on empirical densities as generated by the MCMC sampling, we show the corresponding histograms instead of continuous densities – for ${\mathit{\rho }}_{\mathrm{\Delta }{\widehat{T}}_i}\left(\mathrm{\Delta }{\widehat{t}}_i\right)$, Gaussian kernel density estimates are presented only for the sake of visual clarity. On the level of the test statistics the red dashed line separates the acceptance from the rejection region, based on the null distributions given in black. Qualitatively, the three tests yield the same results. The histograms clearly indicate non-zero probabilities for the test statistic in both regions. Correspondingly, the histograms for the p values stretch well across the significance threshold. The shapes of the histograms resemble an exponential decay towards higher p values. This results from the non-linear mapping of the test statistics to the p values. Despite the pronounced bulk of empirical p values below the significance level, the probability of non-significant p values is still well above 50 % for the three tests (see Table 1). Also, the expected p value exceeds the significance level for all tests. Hence, neither of the two criteria for rejecting the null hypothesis formulated in Sect. 3.2 is met for the proxy pair $({\mathrm{Ca}}^{\mathrm{2}+},{\mathrm{Na}}^{+})$. In contrast, if the observational uncertainties are averaged out on the level of the transition onset lags, all tests yield p values below the significance level, which would indicate that the lags were indeed significant. Hence, the rigorous propagation of uncertainties qualitatively changes the statistical assessment of the uncertain sample of lags $\mathbf{\Delta }\widehat{\mathit{T}}({\mathrm{Ca}}^{\mathrm{2}+},{\mathrm{Na}}^{+})$ (NGRIP). While the expected sample rejects the null hypothesis, rigorous uncertainty propagation leads to acceptance. This holds true for all tests.

Figure 7 Results of the hypothesis tests applied to the uncertain samples of transition onset lags $\mathbf{\Delta }{\widehat{\mathit{T}}}_i^{p,q}$. The violin plots show the Gaussian kernel density estimates of the empirical uncertainty distributions for p values (see Fig. 6) obtained for all tests and for all proxy pairs investigated. Pink bars indicate the corresponding expected p values $\mathrm{E}\left(\widehat{P}\right)$ and yellow bars indicate the p values obtained from testing the expected samples $\mathrm{E}\left(\mathbf{\Delta }\widehat{\mathit{T}}\right)$. All expected p values are above the significance level α=0.05 (red dotted line), while the expected samples appear to be significant consistently across all proxy pairs and all tests. The pie charts indicate the probability of the respective p values to be less then α.

Figure 7 summarizes the results obtained for all proxy pairs under study. Qualitatively, our findings are the same for the other pairs as for the $({\mathrm{Ca}}^{\mathrm{2}+},{\mathrm{Na}}^{+})$ (NGRIP) case discussed in detailed above. All expected p values, as indicated by the pink bars, are above the significance level. Also, the probability for significance is below 60 % for all pairs and all tests as shown by the pie charts. Therefore, for all proxy pairs and for all tests, the formulated decision criteria do thus not allow us to reject the null hypothesis of a population mean greater than or equal to 0. In contrast, all expected samples are significant across all tests with corresponding p values indicated by the yellow bars. The proxy pairs with the lowest expected p values and the highest probability of $\widehat{P}<\mathit{\alpha }$ are $(\mathit{\lambda },{\mathrm{Na}}^{+})$ from NGRIP and $({\mathrm{Ca}}^{\mathrm{2}+},{\mathrm{Na}}^{+})$ from NEEM, as already suggested by the analysis of the uncertain sample mean. For the NGRIP ice core the delay of Na⁺ and δ¹⁸O with respect to Ca²⁺ has a very low probability to be significant (approximately one-third). The pair (λ,δ¹⁸O) ranges in between the latter two.

5 Discussion

Erhardt et al. (2019) have reported an average time lag between the transition onsets in Na⁺ and δ¹⁸O proxy values and their counterparts in Ca²⁺ and λ at the onset of DO events. This statement is based on the assessment of the combined estimate derived from uncertain samples of time lags $\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}$. The samples were obtained by applying a well-suited Bayesian transition onset detection scheme to high-resolution time series of the different proxies. The combined estimate indicates leads of the Ca²⁺ and λ transition onsets with respect to Na⁺ and δ¹⁸O by approximately 1 decade, with the 90 % confidence interval ranging from 0 to approximately 15 years. The combined estimate implicitly assumes that for a given proxy pair all DO events share a common time lag ($\mathrm{\Delta }{\widehat{T}}_i^{p,q}=\mathrm{\Delta }{\widehat{T}}_j^{p,q}$).

We argue that the variability across different DO events cannot be ignored in the assessment of the data. Although the DO events are likely to be caused by the same physical mechanism, changing boundary conditions and other natural climate fluctuations will lead to deviations in the exact timings of the different processes involved in triggering the individual DO events. Figure 2 clearly shows that the different events exhibit different time lags. Provided that the DO events were driven by the same process, physically they constitute different realizations, and they exhibit great variability also in other variables such as the amplitude of the temperature change (Kindler et al., 2014) or the waiting times with respect to the previous event (Ditlevsen et al., 2007; Boers et al., 2018). The random experiment framework introduced in this study allows us to relax the constraint of a common time lag Δt^* shared across all events and reflects the fact that natural variability will cause different expressions of the same mechanism across different DO events. Moreover, this framework relates potential systematic leads and lags in the physical process that drives DO events to a corresponding non-zero mean of a population of lags between proxy variables. This allows for the physically meaningful formulation of a statistical hypothesis and a corresponding null hypothesis. By applying different hypothesis tests we have followed a well-established line of statistical inference. Motivated by the apparent transition onset delays in Na⁺ and δ¹⁸O with respect to the transitions in λ and Ca²⁺, as reported by Erhardt et al. (2019) and confirmed here on the level of uncertain sample means, we tested the null hypothesis that the corresponding populations do not favor the proposed transition sequence. Rejection of this hypothesis would have provided evidence that the observed lag tendency is an imprint of the underlying physical process and therefore a systematic feature. However, generalized versions of three different hypothesis tests consistently fail to reject the null hypothesis under rigorous propagation of the observational uncertainties originating from the MCMC-based transition onset detection. This holds true for all proxy pairs. The fact that the tests rely on different assumptions on the population's shape but nonetheless qualitatively yield the same results makes our assessment robust. We conclude that the possibility that the observed tendencies towards advanced transitions in Ca²⁺ and λ have occurred simply by chance cannot be ruled out. If the common physical interpretation of the studied proxies holds true, our results imply that the hypothesis that the trigger of the DO events is associated directly with the North Atlantic sea-ice cover rather than the atmospheric circulation – be it on a synoptic or hemispheric scale – cannot be ruled out. We emphasize that our results should not be misunderstood as evidence against the alternative hypothesis of a systematic lag. In the presence of a systematic lag (μ<0) the ability of hypothesis tests to reject the null hypothesis of no systematic lag ($(H_{\mathrm{0}}:\mathit{\mu }=\mathrm{0})$) depends on the sample size n, the ratio between the mean lag $\left|\mathit{\mu }\right|$, the variance of the population, and the precision of the measurement. Neither of these quantities is favorable in our case, and thus it is certainly possible that the null hypothesis cannot be rejected despite the alternative being true.

Our main purpose was the consistent treatment of observational uncertainties and we have largely ignored the vibrant debate on the qualitative interpretation of the proxies. Surprisingly, we could not find any literature on the application of hypothesis tests to uncertain samples of the kind discussed here. The theory of fuzzy p values is in fact concerned with uncertainties either in the data or in the hypothesis. However, it is not applicable to measurement uncertainties that are quantifiable in terms of probability density functions (Filzmoser and Viertl, 2004). We have proposed to propagate the uncertainties to the level of the p values and to then consider the expected p values and the share of p values which indicate significance, in order to decide between rejection and acceptance. The p value measures the extremeness of a sample with respect to the null distribution, and we hence regard the expected p value to be a suitable measure for the uncertain samples' extremeness. In cases of a high cost of a wrongly rejected null hypothesis, one might want to have a high degree of certainty that the uncertain sample actually contradicts the null hypothesis and hence a high probability of the uncertain p value being smaller than α. In contrast, if the observational uncertainties are averaged out beforehand, crucial information is lost. The expected sample may either be significant or not, but the uncertainty about the significance can no longer be accurately quantified.

The potential of the availability of data from different sites has probably not been fully leveraged in this study. Naively, one could think of the NEEM and NGRIP $({\mathrm{Ca}}^{\mathrm{2}+},{\mathrm{Na}}^{+})$ lag records as two independent observations of the same entity. However, given the discrepancies in the corresponding sample mean uncertainty distributions, changes in sea-ice cover and atmospheric circulation could in fact have impacted both cores differently. Proxy-enabled modeling studies as presented by Sime et al. (2019) could shed further light on the similarity of the signals at NEEM and NGRIP as a function of change in climatic conditions. Also, a comparison of the NGRIP and NEEM records on an individual event level could provide further insights into how to combine these records statistically. There might be ways to further exploit the advantage of having two recordings of the same signal.

6 Conclusions

We have presented a statistical reinterpretation of the high-resolution proxy records provided and analyzed by Erhardt et al. (2019). The probabilistic transition onset detection also designed by Erhardt et al. (2019) very conveniently quantifies the uncertainty in the transition onset estimation by returning probability densities instead of scalar estimates. While the statistical quantities “combined estimate” (Erhardt et al., 2019) and “uncertain sample mean” (this study) indicate a tendency for a pairwise delayed transition onset in Na⁺ and δ¹⁸O proxy values with respect to Ca²⁺ and λ, a more rigorous treatment of the involved uncertainties shows that these tendencies are not statistically significant. That is, at the significance level α=5 % they do not contradict the null hypothesis that no or the reversed transition sequence is in fact physically favored. Thus, a pairwise systematic lead–lag relation cannot be evidenced for any of the proxies studied here. We have shown that if uncertainties on the level of transition onset lags are averaged out beforehand, the samples of lags indeed appear to be significant, which underpins the importance of rigorous uncertainty propagation in the analysis of paleoclimate proxy data. We have focused on the quantitative uncertainties and have largely ignored qualitative uncertainty stemming from the climatic interpretation of the proxies. However, if the common proxy interpretations hold true, our findings suggest that, for example, the hypothesis of an atmospheric trigger – either of hemispheric or synoptic scale – for the DO events should not be favored over the hypothesis that a change in the North Atlantic sea-ice cover initiates the DO events.

Even though we find that the uncertainty of the transition onset detection combined with the small sample size prevents the deduction of statistically unambiguous statements on the temporal order of events, we think that multi-proxy analysis is a promising approach to investigate the sequential order at the beginning of DO events. In this study, we refrained from analyzing the lags between the different proxies in a combined approach and focused on the marginal populations. However, a combined statistical evaluation – that is, treating the transition onsets of all proxy variables as a four-dimensional random variable – merits further investigation. Also, we propose to statistically combine measurements from NEEM and NGRIP (and potentially further ice cores) of the same proxy pairs. Finally, hierarchical models may be invoked to avoid switching from a Bayesian perspective in the transition onset estimation to a frequentist perspective in the statistical interpretation of the uncertain samples. Finally, effort in conducting modeling studies should be sustained. Especially proxy-enabled modeling bears the potential to improve comparability between model results and paleoclimate records. Together, these lines of research are promising to further constrain the sequence of events that have caused the abrupt climate changes associated with DO events.

Appendix A: Numerical treatment of high-dimensional probability densities

In Sect. 3.1 we introduced the probabilistic transition onset detection designed by Erhardt et al. (2019). Given a single time series, the formulation of a stochastic ramp model induces a posterior probability density for the set of model parameters θ in a Bayesian sense:​​​​​​​

$$\begin{array}{}\text{(A1)}& {\displaystyle }{\displaystyle }\mathit{\pi }\left(\mathit{\theta }\right|\mathcal{D})={\displaystyle \frac{\mathit{\pi }\left(\mathcal{D}\right|\mathit{\theta }\left)\mathit{\pi }\right(\mathit{\theta })}{\mathit{\pi }\left(\mathcal{D}\right)}}.\end{array}$$

However, a classical numerical representation of this density on a discretized grid is inconvenient. Due to its high dimensionality for a reasonable grid spacing the number of data points easily overloads the computational power of ordinary computers. For example, representing each dimension with a minimum of 100 points would amount to a total of 10¹² data points. On top of that, the application of any methods to such a grid is computationally very costly. Here, the MCMC sampler constitutes an efficient solution. By sampling a representative set {θ_j}_j from the posterior probability density it may be used to construct an empirical density in the sense of Eq. (4). For the sake of simplicity in the main text we have formulated the methods in terms of continuous probability densities, although all computations in fact rely on empirical densities obtained from MCMC samples. Here, we show that all steps in the derivation of the methods can be performed equivalently under stringent use of the empirical density. With regards to hypothesis tests, the use of empirical densities for the uncertain transition lag samples $\mathbf{\Delta }{\mathit{T}}_i^{p,q}$ essentially boils down to an application of the tests to each individual value comprised in the respective empirical density.

For a given proxy and a given DO event, in a first step the MCMC algorithm samples from the joint posterior probability density for the models parameter configuration $\mathit{\theta }=(t_{\mathrm{0}},\mathit{\tau },y_{\mathrm{0}},\mathrm{\Delta }y,\mathit{\alpha },\mathit{\sigma })$, giving rise to the empirical density ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathbf{\Theta }}\left(\mathit{\theta }\right)=\frac{\mathrm{1}}{m}\sum \mathit{\delta }(\mathit{\theta }-{\mathit{\theta }}_j)$. Integration over the nuisance parameters then yields the marginal empirical density for the transition onset:

$$\begin{array}{}\text{(A2)}& {\displaystyle }{\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{T_{\mathrm{0}}}^{p,i}\left(t_{\mathrm{0}}^{p,i}\right)={\displaystyle \frac{\mathrm{1}}{m}}\sum _{j=\mathrm{1}}^m\mathit{\delta }\left(t_{\mathrm{0}}^{p,i}-t_{\mathrm{0},j}^{p,i}\right),\end{array}$$

where the index i indicates the DO event and p denotes the proxy variable while j runs over the MCMC sampled values. We use bars to mark empirical densities in contrast to continuous densities. The uncertainty distribution for the lag $\mathrm{\Delta }{T}_i^{p,q}$ between the variables p and q as defined by Eq. (7) may then be approximated as follows (omitting the index i):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\rho }}_{\mathrm{\Delta }T}^{p,q}\left(\mathrm{\Delta }{t}^{p,q}\right)\\ {\displaystyle }& {\displaystyle }=\iint \mathit{\delta }\left(t_{\mathrm{0}}^p-t_{\mathrm{0}}^q-\mathrm{\Delta }{t}^{p,q}\right){\mathit{\rho }}_{T_{\mathrm{0}}}^p\left(t_{\mathrm{0}}^p\right){\mathit{\rho }}_{T_{\mathrm{0}}}^q\left(t_{\mathrm{0}}^q\right)\mathrm{d}{t}_{\mathrm{0}}^p\mathrm{d}{t}_{\mathrm{0}}^q\\ {\displaystyle }& {\displaystyle }\simeq \iint \mathit{\delta }\left(t_{\mathrm{0}}^p-t_{\mathrm{0}}^q-\mathrm{\Delta }{t}^{p,q}\right){\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{T_{\mathrm{0}}}^p\left(t_{\mathrm{0}}^p\right){\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{T_{\mathrm{0}}}^q\left(t_{\mathrm{0}}^q\right)\mathrm{d}{t}_{\mathrm{0}}^p\mathrm{d}{t}_{\mathrm{0}}^q\\ {\displaystyle }& {\displaystyle }=\iint \mathit{\delta }\left(t_{\mathrm{0}}^p-t_{\mathrm{0}}^q-\mathrm{\Delta }{t}^{p,q}\right){\displaystyle \frac{\mathrm{1}}{m}}\sum _{j=\mathrm{1}}^m\mathit{\delta }\left(t_{\mathrm{0}}^p-t_{\mathrm{0},j}^p\right)\\ {\displaystyle }& {\displaystyle }\times {\displaystyle \frac{\mathrm{1}}{m}}\sum _{k=\mathrm{1}}^m\mathit{\delta }\left(t_{\mathrm{0}}^p-t_{\mathrm{0},k}^q\right)\mathrm{d}{t}_{\mathrm{0}}^p\mathrm{d}{t}_{\mathrm{0}}^q\\ \text{(A3)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{m^{\mathrm{2}}}}\sum _{j,k=\mathrm{1}}^m\mathit{\delta }\left(t_{\mathrm{0},j}^p-t_{\mathrm{0},k}^q-\mathrm{\Delta }{t}^{p,q}\right)={\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{\Delta }T}^{p,q}\left(\mathrm{\Delta }{t}^{p,q}\right).\end{array}$$

Thus, the empirical uncertainty distribution for the time lag is induced by the set of all possible differences between members of the two MCMC samples for the respective transition onsets:

$$\begin{array}{}\text{(A4)}& {\displaystyle }{\displaystyle }{\left\{\mathrm{\Delta }{t}_j^{p,q}\right\}}_{j\in [\mathrm{1},m^{\mathrm{2}}]}={\left\{t_{\mathrm{0},k}^p-t_{\mathrm{0},l}^q\right\}}_{p,q\in [\mathrm{1},m]}.\end{array}$$

For this study m=6000 values have been sampled with the MCMC algorithm for each transition under study. This yields $m^{\mathrm{2}}=\mathrm{36}\times {\mathrm{10}}^{\mathrm{6}}$ potential values for the empirical ΔT uncertainty distribution. To keep the computation efficient, the sets of lags were restricted to combinations k=l and thus to 6000 empirical values. We thus approximate

$$\begin{array}{}\text{(A5)}& {\displaystyle }{\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{\Delta }T}^{p,q}\left(\mathrm{\Delta }{t}^{p,q}\right)\simeq {\displaystyle \frac{\mathrm{1}}{m}}\sum _{j=\mathrm{1}}^m\mathit{\delta }\left(t_{\mathrm{0},j}^p-t_{\mathrm{0},j}^q-\mathrm{\Delta }{t}^{p,q}\right).\end{array}$$

This drastic reduction in values certainly requires justification, which we give later by comparing final results of the analysis to those obtained from control runs. The control runs analogously construct the empirical densities for the transition onset lags from 6000 out of the 36×10⁶ possible values, but use randomly shuffled versions of the original sets of transition onset times for the variables p and q:

$$\begin{array}{}\text{(A6)}& {\displaystyle }{\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{\Delta }T}^{p,q,\text{ctrl}}\left(\mathrm{\Delta }{t}^{p,q}\right)\simeq {\displaystyle \frac{\mathrm{1}}{m}}\sum _{j=\mathrm{1}}^m\mathit{\delta }\left(t_{\mathrm{0},s\left(j\right)}^p-t_{\mathrm{0},s^{\prime }\left(j\right)}^q-\mathrm{\Delta }{t}^{p,q}\right).\end{array}$$

Here, s and s^′ denote randomly chosen permutations of the set $\mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots },m\mathit{\}}$.

As in the main text, in the following we denote uncertain quantities with a hat. For a given proxy pair the starting point for the statistical analysis however, is the uncertain sample $\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}=(\mathrm{\Delta }{\widehat{T}}_{\mathrm{1}}^{p,q},\mathrm{\dots },\mathrm{\Delta }{\widehat{T}}_n^{p,q})$, which is characterized by the n-dimensional uncertainty distribution ${\mathit{\rho }}_{\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}}\left(\mathbf{\Delta }{\widehat{\mathit{t}}}^{p,q}\right)=\prod {\mathit{\rho }}_{\mathrm{\Delta }{\widehat{T}}_i^{p,q}}\left(\mathrm{\Delta }{\widehat{t}}_i^{p,q}\right)$. Its empirical counterpart is given by

$$\begin{array}{ll}{\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}}\left(\mathbf{\Delta }{\widehat{\mathit{t}}}^{p,q}\right)& {\displaystyle }=\prod _{i=\mathrm{1}}^n{\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{\Delta }{\widehat{T}}_i^{p,q}}\left(\mathrm{\Delta }{\widehat{t}}_i^{p,q}\right)\\ \text{(A7)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{m^n}}\prod _{i=\mathrm{1}}^n\sum _{j=\mathrm{1}}^m\mathit{\delta }\left(\mathrm{\Delta }{\widehat{t}}_i^{p,q}-\mathrm{\Delta }{t}_{i,j}^{p,q}\right).\end{array}$$

This empirical density is comprised of mⁿ possible values for the n-dimensional random vector $\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}$, and again, a substantial reduction in the representing set is required for practical computation. Defining the reduced empirical density for $\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}$ as

$$\begin{array}{ll}{\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_{\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}}\left(\mathbf{\Delta }{\widehat{\mathit{t}}}^{p,q}\right)& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{m}}\sum _{j=\mathrm{1}}^m\prod _{i=\mathrm{1}}^n\mathit{\delta }\left(\mathrm{\Delta }{\widehat{t}}_i^{p,q}-\mathrm{\Delta }{t}_{i,j}^{p,q}\right)\\ \text{(A8)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{m}}\sum _{j=\mathrm{1}}^m\mathit{\delta }\left(\mathbf{\Delta }{\widehat{\mathit{t}}}^{p,q}-\mathbf{\Delta }{\mathit{t}}_j^{p,q}\right)\end{array}$$

constrains the set that determines ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}_{\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}}\left(\mathbf{\Delta }{\widehat{\mathit{t}}}^{p,q}\right)$ to m values, where those values from different DO events with the same MCMC index j are combined:

$$\begin{array}{}\text{(A9)}& {\displaystyle }{\displaystyle }\mathbf{\Delta }{\mathit{t}}_j^{p,q}=\left(\mathrm{\Delta }{t}_{\mathrm{1},j}^{p,q},\mathrm{\dots },\mathrm{\Delta }{t}_{n,j}^{p,q}\right).\end{array}$$

Again, the validity is checked by randomly permuting the sets $\mathit{\left\{}\mathrm{\Delta }{t}_{i,j}^{p,q}\mathit{\right\}}$ for the individual DO events with respect to the index j before the set reduction in the control runs.

Having found a numerically manageable expression for the empirical uncertainty distribution of the sample $\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}$ it remains to be shown how the hypothesis tests can be formulated on this basis. If {Δt_j}_j denotes the set of n-dimensional vectors forming the empirical uncertainty distribution for the sample of lags obtained from n DO events, then the naive intuition holds true and the corresponding set {ϕ_j=ϕ(Δt_j)}_j represents the empirical uncertainty distribution of the test statistic and correspondingly {p_ϕ(ϕ_j)}_j characterizes the uncertain p value. In the following, we exemplarily derive this relation for the t test – the derivations for the WSR and the bootstrap test are analogous.

Recall the statistic of the t test:

$$\begin{array}{}\text{(A10)}& {\displaystyle }{\displaystyle }z\left(\mathit{d}\right)={\displaystyle \frac{u\left(\mathit{d}\right)-\mathit{\mu }}{s\left(\mathit{d}\right)/\sqrt{\left(n\right)}}}.\end{array}$$

The empirical uncertainty distribution for a sample $\mathbf{\Delta }\widehat{\mathit{T}}$ induces a joint uncertainty distribution for the sample's mean and standard deviation:

$$\begin{array}{ll}{\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\widehat{U},\widehat{S}}(\widehat{u},\widehat{s})& {\displaystyle }=\int \mathit{\delta }\left(u-{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n\mathrm{\Delta }{\widehat{t}}_i\right)\\ {\displaystyle }& {\displaystyle }\times \mathit{\delta }\left(s-{\displaystyle \frac{\mathrm{1}}{n-\mathrm{1}}}\sum _{i=\mathrm{1}}^n(u-\mathrm{\Delta }{\widehat{t}}_i{)}^{\mathrm{2}}\right)\\ {\displaystyle }& {\displaystyle }\times {\displaystyle \frac{\mathrm{1}}{m}}\sum _{j=\mathrm{1}}^m\prod _{i=\mathrm{1}}^n\mathit{\delta }(\mathrm{\Delta }{\widehat{t}}_i-\mathrm{\Delta }{t}_{i,j})\mathrm{d}\mathrm{\Delta }{\widehat{t}}_{\mathrm{1}}\mathrm{\dots }d\mathrm{\Delta }{\widehat{t}}_n\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{m}}\sum _{j=\mathrm{1}}^m\mathit{\delta }\left(\widehat{u}-{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n\mathrm{\Delta }{t}_{i,j}\right)\\ \text{(A11)}& {\displaystyle }& {\displaystyle }\times \mathit{\delta }\left(\widehat{s}-{\displaystyle \frac{\mathrm{1}}{n-\mathrm{1}}}\sum _{i=\mathrm{1}}^n(\widehat{u}-\mathrm{\Delta }{t}_{i,j}{)}^{\mathrm{2}}\right).\end{array}$$

Let $u_j=\frac{\mathrm{1}}{n}{\sum }_{i=\mathrm{1}}^n\mathrm{\Delta }{t}_{i,j}$ and $s_j=\frac{\mathrm{1}}{n-\mathrm{1}}{\sum }_{i=\mathrm{1}}^n(u_j-\mathrm{\Delta }{t}_{i,j}{)}^{\mathrm{2}}$. Then, the empirical uncertainty distribution for $(\widehat{U},\widehat{S})$ can be written as

$$\begin{array}{}\text{(A12)}& {\displaystyle }{\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\widehat{U},\widehat{S}}(\widehat{u},\widehat{s})={\displaystyle \frac{\mathrm{1}}{m}}\sum _{j=\mathrm{1}}^m\mathit{\delta }\left(\widehat{u}-u_j\right)\mathit{\delta }\left(\widehat{s}-s_j\right).\end{array}$$

The (u_j,s_j) that forms the empirical uncertainty distribution is simply the mean and standard deviation of those $\mathbf{\Delta }{\mathit{t}}_j=(\mathrm{\Delta }{t}_{\mathrm{1},j},\mathrm{\Delta }{t}_{\mathrm{2},j},\mathrm{\dots },\mathrm{\Delta }{t}_{n,j})$ that form the vector valued empirical uncertainty distribution for $\mathbf{\Delta }\widehat{\mathit{T}}$. From ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\widehat{U},\widehat{S}}(\widehat{u},\widehat{s})$, the empirical uncertainty distribution for the uncertain test statistic $\widehat{Z}$ can be computed as follows:

$$\begin{array}{ll}{\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\widehat{Z}}\left(\widehat{z}\right)& {\displaystyle }=\int \mathit{\delta }\left(\widehat{z}-{\displaystyle \frac{\widehat{u}-\mathit{\mu }}{\widehat{s}/\sqrt{\left(n\right)}}}\right){\mathit{\rho }}_{\widehat{U},\widehat{S}}(\widehat{u},\widehat{s})\mathrm{d}\widehat{u}\mathrm{d}\widehat{s}\\ \text{(A13)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{m}}\sum _{j=\mathrm{1}}^m\mathit{\delta }(\widehat{z}-\underset{=z_j}{\underbrace{{\displaystyle \frac{u_j-\mathit{\mu }}{s_j/\sqrt{\left(n\right)}}}}}).\end{array}$$

This shows that for a given empirical uncertainty distribution for a sample of time lags ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathbf{\Delta }\widehat{\mathit{T}}}\left(\mathbf{\Delta }\widehat{\mathit{t}}\right)=\frac{\mathrm{1}}{m}{\sum }_{i=\mathrm{1}}^m\mathit{\delta }\left(\mathbf{\Delta }\widehat{\mathit{t}}-\mathbf{\Delta }{\mathit{t}}_j\right)$, the corresponding distribution for the test statistic $\widehat{Z}=z\left(\mathbf{\Delta }\widehat{\mathit{T}}\right)$ is formed by the set $\mathit{\left\{}z\right(\mathbf{\Delta }{\mathit{t}}_j\left)\right|j\in [\mathrm{1},m]\mathit{\}}$ where each Δt_j is a vector in n dimensions. The uncertain (left-handed) p value remains to be derived from ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\widehat{Z}}\left(\widehat{z}\right)$:

$$\begin{array}{ll}{\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{{\widehat{P}}_z}\left({\widehat{p}}_z\right)& {\displaystyle }=\int \mathit{\delta }\left({\widehat{p}}_z-\underset{-\mathrm{\infty }}{\overset{\widehat{z}}{\int }}{t}_{n-\mathrm{1}}\left(z\right)\mathrm{d}z\right){\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\widehat{Z}}\left(\widehat{z}\right)\mathrm{d}\widehat{z}\\ \text{(A14)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{m}}\sum _{j=\mathrm{1}}^m\mathit{\delta }({\widehat{p}}_z-\underset{=p_{z,j}}{\underbrace{\underset{-\mathrm{\infty }}{\overset{z_j}{\int }}{t}_{n-\mathrm{1}}\left(z^{\prime }\right)}}\mathrm{d}{z}^{\prime }).\end{array}$$

Finally, the practical computation of the uncertain p values boils down to an application of the test to all members of the set Δt_j that originates from the MCMC sampling used to approximate the posterior probability density for the ramp parameter configuration Θ. For the WSR test the expression

$$\begin{array}{}\text{(A15)}& {\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{{\widehat{P}}_w}\left({\widehat{p}}_w\right)={\displaystyle \frac{\mathrm{1}}{m}}\sum _{j=\mathrm{1}}^m\mathit{\delta }\left({\widehat{p}}_w-p_{w,j}\right)\end{array}$$

with

$$p_{w,j}=p_w\left(\mathbf{\Delta }{\mathit{t}}_j\right)$$

can be derived analogously. The bootstrap test bears the particularity that each Δt_j induces its own null distribution. Yet, the application of the test to each individual Δt_j induces a set of $p_{v,j}=p_v\left(\mathbf{\Delta }{\mathit{t}}_j\right)$ that determines the empirical density:

$$\begin{array}{}\text{(A16)}& {\displaystyle }{\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{{\widehat{P}}_v}\left({\widehat{p}}_v\right)={\displaystyle \frac{\mathrm{1}}{m}}\sum _{j=\mathrm{1}}^m\mathit{\delta }\left({\widehat{p}}_v-p_{v,j}\right).\end{array}$$

Appendix B: Results of the analysis for the control group

As explained in Sect. A, we drastically reduce the cardinality of the sets that form the empirical densities ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}}\left(\mathbf{\Delta }{\widehat{\mathit{t}}}^{p,q}\right)$ at two points in the analysis. First, for the representation of the uncertain time lag $\mathrm{\Delta }{\widehat{T}}_i^{p,q}$ between the proxies p and q at a given DO event, only 6000 out of the 6000² possible values are utilized. Second, the set of vectors considered in the representation of ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}}\left(\mathbf{\Delta }\widehat{\mathit{t}}\right)=\frac{\mathrm{1}}{\mathrm{6000}}{\sum }_{j=\mathrm{1}}^{\mathrm{6000}}\mathit{\delta }(\mathbf{\Delta }{\widehat{\mathit{t}}}^{p,q}-\mathbf{\Delta }{\mathit{t}}_j^{p,q})$ is comprised of only 6000 out of the 6000¹⁶ theoretically available vectors. To cross-check the robustness of the results obtained within the limits of this approximation, we applied our analysis to a control group of nine alternative realizations of the empirical uncertainty density for $\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}$ for each proxy pair. The control group uncertainty densities are constructed as follows: first, the empirical uncertainty distributions for the event-specific lags $\mathrm{\Delta }{\widehat{T}}_i^{p,q}$ are obtained via Eq. (A6). In a second step, the joint empirical uncertainty distribution for $\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}$ is constructed from randomly shuffled empirical sets $\mathrm{\Delta }{t}_{i,s_i\left(j\right)}^{p,q}$ of each DO event:

$$\begin{array}{}\text{(B1)}& {\displaystyle }{\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_{\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}}^{\text{ctrl}}\left(\mathbf{\Delta }{\widehat{\mathit{t}}}^{p,q}\right)={\displaystyle \frac{\mathrm{1}}{m}}\sum _{j=\mathrm{1}}^m\prod _{i=\mathrm{1}}^n\mathit{\delta }\left(\mathrm{\Delta }{\widehat{t}}_i^{p,q}-\mathrm{\Delta }{t}_{i,s_i\left(j\right)}^{p,q}\right).\end{array}$$

Here s_i denotes an event-specific permutation of the index set $\mathit{\{}\mathrm{1},\mathrm{\dots },\mathrm{6000}\mathit{\}}$. Thus the empirical $\mathrm{\Delta }{t}_{i,j}^{p,q}$ recombines between events and gives rise to a new set of 6000 vectors that constitute 6000 empirical realizations of the uncertain $\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}$.

Table B1 Results obtained from the application of hypothesis tests to the control group. Reported are the mean p values $\mathrm{E}\left(p\right(\mathbf{\Delta }\widehat{\mathit{T}}\left)\right)$ together with the probability of the uncertain sample to yield a p value below the significance level $\mathit{\pi }\left(p\right(\mathbf{\Delta }\widehat{\mathit{T}})<\mathrm{0.05})$ and the p values of the expected samples $p\left(\mathrm{E}\right(\mathbf{\Delta }\widehat{\mathit{T}}\left)\right)$ for all three tests. All results were derived from the corresponding empirical densities ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}}\left(\mathbf{\Delta }{\widehat{\mathit{t}}}^{p,q}\right)$. The column sub-labels z, w, and bs indicate results obtained from the t test, the WSR test, and the bootstrap test. The results presented in the main text are given by the p–q–0 run for each proxy variable. The data presented here are provided in the Supplement to this article.

The results obtained from the control runs show only minor deviations from the results presented in the main text and thus confirm the validity of the reduction in the corresponding sets. Table B1 summarizes the results obtained by the application of the hypothesis tests to the control group.

Appendix C: Computation of the uncertain sample mean

In the main text, we stated that the uncertain sample mean is given by the pairwise convolution of the individual uncertainty distributions that describe the uncertain sample members. Here, we show how the uncertain sample mean can be computed if the individual uncertainty distributions are known.

Consider n random variables which are independently yet not identically distributed:

$$\begin{array}{}\text{(C1)}& {\displaystyle }{\displaystyle }\mathit{X}=(X_{\mathrm{1}},\mathrm{\dots },X_n)\text{ with }{X}_i\sim {\mathit{\rho }}_{X_i}\left(x_i\right)\mathrm{d}{x}_i\end{array}$$

in analogy to the

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathbf{\Delta }{\widehat{\mathit{T}}}^{p,q}=\left(\mathrm{\Delta }{\widehat{T}}_{\mathrm{1}}^{p,q},\mathrm{\dots },\mathrm{\Delta }{\widehat{T}}_n^{p,q}\right)\\ \text{(C2)}& {\displaystyle }& {\displaystyle }\text{with }\mathrm{\Delta }{\widehat{T}}_i^{p,q}\sim {\mathit{\rho }}_{\mathrm{\Delta }{\widehat{T}}_i^{p,q}}\left(\mathrm{\Delta }{\widehat{t}}_i^{p,q}\right)\end{array}$$

from the main text. Further, let

$$\begin{array}{}\text{(C3)}& {\displaystyle }{\displaystyle }U={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{X}_i\end{array}$$

denote the mean of the sample of random variables, which is in turn a random variable by itself. In order to compute the distribution ρ_U(u) du we introduce the variable V=nU and the sequence of variables

$$\begin{array}{}\text{(C4)}& {\displaystyle }{\displaystyle }{V}_j=\sum _{i=\mathrm{1}}^j{X}_i,\end{array}$$

such that V_n=V. From Eq. (C4) it follows that

$$\begin{array}{}\text{(C5)}& {\displaystyle }{\displaystyle }{V}_{j+\mathrm{1}}=V_j+X_{j+\mathrm{1}}\end{array}$$

and hence

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\rho }}_{V_{j+\mathrm{1}}}\left(v_{j+\mathrm{1}}\right)\mathrm{d}{v}_{j+\mathrm{1}}\\ {\displaystyle }& {\displaystyle }=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\mathit{\rho }}_{V_j}\left(v_j\right){\mathit{\rho }}_{X_{j+\mathrm{1}}}\left(x_{j+\mathrm{1}}\right)\\ {\displaystyle }& {\displaystyle }\times \mathit{\delta }(v_{j+\mathrm{1}}-v_j-x_{j+\mathrm{1}})\mathrm{d}{x}_{j+\mathrm{1}}\mathrm{d}{v}_j\mathrm{d}{v}_{j+\mathrm{1}}\\ \text{(C6)}& {\displaystyle }& {\displaystyle }=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\mathit{\rho }}_{V_j}\left(v_j\right){\mathit{\rho }}_{X_{j+\mathrm{1}}}(v_{j+\mathrm{1}}-v_j)\mathrm{d}{v}_j\mathrm{d}{v}_{j+\mathrm{1}}.\end{array}$$

Self-iteration of Eq. (C6) yields

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\rho }}_{V_{j+\mathrm{1}}}\left(v_{j+\mathrm{1}}\right)\mathrm{d}{v}_{j+\mathrm{1}}\\ {\displaystyle }& {\displaystyle }=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\underset{={\mathit{\rho }}_{V_j}\left(v_j\right)}{\underbrace{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\mathit{\rho }}_{V_{j-\mathrm{1}}}\left(v_{j-\mathrm{1}}\right){\mathit{\rho }}_{X_j}(v_j-v_{j-\mathrm{1}})\mathrm{d}{v}_{j-\mathrm{1}}}}\\ {\displaystyle }& {\displaystyle }\times {\mathit{\rho }}_{X_{j+\mathrm{1}}}(v_{j+\mathrm{1}}-v_j)\mathrm{d}{v}_j\mathrm{d}{v}_{j+\mathrm{1}}\\ {\displaystyle }& {\displaystyle }=\mathrm{\dots }\\ \text{(C7)}& {\displaystyle }& {\displaystyle }=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\mathrm{\dots }\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\prod _{i=\mathrm{1}}^{j+\mathrm{1}}{\mathit{\rho }}_{X_i}(v_i-v_{i-\mathrm{1}})\mathrm{d}{v}_{i-\mathrm{1}}\mathrm{d}{v}_{j+\mathrm{1}},\end{array}$$

where v₀=0. With $V_n/n=U$, the distribution for the uncertain sample mean reads

$$\begin{array}{}\text{(C8)}& {\displaystyle }{\mathit{\rho }}_{V_n}\left(v_n\right)\mathrm{d}{V}_n={\mathit{\rho }}_{V_n}\left(nu\right)n\mathrm{d}u={\mathit{\rho }}_U\left(u\right)\mathrm{d}u,\end{array}$$

and thus

$$\begin{array}{}\text{(C9)}& {\displaystyle }{\mathit{\rho }}_U\left(u\right)\mathrm{d}u=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\prod _{i=\mathrm{1}}^n{\mathit{\rho }}_{X_i}(v_i-v_{i-\mathrm{1}})\mathrm{d}{v}_{i-\mathrm{1}}n\mathrm{d}u,\end{array}$$

with v₀=0 and v_n=nu.