Modular and Adaptive Conformal Prediction for Sequential Models via Residual Decomposition

William Zhang
wzhang64@mit.edu
Operations Research Center, Massachusetts Institute of Technology, Cambridge, MA 02139
   Saurabh Amin11footnotemark: 1
amins@mit.edu
   Georgia Perakis11footnotemark: 1
georgiap@mit.edu
Abstract

Conformal prediction offers finite-sample coverage guarantees under minimal assumptions. However, existing methods treat the entire modeling process as a black box, overlooking opportunities to exploit modular structure. We introduce a conformal prediction framework for two-stage sequential models, where an upstream predictor generates intermediate representations for a downstream model. By decomposing the overall prediction residual into stage-specific components, our method enables practitioners to attribute uncertainty to specific pipeline stages. We develop a risk-controlled parameter selection procedure using family-wise error rate (FWER) control to calibrate stage-wise scaling parameters, and propose an adaptive extension for non-stationary settings that preserves long-run coverage guarantees. Experiments on synthetic distribution shifts, as well as real-world supply chain and stock market data, demonstrate that our approach maintains coverage under conditions that degrade standard conformal methods, while providing interpretable stage-wise uncertainty attribution. This framework offers diagnostic advantages and robust coverage that standard conformal methods lack.

1 Introduction

Modern machine learning systems increasingly rely on modular pipelines, where upstream predictors generate intermediate representations for downstream models. Such pipelines appear across diverse applications: macroeconomic forecasting, where supply chain indicators inform market predictions, and medical diagnosis, where imaging features guide treatment decisions [11, 25]. In these settings, uncertainty quantification is essential and should reflect how error propagates across stages. Conformal prediction [28] provides principled prediction intervals with finite-sample coverage guarantees under exchangeability. However, existing conformal prediction methods treat these multi-stage systems as monolithic black boxes, overlooking their modular structure and missing opportunities for targeted error attribution.

Recent work expands the flexibility of conformal prediction through reweighting and adaptive calibration [7, 1, 20, 12, 4, 14, 15], structured model-aware adaptations [31, 32, 29], and risk-control [8, 2], but largely focuses on single-stage or black-box models, rather than taking a modular perspective.

This perspective is especially pertinent under distribution shift, which often affects stages asymmetrically—e.g., upstream sensors may drift while downstream mappings remain stable. Standard conformal methods cannot disentangle such effects, which may force practitioners to retrain the entire pipeline when targeted interventions would suffice. To address this gap, we propose a stage-wise abstraction that isolates these effects, embedding stage-wise uncertainty into conformal prediction, yielding robust, interpretable intervals under distribution shift.

By decomposing residuals into stage-wise components, our method constructs prediction intervals through a linear combination of stage-specific quantiles, weighted by scaling parameters selected via family-wise error rate (FWER) control over a calibration set. This yields valid, interpretable intervals without requiring internal model access—only structural knowledge of the pipeline. Our approach offers two key advantages over existing conformal methods: (i) it provides diagnostic transparency, identifying which stages contribute most to uncertainty, and (ii) it maintains coverage under distribution shifts affecting individual components where standard approaches degrade.

For intuition, we briefly introduce the two-stage setting with triplets (w,x,y)(w,x,y) and latent structure x=μ1(w)+ε1x=\mu_{1}(w)+\varepsilon_{1}, y=μ2(x)+ε2y=\mu_{2}(x)+\varepsilon_{2}, where ε1,ε2\varepsilon_{1},\varepsilon_{2} denote additive noise, and the model learns estimators μ^1\hat{\mu}_{1}, μ^2\hat{\mu}_{2}. For example, in automobile supply chains, ww could denote semiconductor prices, with μ^1(w)\hat{\mu}_{1}(w) estimating new vehicle demand (xx), and μ^2(x)\hat{\mu}_{2}(x) predicting used vehicle prices (yy). Our method decomposes the residual R=|yμ^2(μ^1(w))|R=|y-\hat{\mu}_{2}(\hat{\mu}_{1}(w))| into upstream (μ^1\hat{\mu}_{1}) error component ΔR1\Delta R_{1} and downstream (μ^2\hat{\mu}_{2}) residual R2R_{2}, capturing the uncertainty of each stage.

We also extend our framework to adaptive settings for which we update scaling parameters based on component-wise empirical coverage, improving responsiveness to (i) upstream, (ii) downstream, and (iii) end-to-end distribution shifts. We illustrate our method on synthetic shifts and real-world supply chain and financial data, showing improved robustness and interpretability over adaptive conformal baselines [14, 15, 1, 4]. For simplicity, we focus on two-stage models, but our methodology can be extended to multi-stage models (Appendix C).

Contributions:
  • We propose a residual decomposition framework for sequential multi-stage models, that partitions prediction error into interpretable upstream and downstream components, enabling stage-wise uncertainty attribution (Section˜3).

  • We develop a risk-controlled parameter selection procedure using FWER-based hypothesis testing to construct valid prediction intervals from decomposed residuals, with formal coverage guarantees (Section˜5).

  • We introduce an adaptive algorithm that dynamically adjusts scaling parameters based on component-wise coverage feedback, preserving long-run coverage guarantees while providing interpretable diagnostics for distribution shifts (Section˜6).

  • We demonstrate the framework’s effectiveness on synthetic shifts and real-world economic forecasting, showing maintained coverage under conditions that degrade existing conformal methods, with diagnostic capabilities that enable targeted model interventions (Section˜7).

Paper Outline. In Sections 2 and 3, we formalize the two-stage prediction setting and introduce our residual decomposition. Section˜4 describes our method for constructing stage-aware prediction intervals, followed by the FWER calibration procedure in Section˜5. We extend our method to an adaptive version Section˜6, and evaluate performance under distribution shift in Section˜7.

2 Problem setting and assumptions

We consider a sequential two-stage prediction problem. Each data point is characterized by a triplet z=(w,x,y)z=(w,x,y) where w𝒲w\in\mathcal{W} is the input to the first-stage model (upstream features), x𝒳x\in\mathcal{X} is an intermediate representation, and y𝒴y\in\mathcal{Y} is the final prediction target. We learn predictors μ^1:𝒲𝒳\hat{\mu}_{1}:\mathcal{W}\to\mathcal{X} and μ^2:𝒳𝒴\hat{\mu}_{2}:\mathcal{X}\to\mathcal{Y} which are composed to form end-to-end predictor μ^2(μ^1(w))=μ^2(x^)\hat{\mu}_{2}(\hat{\mu}_{1}(w))=\hat{\mu}_{2}(\hat{x}), where x^=μ^1(w)\hat{x}=\hat{\mu}_{1}(w).

To learn these models and perform conformal prediction, we assume access to three disjoint subsets: (i) a training set Strain={(wi,xi,yi)}i=1nS^{\text{train}}=\{(w_{i},x_{i},y_{i})\}_{i=1}^{n} used to fit both stages of the model via (w,x)(w,x) and (x,y)(x,y) pairs for each stage respectively; (ii) a conformal set Sconf={(wi,xi,yi)}i=n+1n+mS^{\text{conf}}=\{(w_{i},x_{i},y_{i})\}_{i=n+1}^{n+m} for computing nonconformity scores; and (iii) a calibration set Scal={(wi,xi,yi)}i=n+m+1n+m+lS^{\text{cal}}=\{(w_{i},x_{i},y_{i})\}_{i=n+m+1}^{n+m+l} for parameter selection. At prediction time, only the upstream input wtestw_{\text{test}} is observed, while the intermediate value xtestx_{\text{test}} and target ytesty_{\text{test}} are unobserved and must be predicted.

We list some assumptions. Exchangeability. Unless stated otherwise, we assume the data (w,x,y)(w,x,y) in Strain,Sconf,ScalS^{\text{train}},S^{\text{conf}},S^{\text{cal}} to be exchangeable, as well as the test point ztest=(wtest,xtest,ytest)z_{\text{test}}=(w_{\text{test}},x_{\text{test}},y_{\text{test}}). Learning algorithms. We assume that the algorithms that learn μ^1,μ^2\hat{\mu}_{1},\hat{\mu}_{2} are deterministic, i.e. given the same training data, they produce identical predictors; and we also assume they are symmetric, i.e. invariant to permutations of the data. We also assume that the intermediate variable xx is observable for the given datasets, enabling residual decomposition. Distribution shifts. Distribution shifts violate exchangeability and can impact each prediction stage differently. Let PP denote the distribution of a data point in the training set and PP^{\prime} denote the distribution of the test point. We consider three types of shift: Upstream covariate shift: P(w)P(w)P(w)\neq P^{\prime}(w); Upstream concept shift: P(x|w)P(x|w)P(x|w)\neq P^{\prime}(x|w); and Downstream concept shift: P(y|x)P(y|x)P(y|x)\neq P^{\prime}(y|x).

Objective.

Given test upstream input wtestw_{\text{test}}, the goal is to construct prediction interval C^α(wtest)\hat{C}_{\alpha}(w_{\text{test}}) such that (ytestC^α(wtest))1α\mathbb{P}(y_{\text{test}}\in\hat{C}_{\alpha}(w_{\text{test}}))\geq 1-\alpha, where α(0,1)\alpha\in(0,1) is the target miscoverage level. Under exchangeability, this is the standard conformal prediction objective. However, under distribution shifts, we seek to maintain robust coverage while providing interpretable attribution of uncertainty to specific pipeline stages.

The key challenges in this setting are: (i) Attribution: Understanding which stage contributes more to prediction uncertainty; (ii) Adaptivity: Maintaining coverage under shifts affecting different stages; (iii) Transparency: Providing actionable insights for model improvement or retraining decisions. Our approach addresses these challenges by decomposing the end-to-end prediction error into stage-specific components for targeted uncertainty quantification and robust interval construction. While we define these concepts for two-stage models, we discuss auxiliary inputs, multiple upstream models, and deeper sequential pipelines in Appendix C. We also provide a notation table in Appendix B.2.

3 Two-stage residual decomposition for conformal prediction

To address the aforementioned challenge of attribution, we partition the total prediction residual R(w,y)=|yμ^2(μ^1(w))|R(w,y)=|y-\hat{\mu}_{2}(\hat{\mu}_{1}(w))| into upstream and downstream components. This decomposition enables stage-wise attribution of error, in contrast to standard black-box conformal methods. We provide a visualization in Figure˜1.

Refer to caption
Figure 1: Visualization of a sequential two-stage model with residual components R,ΔR1,R2R,\Delta R_{1},R_{2}.
Definition 1 (Second-stage residual).

Given a point (w,x,y)(w,x,y) and downstream predictor μ^2:𝒳𝒴\hat{\mu}_{2}:\mathcal{X}\to\mathcal{Y}, the second-stage residual is

R2(x,y)=|yμ^2(x)|.R_{2}(x,y)=|y-\hat{\mu}_{2}(x)|.

This captures downstream prediction error when given the true intermediate input xx.

Definition 2 (First-stage delta).

Let x^=μ^1(w)\hat{x}=\hat{\mu}_{1}(w) be the output of the first-stage predictor. The first-stage delta is defined as

ΔR1(w,x,y)=||yμ^2(x)||yμ^2(x^)||.\Delta R_{1}(w,x,y)=\left|\,|y-\hat{\mu}_{2}(x)|-|y-\hat{\mu}_{2}(\hat{x})|\,\right|.

This quantifies the change in downstream prediction error induced by replacing the true intermediate xx with its prediction x^\hat{x}.

For notational purposes, we drop the dependence on (w,x,y)(w,x,y), and refer to these terms as R,ΔR1R,\Delta R_{1}, and R2R_{2}. Intuitively, ΔR1\Delta R_{1} reflects the error of the upstream predictor, while R2R_{2} isolates downstream error without upstream influence. This decomposition satisfies a fundamental upper-bound property.

Proposition 1.

For any point (w,x,y)(w,x,y), the total residual satisfies:

R=|yμ^2(μ^1(w))|ΔR1+R2.R=|y-\hat{\mu}_{2}(\hat{\mu}_{1}(w))|\leq\Delta R_{1}+R_{2}.

By definition, ΔR1=||yμ^2(x)||yμ^2(x^)||\Delta R_{1}=||y-\hat{\mu}_{2}(x)|-|y-\hat{\mu}_{2}(\hat{x})||. Let A=|yμ^2(x)|=R2A=|y-\hat{\mu}_{2}(x)|=R_{2} and B=|yμ^2(x^)|=RB=|y-\hat{\mu}_{2}(\hat{x})|=R. By the reverse triangle inequality, BA+|AB|=R2+ΔR1B\leq A+|A-B|=R_{2}+\Delta R_{1}. This upper bound property ensures that the sum of components provides a conservative estimate of the total error, which is crucial for the coverage guarantees developed subsequently.

This decomposition provides a clean interpretation: R2R_{2} measures the inherent uncertainty of the downstream model, while ΔR1\Delta R_{1} measures how upstream prediction errors affect the final prediction error magnitude. When ΔR1\Delta R_{1} is small relative to R2R_{2}, the upstream predictor performs well and downstream uncertainty dominates. Conversely, when ΔR1\Delta R_{1} is large relative to R2R_{2}, upstream errors drive prediction uncertainty. This attribution enables practitioners to identify which stage requires improvement and guide retraining decisions. Furthermore, these components provide insights for handling distribution shifts: under upstream shifts, ΔR1\Delta R_{1} typically increases as the first-stage predictor encounters out-of-distribution inputs, while under downstream shifts, R2R_{2} increases. These changes can be experimentally visualized in Appendix Figure˜12(b). Thus, the varied responses enable targeted adaptive strategies for different distribution shifts.

These components exhibit intuitive relationships with the total residual: when R2R_{2} is small, ΔR1\Delta R_{1} closely approximates RR, while when the upstream model is accurate and μ^2\hat{\mu}_{2} is smooth, R2R_{2} closely approximates RR (see Appendix A.2). Importantly, at least one component must represent a majority of the error, ensuring meaningful stage-wise attribution. Next, we describe two complementary approaches to combining these residual components into prediction intervals.

4 Constructing prediction intervals

We describe two approaches that incorporate ΔR1,R2\Delta R_{1},R_{2}, utilizing the conformal set SconfS^{\text{conf}}, to compute component-wise quantiles, with data-driven parameter selection using ScalS^{\text{cal}} addressed in Section˜5. Note that there exists a rich space of heuristics for combining these residual components beyond those listed below—see Appendix C.1. For proofs of the following results, see Appendix A.1.

4.1 Separate component quantiles

We compute sets of residual components on SconfS^{\text{conf}}: {ΔR1i}i=1m\{\Delta R_{1}^{i}\}_{i=1}^{m} and {R2i}i=1m\{R_{2}^{i}\}_{i=1}^{m}, which we abbreviate as {ΔR1}\{\Delta R_{1}\} and {R2}\{R_{2}\}. Prediction intervals are constructed by summing quantiles computed separately from each set, with the quantile levels controlled by cc and dd for the respective stages.

Definition 3 (Separate component quantiles).

Let c,d(0,1)c,d\in(0,1) be quantile levels. For test input wtestw_{\text{test}}, the prediction interval is

C^c,d(wtest)=μ^2(μ^1(wtest))±(Q1c({ΔR1})+Q1d({R2})).\hat{C}_{c,d}(w_{\text{test}})=\hat{\mu}_{2}(\hat{\mu}_{1}(w_{\text{test}}))\pm(Q_{1-c}(\{\Delta R_{1}\})+Q_{1-d}(\{R_{2}\})).

This construction inherits coverage guarantees from standard conformal prediction:

Theorem 1 (Coverage of separate component quantiles).

Under the assumption of exchangeability, for c,d(0,1)c,d\in(0,1), the prediction interval C^c,d(wtest)\hat{C}_{c,d}(w_{\text{test}}) satisfies

(ytestC^c,d(wtest))1cd.\mathbb{P}(y_{\text{test}}\in\hat{C}_{c,d}(w_{\text{test}}))\geq 1-c-d.

4.2 Scaled component quantiles

Our second approach fixes a quantile level α(0,1)\alpha\in(0,1) for both components and selects scaling parameters a,b[0,1]a,b\in[0,1] to weight their respective contributions. This provides interpretable control over stage-wise uncertainty attribution.

Definition 4 (Scaled component quantiles).

For a fixed quantile α(0,1)\alpha\in(0,1) and scaling coefficients a,b[0,1]a,b\in[0,1], the prediction interval is

C^α,a,b(wtest)=μ^2(μ^1(wtest))±(aQ1α({ΔR1})+bQ1α({R2}))\hat{C}_{\alpha,a,b}(w_{\text{test}})=\hat{\mu}_{2}(\hat{\mu}_{1}(w_{\text{test}}))\pm(a\cdot Q_{1-\alpha}(\{\Delta R_{1}\})+b\cdot Q_{1-\alpha}(\{R_{2}\}))

where the quantiles are computed over SconfS^{\text{conf}}.

The choice of scaling parameters aa and bb provides interpretable control: setting a=0a=0 ignores upstream uncertainty while b=0b=0 focuses only on upstream effects. However, coverage guarantees for arbitrary choices of (a,b)(a,b) require careful analysis. For fixed weights, we have the following:

Corollary 1 (Coverage with a=b=1a=b=1).

Under exchangeability, for a=b=1a=b=1 and α(0,1)\alpha\in(0,1), the interval C^α,a,b(wtest)\hat{C}_{\alpha,a,b}(w_{\text{test}}) satisfies

(ytestC^α,a,b(wtest))12α.\mathbb{P}(y_{\text{test}}\in\hat{C}_{\alpha,a,b}(w_{\text{test}}))\geq 1-2\alpha.
Proof.

This follows directly from Theorem˜1 with c=d=αc=d=\alpha. ∎

For appropriately chosen quantiles and scaling weights, we observe that both can methods yield similar intervals. To provide maximum flexibility for both theoretical analysis and implementation, we can combine both approaches into a general framework:

C^a,b,c,d(wtest)=μ^2(μ^1(wtest))±(aQ1c({ΔR1})+bQ1d({R2})),\hat{C}_{a,b,c,d}(w_{\text{test}})=\hat{\mu}_{2}(\hat{\mu}_{1}(w_{\text{test}}))\pm(a\cdot Q_{1-c}(\{\Delta R_{1}\})+b\cdot Q_{1-d}(\{R_{2}\})),

This unified form allows independent control of both quantile levels and scaling parameters, enabling fine-tuned balance between coverage guarantees and stage-wise identifiability.

4.3 Coverage for scaled parameters

While coverage for intervals of the form C^c,d\hat{C}_{c,d} (Definition 3) follows from standard conformal analysis, establishing similar guarantees for intervals with scaled residuals C^α,a,b\hat{C}_{\alpha,a,b} (Definition 4) presents significant challenges. The scaling parameters (a,b)(a,b) have no direct mapping to quantile levels, making it difficult to derive explicit guarantees. Despite this, we can establish that valid scaling parameters exist in principle. Under mild regularity conditions, there always exist optimal scaling parameters (a,b)(a^{*},b^{*}) that yield exact marginal coverage:

Proposition 2 (Existence of ideal scaling parameters).

For desired coverage level 1α1-\alpha, a,b[0,1]\exists a^{*},b^{*}\in[0,1] such that the interval with those scaling parameters satisfies the marginal coverage guarantee

(ytestC^α,a,b(wtest))=1α,\mathbb{P}\left(y_{\text{test}}\in\hat{C}_{\alpha,a^{*},b^{*}}(w_{\text{test}})\right)=1-\alpha,

provided the distribution of the residual RR has no point masses.

However, finding optimal scaling parameters (a,b)(a^{*},b^{*}) requires knowledge of the residual distribution, which is unavailable in practice. This creates a fundamental trade-off: scaled intervals offer interpretable stage-wise control but lack accessible guarantees, while separate quantile intervals provide coverage under minimal assumptions but offer less direct interpretability. To resolve this trade-off, we adopt a data-driven approach [8, 3, 2] to selecting scaling parameters (a,b)(a,b) from a candidate set using the calibration set ScalS^{\text{cal}} for the unified interval.

5 Risk-controlling approach with residual components

Since finding (a,b)(a^{*},b^{*}) is impossible in practice, we reframe the problem as selecting (a,b)(a,b) coefficient pairs that satisfy the nominal coverage level 1α1-\alpha. We approach this from the perspective of multiple hypothesis testing; for each (a,b)(a,b) we test the null hypothesis that the miscoverage rate exceeds α\alpha.

5.1 Testing miscoverage via empirical risk

We aim to identify scaling parameters λ=(a,b)[0,1]2\lambda=(a,b)\in[0,1]^{2} for which the resulting prediction interval C^λ(w)\hat{C}_{\lambda}(w) achieves coverage at least 1α1-\alpha. Since theoretical guarantees for arbitrary λ\lambda are unavailable, we perform a hypothesis test for whether its miscoverage rate exceeds α\alpha. We fix a finite candidate set Λ[0,1]2\Lambda\subseteq[0,1]^{2} of scaling pairs and define the corresponding prediction interval for each λ=(a,b)Λ\lambda=(a,b)\in\Lambda:

C^λ(w)={y:|yμ^2(μ^1(w))|aQ1c({ΔR1})+bQ1d({R2})},\hat{C}_{\lambda}(w)=\left\{y:\left|y-\hat{\mu}_{2}(\hat{\mu}_{1}(w))\right|\leq a\cdot Q_{1-c}(\{\Delta R_{1}\})+b\cdot Q_{1-d}(\{R_{2}\})\right\},

using the general framework introduced earlier that combines Definition˜4 and Definition˜3. Here, the quantile parameters cc and dd are fixed in advance with quantiles taken over SconfS^{\text{conf}}, while the miscoverage testing is performed using ScalS^{\text{cal}} to identify suitable λ\lambda.

5.2 Computing pp-values from calibration data

Let l=|Scal|l=|S^{\text{cal}}| denote the size of the calibration set. For each choice of scaling parameters λΛ\lambda\in\Lambda, we define the empirical miscoverage rate ^(λ)=1li=1l𝟙{yiC^λ(wi)}\hat{\mathcal{R}}(\lambda)=\frac{1}{l}\sum_{i=1}^{l}\mathbbm{1}_{\{y_{i}\notin\hat{C}_{\lambda}(w_{i})\}}. Under the stronger assumption that the calibration points are IID, and that the true miscoverage rate for a given λ\lambda is constant, the number of missed points follows a Binomial distribution with parameters ll and the underlying miscoverage probability. For each λ\lambda, define the null hypothesis

0(λ):(yC^λ(w))>α.\mathcal{H}_{0}(\lambda):\mathbb{P}(y\notin\hat{C}_{\lambda}(w))>\alpha.

Thus, for each λ\lambda we define a pp-value pλ=(Bin(l,α)l^(λ))p_{\lambda}=\mathbb{P}\left(\text{Bin}(l,\alpha)\leq l\hat{\mathcal{R}}(\lambda)\right). The pp-values pλp_{\lambda} are super-uniform under 0\mathcal{H}_{0} (Appendix A.2); that is, for any u[0,1]u\in[0,1], we have (pλu)u\mathbb{P}(p_{\lambda}\leq u)\leq u, which is crucial for FWER-controlling guarantees [8]. Thus, we apply FWER-controlling multiple testing algorithms (Appendix B.4) to the collection of pλp_{\lambda} to obtain the set of valid scaling parameters ΛvalΛ\Lambda_{\text{val}}\subseteq\Lambda. This ensures that, with probability at least 1δ1-\delta for some δ>0\delta>0, no λ\lambda with a miscoverage rate exceeding α\alpha is accepted. Note that the λΛval\lambda\in\Lambda_{\text{val}} are more conservative as they require evidence of coverage of at least 1α1-\alpha. In contrast, other prediction interval methods yield coverage that is merely close to 1α1-\alpha. This conservatism may result in unnecessarily wide intervals, particularly when coverage is less critical than efficiency, such as in IID settings. To remedy this, we can include tolerance parameter τ\tau and calculate pλp_{\lambda} with Bin(l,α+τ)\text{Bin}(l,\alpha+\tau), trading some guarantees for practical performance. As demonstrated in our experiments (Appendix D.2.2), even small values of τ\tau significantly improve efficiency while maintaining empirical coverage under IID settings.

5.3 Coverage guarantees via risk-controlled scaling

Thus, given the pp-values pλp_{\lambda} and a FWER-controlling multiple testing algorithm (Appendix B.4), we identify the set of validated scaling parameters ΛvalΛ\Lambda_{\text{val}}\subseteq\Lambda. The following result states any intervals associated with Λval\Lambda_{\text{val}} achieve coverage with high probability:

Theorem 2 (Risk control via FWER calibration).

Let ΛvalΛ\Lambda_{\text{val}}\subseteq\Lambda be the set selected by a FWER-controlling algorithm at level δ\delta, based on pp-values computed over the calibration set ScalS^{\text{cal}} with tolerance τ>0\tau>0. Then, for any λ^Λval\hat{\lambda}\in\Lambda_{\text{val}}, we have

((ytestC^λ^(wtest)|Scal)1ατ)1δ,\mathbb{P}\left(\mathbb{P}\left(y_{\text{test}}\in\hat{C}_{\hat{\lambda}}(w_{\text{test}})\,\big|\,S^{\text{cal}}\right)\geq 1-\alpha-\tau\right)\geq 1-\delta,

where the outer probability is over the randomness of ScalS^{\text{cal}}, and the inner probability is over the test point (wtest,xtest,ytest)(w_{\text{test}},x_{\text{test}},y_{\text{test}}).

This guarantee follows from the FWER-controlling property: by reducing the probability of false positives, i.e. selecting scaling parameters with poor coverage, we ensure that all selected parameters satisfy the coverage requirement with high probability. Thus, with probability at least 1δ1-\delta over ScalS^{\text{cal}}, any selected interval C^λ^\hat{C}_{\hat{\lambda}} has coverage rate at least 1ατ1-\alpha-\tau. Note that while ll does not explicitly appear in the guarantee, larger calibration sets lead to more precise pp-value estimates, typically resulting in a larger Λval\Lambda_{\text{val}} and less conservative interval selection. In practice, any pair (a,b)Λval(a,b)\in\Lambda_{\text{val}} can be used to construct prediction intervals. Although the quantile levels c,dc,d from Definition˜3 could also be tuned jointly with (a,b)(a,b), we fix c,dc,d for simplicity, which is typically sufficient in our experiments.

We note that Λval\Lambda_{\text{val}} can be empty, producing no interval. This is intentional, as it can indicate that the shift is too large and retraining is needed, rather than producing an unreliable interval. In this setting, c,dc,d offer a clear interpretation as stage-wise sensitivity to shifts (Section˜7.1). Appendix C.2 also shows how the IID assumption on ScalS^{\text{cal}} can be relaxed to allow stationary ϕ\phi-mixing sequences.

6 Adaptive risk control with residual decomposition

We consider an adaptive variant of our method for nonstationary data by updating the prediction intervals over time. The sets SconfS^{\text{conf}} and ScalS^{\text{cal}} are defined using a sliding window over the most recent observations, with a user-specified window-length kk. We construct intervals combining Definition˜4 and Definition˜3: C^α,a,b,c,d(w)=μ^2(μ^1(w))±(aQ1c({ΔR1})+bQ1d({R2})),\hat{C}_{\alpha,a,b,c,d}(w)=\hat{\mu}_{2}(\hat{\mu}_{1}(w))\pm(a\,Q_{1-c}(\{\Delta R_{1}\})+b\,Q_{1-d}(\{R_{2}\})), parameterized by scaling coefficients a,ba,b, quantile levels c,dc,d, and target coverage level α\alpha, where quantiles are taken over SconfS^{\text{conf}}, and Λvalt\Lambda_{\text{val}}^{t} is recalculated at each time step. With each new point, we dynamically update the (at,bt,ct,dt,αt)(a_{t},b_{t},c_{t},d_{t},\alpha_{t}) using adaptive rules based on recent performance: when coverage drops below target, we decrease αt\alpha_{t} and vice versa; when specific components show persistent errors we adjust their corresponding scaling parameters; and when Λvalt\Lambda_{\text{val}}^{t} is too restrictive we adjust the quantile levels to expand future options. The algorithmic details are provided in Appendix Section˜B.5.

Crucially, the SelectLambda algorithm implements the adaptive adjustments to (at,bt,ct,dt)(a_{t},b_{t},c_{t},d_{t}): it first identifies which component constitutes more of the error by comparing recent averages ΔR1¯=1ki=tkt1ΔR1(i)\bar{\Delta R_{1}}=\frac{1}{k}\sum_{i=t-k}^{t-1}\Delta R_{1}^{(i)} and R2¯=1ki=tkt1R2(i)\bar{R_{2}}=\frac{1}{k}\sum_{i=t-k}^{t-1}R_{2}^{(i)} over the sliding window. If coverage fails and upstream errors dominate (ΔR¯1>R¯2\Delta\bar{R}_{1}>\bar{R}_{2}), it seeks to increase ata_{t} within the validated set Λvalt\Lambda_{\text{val}}^{t}. Conversely, when downstream errors dominate (R¯2>ΔR¯1\bar{R}_{2}>\Delta\bar{R}_{1}), it prioritizes btb_{t}. If the desired scaling adjustment is unavailable in Λvalt\Lambda_{\text{val}}^{t}, the algorithm returns signals Δct,Δdt{1,0,1}\Delta c_{t},\Delta d_{t}\in\{-1,0,1\} to adjust ct+1,dt+1c_{t+1},d_{t+1}, effectively “adjusting" the constraints to allow more suitable scaling options.

Even with additional parameters, this approach preserves the long-run coverage guarantee from prior work: limT1Tt=1Tcovta.s.1α\lim_{T\to\infty}\frac{1}{T}\sum_{t=1}^{T}\text{cov}_{t}\xrightarrow{\text{a.s.}}1-\alpha, where covt\text{cov}_{t} denotes coverage at time tt. This is because the core αt\alpha_{t} update rule remains identical to the adaptive conformal method of Gibbs and Candes [14], resulting in the convergence guarantee formally stated in Appendix A.1.1.

7 Experiments

We evaluate our method on synthetic and real-world forecasting tasks to assess its ability to (i) maintain coverage, (ii) attribute predictive error to specific model stages, and (iii) adapt to distribution shifts in modular pipelines. Our experiments are designed to isolate upstream, downstream, and full-pipeline shifts, highlighting how stage-aware intervals improve transparency and robustness over existing conformal baselines. Additional experiments and details are in Appendix D.

7.1 Non-adaptive methods

We evaluate our non-adaptive method (Section˜4), using the FWER-based procedure for our unified interval which we denote as SRa,b\text{SR}_{a,b} to generate a validated set Λval\Lambda_{\text{val}} of (a,b)(a,b), using fixed quantile levels c,dc,d. We then select from this set the pair (a,b)(a,b) that yields coverage closest to the nominal level α=0.1\alpha=0.1 on SconfS^{\text{conf}}. We compare against two nonadaptive baselines: standard split conformal prediction (SC, Vovk et al. [28]) and weighted split conformal prediction (WSC, Barber et al. [7]), which accounts for non-exchangeability.

We begin with a synthetic dataset generated by the structural equations x=3w+ν1x=3w+\nu_{1}, y=4x+ν2y=4x+\nu_{2}, where w,ν1,ν2𝒩(0,0.1)w,\nu_{1},\nu_{2}\sim\mathcal{N}(0,0.1). After an initial stationary period, the data undergoes a concept shift in either the upstream or downstream with increasing noise. We vary the rate of these shifts to include both gradual and rapid changes. Results for both upstream and downstream shifts are shown in Table˜1.

Table 1: Average coverage and interval width for upstream and downstream shifts under gradual and rapid shifts at α=0.1\alpha=0.1, over 50 samples, with std for SR, SC, WSC methods. NA indicates that our method abstains from providing an interval
(a) Upstream shift coverage and width
Shift Metric SRa,b (c=d=0.01) SRa,b (c=0.05,d=0.01) SC WSC
Gradual Coverage 0.8419±0.010.8419\pm 0.01 0.8233±0.01\pm 0.01 0.7280±0.02\pm 0.02 0.7488±0.02\pm 0.02
Width 2.9884±0.10\pm 0.10 2.4272±0.10\pm 0.10 2.3726±0.10\pm 0.10 2.4849±0.09\pm 0.09
Rapid Coverage 0.8041±0.01\pm 0.01 NA 0.6910±0.01\pm 0.01 0.7218±0.01\pm 0.01
Width 3.9703±0.17\pm 0.17 NA 2.9544±0.11\pm 0.11 3.1969±0.15\pm 0.15
(b) Downstream shift coverage and width
Shift Metric SRa,b (c=d=0.01) SRa,b (c=0.01,d=0.05) SC WSC
Gradual Coverage 0.9309±0.02\pm 0.02 0.9133±0.02\pm 0.02 0.8695±0.03\pm 0.03 0.8748±0.02\pm 0.02
Width 2.4274±0.12\pm 0.12 2.3738±0.12\pm 0.12 2.0143±0.06\pm 0.06 2.0472±0.06\pm 0.06
Rapid Coverage 0.8981±0.01\pm 0.01 NA 0.8419±0.02\pm 0.02 0.8506±0.02\pm 0.02
Width 2.4878±0.11\pm 0.11 NA 2.1035±0.08\pm 0.08 2.1498±0.08\pm 0.08

Our method achieves higher coverage than the other methods, which suffer degradation under distributional shifts, particularly under upstream. For example, our method drops from 0.84 to 0.80 coverage vs. SC dropping from 0.73 to 0.69 for gradual vs. rapid upstream shifts. These coverage improvements are at the cost of wider intervals and occasional abstention; this trade-off is worthwhile when attribution and robustness are prioritized. We note that selecting larger cc (resp. dd) increases the sensitivity of the method to upstream (resp. downstream shift), resulting in abstention under rapid shifts. This is intentional, with abstention signaling that retraining may be necessary for the given stage; black-box methods lack this ability, as they cannot isolate which part of the model fails, resulting in unnecessary retraining.

Thus, cc and dd control sensitivity: by selecting larger or smaller cc and dd, practitioners can balance robustness and abstention at each stage, offering flexibility beyond traditional conformal approaches. Next, we evaluate the adaptive version of our method under more extreme distributional shifts.

7.2 Adaptive methods under structured distribution shifts

We evaluate our adaptive method (Section˜6) in settings with distributional shocks, comparing against recent adaptive baselines: Adaptive Conformal Inference (ACI, Gibbs and Candes [14]), Proportional-Integral-Derivative control (PID, Angelopoulos et al. [1]), and Online Conformal Inference with Decaying step sizes (OCID, Angelopoulos et al. [4]). We also consider DtACI (Gibbs and Candès [15]), an extension of ACI, and report its results in Appendix D.3.2 due to weaker performance. While DtACI is designed for online adaptation, it may be less effective in fixed-horizon forecasting tasks where the lag between prediction and feedback complicates adaptation. We also consider covariate shifts, hyperparameter sweeps, and a real-world stock price dataset, with full experimental details in Appendix D.

Controlled two-stage regression with targeted shifts.

We simulate a two-stage regression pipeline, with controlled stochastic relationships between upstream input ww, intermediate representation xx, and downstream output yy. The initial data follows an i.i.d. structure with x=μ1(w)=3w+ν1,y=μ2(x)=4x+ν2x=\mu_{1}(w)=3w+\nu_{1},\quad y=\mu_{2}(x)=4x+\nu_{2} where ν1,ν2𝒩(0,1)\nu_{1},\nu_{2}\sim\mathcal{N}(0,1). Both stages are modeled using ordinary least squares. At test time, we simulate a temporal sequence of three phases: (i) Upstream concept shift: μ1\mu_{1} becomes μ1shift(w)=8w+1+ν1\mu_{1}^{\text{shift}}(w)=8w+1+\nu_{1}; (ii) Reversion: The upstream returns to μ1(w)\mu_{1}(w); (iii) Downstream concept shift: μ2\mu_{2} becomes μ2shift(x)=7x+5+ν2\mu_{2}^{\text{shift}}(x)=7x+5+\nu_{2}. We report coverage and widths in Figure˜2, capturing the robustness of our methods to stage-specific shocks without overcompensating on width, particularly for upstream shocks.

Refer to caption
(a) Widths of the four intervals for α=0.1\alpha=0.1, δ=0.1\delta=0.1, γ=0.01\gamma=0.01, η=0.01\eta=0.01, k=100k=100
Refer to caption
(b) Coverage averaged over sliding window of 200 points
Figure 2: Comparison of interval width and coverage under synthetic distribution shifts
Automobile indicators dataset.

We evaluate our method on a real-world dataset consisting of monthly economic indicators for the U.S. automobile supply chain and used vehicle valuations. Upstream features (w)(w) are sourced from FRED [13], including metrics such as price indices, inventory levels, and manufacturing prices, which are then forecasted (x)(x) at a 6-month horizon. Downstream outputs (y)(y), aggregated used car prices, are obtained from Manheim [17]. This setting fits the two-stage modeling framework: upstream economic conditions influence downstream pricing, and training a separate upstream predictor leverages more widely available historical supply chain data.

A notable feature of this dataset is the sharp shock in 2020 due to the COVID-19 pandemic, which disrupted global supply chains, leading to prolonged production halts, resulting in a steep and persistent rise in used vehicle prices [19]. The nature of the shock may also affect the model asymmetrically, leading to exacerbated upstream inaccuracy, which we observe via the scaling parameter changes in Appendix D.3.5. For the upstream forecasting task, we use an ensemble of VAR and ARIMA models to predict supply-chain indicators. The model is trained on a rolling-window basis, with residual scores computed at each time step. A 40-month sliding window forms the held-out set from which we derive SconfS^{\text{conf}} and ScalS^{\text{cal}}.

We present the resulting prediction intervals for our method, ACI, PID, and OCID in Figure˜3. Our method effectively responds to the shock in 2020, maintaining coverage relative to other methods. OCID is also able to capture the shock, but shows inefficiency during 2014-2019. Notably, our approach maintains a comparable average width (21.085 vs 15.696 (ACI), 18.469 (PID), and 31.058 (OCID)) and can quickly adjust to shocks by consistently considering conservative values of λ\lambda in Λval\Lambda_{\text{val}}. This suggests that our method is well-suited for real-world forecasting and shines under asymmetric structural shocks. In contrast, for covariate shifts (Appendix D) that affect the entire pipeline, the benefits are less significant.

Refer to caption
(a) SR prediction intervals from 2013-2023
Refer to caption
(b) ACI prediction intervals from 2013-2023
Refer to caption
(c) PID prediction intervals from 2013-2023
Refer to caption
(d) OCID prediction intervals from 2013-2023
Figure 3: Coverage of prediction intervals for α=0.2\alpha=0.2, δ=0.3\delta=0.3, γ=0.01\gamma=0.01, η=0.01\eta=0.01, k=40k=40

8 Conclusion

We proposed a conformal prediction framework for sequential models that decomposes prediction residuals into interpretable stage-specific components. Our method combines this decomposition with risk-controlled parameter selection to construct prediction intervals that provide both coverage guarantees and uncertainty attribution. The approach allows practitioners to identify which pipeline stage contributes to prediction uncertainty and provides diagnostic tools for targeted model improvement. Empirical evaluation on synthetic distribution shifts and real-world economic forecasting shows that our method maintains coverage under shifts that degrade standard conformal approaches. While this robustness comes at the cost of wider intervals and occasional abstention when shifts are severe, these trade-offs reflect the method’s conservative design that prioritizes reliable coverage over optimistic predictions. The stage-wise decomposition proves particularly valuable for asymmetric shifts affecting different pipeline components, enabling targeted adaptive responses that standard methods cannot provide.

References

  • Angelopoulos et al. [2023] Anastasios Angelopoulos, Emmanuel Candes, and Ryan J Tibshirani. Conformal pid control for time series prediction. Advances in neural information processing systems, 36:23047–23074, 2023.
  • Angelopoulos et al. [2021] Anastasios N Angelopoulos, Stephen Bates, Emmanuel J Candès, Michael I Jordan, and Lihua Lei. Learn then test: Calibrating predictive algorithms to achieve risk control. arXiv preprint arXiv:2110.01052, 2021.
  • Angelopoulos et al. [2022] Anastasios N Angelopoulos, Stephen Bates, Adam Fisch, Lihua Lei, and Tal Schuster. Conformal risk control. arXiv preprint arXiv:2208.02814, 2022.
  • Angelopoulos et al. [2024] Anastasios N Angelopoulos, Rina Foygel Barber, and Stephen Bates. Online conformal prediction with decaying step sizes. arXiv preprint arXiv:2402.01139, 2024.
  • Athreya and Pantula [1986] Krishna B Athreya and Sastry G Pantula. A note on strong mixing of arma processes. Statistics & probability letters, 4(4):187–190, 1986.
  • Barber [2024] Rina Foygel Barber. Hoeffding and bernstein inequalities for weighted sums of exchangeable random variables. arXiv preprint arXiv:2404.06457, 2024.
  • Barber et al. [2023] Rina Foygel Barber, Emmanuel J Candes, Aaditya Ramdas, and Ryan J Tibshirani. Conformal prediction beyond exchangeability. The Annals of Statistics, 51(2):816–845, 2023.
  • Bates et al. [2021] Stephen Bates, Anastasios Angelopoulos, Lihua Lei, Jitendra Malik, and Michael Jordan. Distribution-free, risk-controlling prediction sets. Journal of the ACM (JACM), 68(6):1–34, 2021.
  • Bauer [1991] Peter Bauer. Multiple testing in clinical trials. Statistics in medicine, 10(6):871–890, 1991.
  • Bonferroni [1937] Carlo Emilio Bonferroni. On the multiple hypotheses test, 1937. origin of the Bonferroni correction.
  • Crane and Crotty [1967] Dwight B Crane and James R Crotty. A two-stage forecasting model: Exponential smoothing and multiple regression. Management Science, 13(8):B–501, 1967.
  • Farinhas et al. [2023] António Farinhas, Chrysoula Zerva, Dennis Ulmer, and André FT Martins. Non-exchangeable conformal risk control. arXiv preprint arXiv:2310.01262, 2023.
  • FRED [2025] FRED. automobile indicators, 2025. data retrieved from FRED, https://fred.stlouisfed.org/series/CUSR0000SETA02.
  • Gibbs and Candes [2021] Isaac Gibbs and Emmanuel Candes. Adaptive conformal inference under distribution shift. Advances in Neural Information Processing Systems, 34:1660–1672, 2021.
  • Gibbs and Candès [2024] Isaac Gibbs and Emmanuel J Candès. Conformal inference for online prediction with arbitrary distribution shifts. Journal of Machine Learning Research, 25(162):1–36, 2024.
  • Kontorovich and Ramanan [2008] Leonid Kontorovich and Kavita Ramanan. Concentration inequalities for dependent random variables via the martingale method. Annals of Probability, 2008.
  • Manheim [2025] Manheim. Used vehicle value index, 2025. data retrieved from Cox Automative, https://site.manheim.com/en/services/consulting/used-vehicle-value-index.html.
  • Mohri and Rostamizadeh [2010] Mehryar Mohri and Afshin Rostamizadeh. Stability bounds for stationary φ\varphi-mixing and β\beta-mixing processes. Journal of Machine Learning Research, 11(2), 2010.
  • Mullin [2022] John Mullin. Supply chain disruptions, inflation, and the fed. Econ Focus, 22(3Q):14–17, 2022.
  • Oliveira et al. [2024] Roberto I Oliveira, Paulo Orenstein, Thiago Ramos, and João Vitor Romano. Split conformal prediction and non-exchangeable data. Journal of Machine Learning Research, 25(225):1–38, 2024.
  • Quach et al. [2023] Victor Quach, Adam Fisch, Tal Schuster, Adam Yala, Jae Ho Sohn, Tommi S Jaakkola, and Regina Barzilay. Conformal language modeling. arXiv preprint arXiv:2306.10193, 2023.
  • Romano et al. [2019] Yaniv Romano, Evan Patterson, and Emmanuel Candes. Conformalized quantile regression. Advances in neural information processing systems, 32, 2019.
  • Serfling [1974] Robert J Serfling. Probability inequalities for the sum in sampling without replacement. The Annals of Statistics, pages 39–48, 1974.
  • Shafer and Vovk [2008] Glenn Shafer and Vladimir Vovk. A tutorial on conformal prediction. Journal of Machine Learning Research, 9(3), 2008.
  • Soybilgen and Yazgan [2021] Barış Soybilgen and Ege Yazgan. Nowcasting us gdp using tree-based ensemble models and dynamic factors. Computational Economics, 57(1):387–417, 2021.
  • Tibshirani [2023] Ryan J. Tibshirani. Conformal prediction, 2023. URL https://www.stat.berkeley.edu/˜ryantibs/statlearn-s23/lectures/conformal.pdf. Lecture notes for Advanced Topics in Statistical Learning, Spring 2023, University of California, Berkeley.
  • Vovk [2012] Vladimir Vovk. Conditional validity of inductive conformal predictors. In Asian conference on machine learning, pages 475–490. PMLR, 2012.
  • Vovk et al. [2005] Vladimir Vovk, Alexander Gammerman, and Glenn Shafer. Algorithmic learning in a random world, volume 29. Springer, 2005.
  • Wu et al. [2025] Junxi Wu, Dongjian Hu, Yajie Bao, Shu-Tao Xia, and Changliang Zou. Error-quantified conformal inference for time series. arXiv preprint arXiv:2502.00818, 2025.
  • Zaffran et al. [2022] Margaux Zaffran, Olivier Féron, Yannig Goude, Julie Josse, and Aymeric Dieuleveut. Adaptive conformal predictions for time series. In International Conference on Machine Learning, pages 25834–25866. PMLR, 2022.
  • Zargarbashi et al. [2023] Soroush H Zargarbashi, Simone Antonelli, and Aleksandar Bojchevski. Conformal prediction sets for graph neural networks. In International Conference on Machine Learning, pages 12292–12318. PMLR, 2023.
  • Zhang et al. [2025] Botong Zhang, Shuo Li, and Osbert Bastani. Conformal structured prediction. In ICLR 2025 Workshop on Building Trust in Language Models and Applications, 2025. URL https://openreview.net/forum?id=mKfmLQXP6J.

Appendix A Appendix

A.1 Main Results/Proofs

Proof of Theorem 1
Proof.

We consider the event

ytestC^c,d(wtest)\displaystyle y_{test}\notin\hat{C}_{c,d}(w_{test})

where C^c,d\hat{C}_{c,d} is the prediction interval defined by μ^2(μ^1(wtest))±Q1c({ΔR1})+Q1d({R2})\hat{\mu}_{2}(\hat{\mu}_{1}(w_{test}))\pm Q_{1-c}(\{\Delta R_{1}\})+Q_{1-d}(\{R_{2}\}) with the quantiles being taken over the residual components calculated on held-out conformal set SconfS^{\text{conf}}. This event is when the coverage of the interval fails. See that

ytestC^c,d|ytestμ^2(μ^1(wtest))|Q1c({ΔR1})+Q1d({R2})\displaystyle y_{test}\notin\hat{C}_{c,d}\implies|y_{test}-\hat{\mu}_{2}(\hat{\mu}_{1}(w_{test}))|\geq Q_{1-c}(\{\Delta R_{1}\})+Q_{1-d}(\{R_{2}\})

But then by definition

|ytestμ^2(μ^1(wtest))|Q1c({ΔR1})+Q1d({R2})\displaystyle|y_{test}-\hat{\mu}_{2}(\hat{\mu}_{1}(w_{test}))|\geq Q_{1-c}(\{\Delta R_{1}\})+Q_{1-d}(\{R_{2}\})\implies
||ytestμ^2(xtest)||ytestμ^2(μ^1(wtest))||+|ytestμ^2(xtest)|Q1c({ΔR1})+Q1d({R2})\displaystyle||y_{test}-\hat{\mu}_{2}(x_{test})|-|y_{test}-\hat{\mu}_{2}(\hat{\mu}_{1}(w_{test}))||+|y_{test}-\hat{\mu}_{2}(x_{test})|\geq Q_{1-c}(\{\Delta R_{1}\})+Q_{1-d}(\{R_{2}\})

Thus, the split residual on the testtest point also falls outside the region. However, this implies that at least one of the following events occur:

A1:|ytestμ^2(xtest)|Q1d({R2})A_{1}:|y_{test}-\hat{\mu}_{2}(x_{test})|\geq Q_{1-d}(\{R_{2}\})

or

A2:||ytestμ^2(xtest)||ytestμ^2(μ^1(wtest))||Q1c({ΔR1})A_{2}:||y_{test}-\hat{\mu}_{2}(x_{test})|-|y_{test}-\hat{\mu}_{2}(\hat{\mu}_{1}(w_{test}))||\geq Q_{1-c}(\{\Delta R_{1}\})

Then each of these pieces hold due to exchangeable properties as the probability of the new residual piece being higher in rank than the 1c1-c quantile (resp. 1d1-d) is cc (resp. dd). For R2R_{2}, the exchangeability is clear as it is directly the prediction interval for μ2\mu_{2} on data (x,y)(x,y). For the latter, it is still clearly exchangeable as well: the functions μ^1,μ^2\hat{\mu}_{1},\hat{\mu}_{2} are symmetric and deterministic and are pretrained, thus being fixed on the held-out set. Thus ΔR1\Delta R_{1} can be interpreted as the result of some deterministic function ψ\psi on each point (wi,xi,yi)(w_{i},x_{i},y_{i}), which implies that

P(ΔR11\displaystyle P(\Delta R_{1}^{1} =r1,,ΔR1m=rm)\displaystyle=r_{1},\dots,\Delta R_{1}^{m}=r_{m})
=P((w1,x1,y1)ψ1(r1),,(wm,xm,ym)ψ1(rm))\displaystyle=P((w_{1},x_{1},y_{1})\in\psi^{-1}(r_{1}),\dots,(w_{m},x_{m},y_{m})\in\psi^{-1}(r_{m}))
=P((wπ(1),xπ(1),yπ(1))ψ1(rπ(1)),,(wπ(m),xπ(m),yπ(m))ψ1(rπ(m)))\displaystyle=P((w_{\pi(1)},x_{\pi(1)},y_{\pi(1)})\in\psi^{-1}(r_{\pi(1)}),\dots,(w_{\pi(m)},x_{\pi(m)},y_{\pi(m)})\in\psi^{-1}(r_{\pi(m)}))
=P(ΔR1π(1)=rπ(1),,ΔR1π(m)=rπ(m))\displaystyle=P(\Delta R_{1}^{\pi(1)}=r_{\pi(1)},\dots,\Delta R_{1}^{\pi(m)}=r_{\pi(m)})

where ΔR1i\Delta R_{1}^{i} denotes the first residual component of the ii-th point of SconfS^{\text{conf}}. Then we have exchangeability of ΔR1\Delta R_{1}, thus the probability of each of the noncontainment events for each residual component are bounded by cc and dd respectively.

Therefore

P(ytestC^c,d)c+dP(A1A2)c+d.P(y_{test}\notin\hat{C}_{c,d})\leq c+d-P(A_{1}\cap A_{2})\leq c+d.

In fact, as can be seen in the above equation, the stated result is looser than necessary. ∎

Proof of Proposition 1
Proof.

Since the traditional residual RR is assumed to be continuous, the function that maps the scaling parameters (a,b)(a,b) to the marginal coverage probability of the interval C^α,a,b\hat{C}_{\alpha,a,b} is itself continuous. When a=b=0a=b=0, the interval degenerates to a single point and thus has zero coverage (assuming the problem is nontrivial). When a=b=1a=b=1, the interval is wide enough to ensure coverage of at least 12α1-2\alpha by Corollary 1. By the intermediate value theorem, there must exist a pair (a,b)[0,1]2(a^{*},b^{*})\in[0,1]^{2} such that the interval achieves exactly 1α1-\alpha coverage. As an aside, one can further view the problem as an optimization problem

mina,baQ1α({ΔR1})+bQ1α({R2})\displaystyle\min_{a,b}aQ_{1-\alpha}(\{\Delta R_{1}\})+bQ_{1-\alpha}(\{R_{2}\})
s.t.P(RaQ1α({ΔR1})+bQ1α({R2}))α\displaystyle s.t.\quad P(R\geq aQ_{1-\alpha}(\{\Delta R_{1}\})+bQ_{1-\alpha}(\{R_{2}\}))\leq\alpha

which by KKT conditions implies that any a,ba^{*},b^{*} must satisfy

1Q1α({ΔR1})P(RaQ1α({ΔR1})+bQ1α({R2})))a\displaystyle\frac{1}{Q_{1-\alpha}(\{\Delta R_{1}\})}\frac{\partial P(R\geq aQ_{1-\alpha}(\{\Delta R_{1}\})+bQ_{1-\alpha}(\{R_{2}\})))}{\partial a}
=1Q1α({R2})P(RaQ1α({ΔR1})+bQ1α({R2})))b\displaystyle=\frac{1}{Q_{1-\alpha}(\{R_{2}\})}\frac{\partial P(R\geq aQ_{1-\alpha}(\{\Delta R_{1}\})+bQ_{1-\alpha}(\{R_{2}\})))}{\partial b}

which provides that at any optimal solution a,ba^{*},b^{*} there must be the same balance between the magnitude of components and the impact (derivative) of the scaling parameter on true coverage probability for each parameter. Thus large magnitude implies large impact, confirming intuitive understanding. ∎

Proof of Theorem 2
Proof.

By Lemma˜A.2, the pp-values are super-uniform under the null hypothesis. Then we may apply Theorem 1 of [2] (Theorem˜3), given a FWER algorithm 𝒜\mathcal{A} with parameter δ>0\delta>0 to obtain the result. We note that the size of ScalS^{\text{cal}} does not explicitly appear in the bound, however, it directly affects the pp-values which can be seen via Hoeffding’s inequality:

P(αBin(l,α)αl^(λ))exp(2(αl^(λ))2l)P(\alpha-\text{Bin}(l,\alpha)\geq\alpha-l\hat{\mathcal{R}}(\lambda))\leq\exp\left(\frac{-2(\alpha-l\hat{\mathcal{R}}(\lambda))^{2}}{l}\right)

where ^\hat{\mathcal{R}} is the empirical error rate on ScalS^{\text{cal}}

A.1.1 Long-Run Coverage

Proposition 3 (Long-Run Coverage of Adaptive Method).

Let γ>0\gamma>0 be the step size, α(0,1)\alpha\in(0,1) the nominal coverage level, and let covt{0,1}\text{cov}_{t}\in\{0,1\} denote the coverage indicator at time tt. Then the adaptive algorithm satisfies:

limT1Tt=1Tcovta.s.1α.\lim_{T\to\infty}\frac{1}{T}\sum_{t=1}^{T}\text{cov}_{t}\xrightarrow{\text{a.s.}}1-\alpha.

This result ensures that in the long run, the algorithm achieves the desired coverage level 1α1-\alpha, without requiring distributional assumptions on the data.

The additional parameters (at,bt,ct,dt)(a_{t},b_{t},c_{t},d_{t}) are updated as deterministic functions of the algorithm’s history, ensuring they do not interfere with the update steps that underlie the coverage guarantee.

A finite-sample convergence bound (Appendix Proposition˜6) extends existing Markov chain analysis to include our additional parameters, with the key insight that the monotonicity properties required still hold under our parameter update rules.

Proof of Proposition 3

Please refer to Algorithm˜1 for notation. We first show a preliminary lemma.

Lemma A.1 (Boundedness of ctc_{t} and dtd_{t}).

Let η>0\eta>0 be a constant. Then, for all tt\in\mathbb{N}, the sequences ctc_{t} and dtd_{t} are bounded within the interval [η,1+η][-\eta,1+\eta]. That is,

t,ηct1+ηandηdt1+η.\forall t\in\mathbb{N},\quad-\eta\leq c_{t}\leq 1+\eta\quad\text{and}\quad-\eta\leq d_{t}\leq 1+\eta.
Proof.

To establish the lower bound for ctc_{t}, observe that the argument for dtd_{t} follows similarly, and the upper-bound can be treated identically.

Suppose that ct<0c_{t}<0 for some time step tt. In this case, by the definition of the residual component ΔR1\Delta R_{1}, we have:

ΔR1=covt=1,\Delta R_{1}=\infty\implies\text{cov}_{t}=1,

where covt\text{cov}_{t} denotes the coverage of the model at time tt. Additionally, since covtΔR1=0\text{cov}_{t}^{\Delta R_{1}}=0 by the construction of the residual, this implies that in the output of the algorithm SelectLambda, the update step Δct0\Delta c_{t}\geq 0, meaning that:

ct+1ct.c_{t+1}\geq c_{t}.

Thus, if ctc_{t} were negative at some time step, the next value is non-decreasing. Consequently, the sequence ctc_{t} cannot decrease below by η-\eta. Therefore, the minimum value attainable by ctc_{t} is η-\eta.

A similar argument holds for the upper-bound.

Hence, we conclude that:

ηct1+ηandηdt1+ηt.-\eta\leq c_{t}\leq 1+\eta\quad\text{and}\quad-\eta\leq d_{t}\leq 1+\eta\quad\forall t\in\mathbb{N}.

This holds similarly for dtd_{t}. ∎

Next, we show the asymptotic result, Section˜A.1.1 via the same arguments as Proposition 4.1 of [14].

Proof.

Although our algorithm includes additional parameters (at,bt,ct,dt)(a_{t},b_{t},c_{t},d_{t}), the update rule for αt\alpha_{t} remains unchanged from [14]:

αt=αt1γ((1covt1)α),\alpha_{t}=\alpha_{t-1}-\gamma((1-\mathrm{cov}_{t-1})-\alpha),

where covt\mathrm{cov}_{t} is the coverage indicator at time tt and α\alpha is the target coverage level.

As in [14], we show that αt[γ,1+γ]\alpha_{t}\in[-\gamma,1+\gamma] with probability 1. For the lower bound (the upper-bound follows symmetrically), suppose αt<0\alpha_{t}<0 at some time tt. Then the candidate set Λval=\Lambda_{\mathrm{val}}=\emptyset, so covt=1\mathrm{cov}_{t}=1, and the update becomes:

αt+1=αt+γα>αt.\alpha_{t+1}=\alpha_{t}+\gamma\alpha>\alpha_{t}.

Hence, αt\alpha_{t} is pushed upward, preventing it from decreasing without bound. Therefore, αtγ\alpha_{t}\geq-\gamma for all tt.

Next, unfolding the recursion gives:

αT=α1t=1T1γ((1covt)α).\alpha_{T}=\alpha_{1}-\sum_{t=1}^{T-1}\gamma((1-\mathrm{cov}_{t})-\alpha).

Rearranging terms and using the bound, we obtain:

|1Tt=1T(1covt)α|max(α1,1α1)+γγT.\left|\frac{1}{T}\sum_{t=1}^{T}(1-\mathrm{cov}_{t})-\alpha\right|\leq\frac{\max(\alpha_{1},1-\alpha_{1})+\gamma}{\gamma T}.

Equivalently:

|1Tt=1Tcovt(1α)|max(α1,1α1)+γγT,\left|\frac{1}{T}\sum_{t=1}^{T}\mathrm{cov}_{t}-(1-\alpha)\right|\leq\frac{\max(\alpha_{1},1-\alpha_{1})+\gamma}{\gamma T},

which converges to zero as TT\to\infty, yielding the desired result. ∎

A.2 Auxiliary Results

Proposition 4.

For a given point (w,x,y)(w,x,y), if R2ϵ<RR_{2}\leq\epsilon<R for some small ϵ>0\epsilon>0 i.e. μ2(x)\mu_{2}(x) is close to the true value yy, then |ΔR1R|ϵ|\Delta R_{1}-R|\leq\epsilon

Proof.

This follows immediately by definition of ΔR1\Delta R_{1}.

|ΔR1R|\displaystyle|\Delta R_{1}-R| =||R2R|R|\displaystyle=||R_{2}-R|-R|
ϵ\displaystyle\leq\epsilon

Next, we show the equivalent for R2R_{2} under additional assumptions

Proposition 5.

For a given point (w,x,y)(w,x,y), |R2R|<ϵ|R_{2}-R|<\epsilon for small ϵ>0\epsilon>0 under certain assumptions, such as if μ^2\hat{\mu}_{2} is a Lipschitz continuous function for some parameter LL, and the first stage has small prediction error, |xμ^1(w)|<δ|x-\hat{\mu}_{1}(w)|<\delta. Then |R2R|Lδ|R_{2}-R|\leq L\delta

Proof.
|R2R|\displaystyle|R_{2}-R| |μ^2(x)μ^2(μ^1(w))|\displaystyle\leq|\hat{\mu}_{2}(x)-\hat{\mu}_{2}(\hat{\mu}_{1}(w))|
L|xμ^1(w)|\displaystyle\leq L|x-\hat{\mu}_{1}(w)|
Lδ\displaystyle\leq L\delta

The main takeaway of these two results is that they confirm the intuition for when ΔR1\Delta R_{1} and R2R_{2} are roughly similar to RR, particularly when one model is accurate (and in the case of R2R_{2} when the learned function is smooth). This also implies that both ΔR1,R2\Delta R_{1},R_{2} being small is impossible, thus one must represent a majority of the error. Note that ΔR1\Delta R_{1} and R2R_{2} can also serve as rough lower bounds on RR, though each can individually exceed RR under adversarial noise or model miscalibration.

Lemma A.2.

The pλp_{\lambda} defined as pλ=(Bin(n,α)n^(λ))p_{\lambda}=\mathbb{P}\left(\text{Bin}(n,\alpha)\leq n\hat{\mathcal{R}}(\lambda)\right) are super-uniform.

Proof.

See Appendix D.1 of [21]. ∎

Theorem 3 (Theorem 1 of [2]).

Suppose pλp_{\lambda} has a distribution stochastically dominating the uniform distribution. Let 𝒜\mathcal{A} be a FWER controlling algorithm at level δ>0\delta>0. Let (λ)\mathcal{R}(\lambda) be the expected risk for the choice of λ\lambda and α>0\alpha>0 be the maximal error rate. Then Λval=𝒜({pλ}λΛ)\Lambda_{\text{val}}=\mathcal{A}(\{p_{\lambda}\}_{\lambda\in\Lambda}) satisfies the following

(supλΛval{(λ)}α)δ\displaystyle\mathbbm{P}\left(\sup_{\lambda\in\Lambda_{\text{val}}}\{\mathcal{R}(\lambda)\}\geq\alpha\right)\leq\delta

A.3 Extensions of Theorem 1

We can also state Theorem˜1 for the setting with auxiliary data (Section˜C.3).

Corollary 2.

For c,d(0,1)c,d\in(0,1) and new point (wtest,xtest,xtest,ytest)(w_{\text{test}},x_{\text{test}},x_{\text{test}}^{\prime},y_{\text{test}}), and C^c,d\hat{C}_{c,d} interval defined as in Section˜C.3 to include auxiliary data, we have

P(ytestC^c,d(wtest))1cd.\displaystyle P(y_{\text{test}}\in\hat{C}_{c,d}(w_{\text{test}}))\geq 1-c-d.

which follows by the same proof as above; auxiliary data does not change the proof (assuming that with auxiliary features, the data is still exchangeable). Lastly, we note that the prior results can easily be adapted to nonexchangable data settings using weighted quantiles and a weighting scheme as mentioned in [7]. Thus we have

Corollary 3.

Let |Sconf|=m|S^{\text{conf}}|=m and p1,pmp_{1}\dots,p_{m} be weights for each point. For c,d(0,1)c,d\in(0,1) and new point (wtest,xtest,ytest)(w_{\text{test}},x_{\text{test}},y_{\text{test}}), and C^c,d\hat{C}_{c,d} interval defined as in Definition˜3, we have

P(ytestC^c,d)1cdψ1ψ2.\displaystyle P(y_{\text{test}}\in\hat{C}_{c,d})\geq 1-c-d-\psi_{1}-\psi_{2}.

where ψ1\psi_{1} and ψ2\psi_{2} are coverage penalties from nonexchangability.

Proof.

The result follows analogously to the proof of Section˜A.1, with modifications to account for the use of weighted quantiles and associated penalty terms. Specifically, we let ψ1\psi_{1} and ψ2\psi_{2} denote penalties applied to the empirical weighted quantiles, as introduced in Theorem 2 of [7]. These penalties are of the form i=1mp~idTV(R,Ri))\sum_{i=1}^{m}\tilde{p}_{i}d_{TV}(\vec{R},\vec{R}^{i})) where p~i=pi1+i=1mpi\tilde{p}_{i}=\frac{p_{i}}{1+\sum_{i=1}^{m}p_{i}}, R\vec{R} represents the vector of full residuals and Ri\vec{R}^{i} is the same vector with the ii-th point swapped with (wtest,ytest,ztest)(w_{\text{test}},y_{\text{test}},z_{\text{test}}).

When applying weighted conformal prediction, the standard quantile threshold used in the unweighted case is replaced with a weighted quantile. By Theorem 2 of [7], these ψ1,ψ2\psi_{1},\psi_{2} ensure that the resulting prediction set retains valid marginal coverage.

Therefore, by adapting the same steps as in Section˜A.1—but replacing the empirical quantiles with weighted quantiles and applying Theorem 2 of [7], we obtain the stated result. ∎

We also introduce a result for coverage for fixed scaling values of a,ba,b under some assumptions.

Corollary 4.

Assume that the CDF of the residual for an out-of-sample point is Lipschitz continuous with constant LL. Furthermore, if the maximal value of {ΔR1}\{\Delta R_{1}\} (resp. {R2}\{R_{2}\}) is bounded by ϵ>0\epsilon>0, then if we set a=0,b=1a=0,b=1 (resp. a=1,b=0a=1,b=0) with c=d=αc=d=\alpha, we have coverage with

P(ytestC^0,1,α(wtest))2α+Lϵ.P(y_{\text{test}}\notin\hat{C}_{0,1,\alpha}(w_{\text{test}}))\leq 2\alpha+L\epsilon.
Proof.

Consider C^1,1,α(wtest)=Q1α({ΔR1})+Q1α({R2})\hat{C}_{1,1,\alpha}(w_{\text{test}})=Q_{1-\alpha}(\{\Delta R_{1}\})+Q_{1-\alpha}(\{R_{2}\}), which has coverage guarantee

P(ytestC^1,1,α(wtest))2α.P(y_{\text{test}}\notin\hat{C}_{1,1,\alpha}(w_{\text{test}}))\leq 2\alpha.

By assumption, Q1α({ΔR1})<ϵQ_{1-\alpha}(\{\Delta R_{1}\})<\epsilon, and thus instead using C^0,1,α(wtest)\hat{C}_{0,1,\alpha}(w_{\text{test}}) results in a change of at most ϵ\epsilon, which by Lipschitz continuity of the CDF produces a resultant change of LϵL\epsilon in the probability guarantee. ∎

The existence of the Lipschitz condition is quite reasonable, for example, we consider μ2\mu_{2} is a linear function and the noise term is normally distributed, the residuals themselves are half-normally distributed and thus satisfy the above requirement. We pay a small theoretical guarantee when shrinking the interval, under the assumption LϵL\epsilon is small, which is not necessarily true (and LL is generally unknown). However, this result gives the intuition that if one of the residual pieces contributes little to the overall error, we can remove it for little coverage cost and serves as part of the motivation for the adaptive algorithm. We find that for certain settings with unbalanced residual components, even simple heuristics such as max(Q1c({ΔR1}),Q1d({R2}))\max(Q_{1-c}(\{\Delta R_{1}\}),Q_{1-d}(\{R_{2}\})) where c=d=αc=d=\alpha achieves similar coverage to the full residuals.

We can obtain non-asymptotic bounds for the specific setting of a hidden Markov Model, as in Theorem 4.1 of [14] once again, with little changes due to the shared update step.

Proposition 6.

Suppose that the data is coming from a hidden Markov model, with underlying states At𝒜A_{t}\in\mathcal{A} such that the data (wt,xt,yt)PAt(w_{t},x_{t},y_{t})\sim P_{A_{t}}. Observe that at,bt,ct,dt,αt,Ata_{t},b_{t},c_{t},d_{t},\alpha_{t},A_{t} form a Markov chain on [0,1]2×[η,1+η]2×[γ,1+γ]×𝒜[0,1]^{2}\times[-\eta,1+\eta]^{2}\times[-\gamma,1+\gamma]\times\mathcal{A} with stationary distribution π\pi. Assume that the algorithm has passed its burn-in period and is now in a stationary setting. Suppose that the spectral gap of {At}=1η>0\{A_{t}\}=1-\eta>0. Let errt=1covt\mathrm{err}_{t}=1-\text{cov}_{t}. Then let B=supa𝒜|E[errt|At=a]α|B=\sup_{a\in\mathcal{A}}|E[\mathrm{err}_{t}|A_{t}=a]-\alpha| and σB2=E[(E[errt|At]α)2]\sigma_{B}^{2}=E[(E[\mathrm{err}_{t}|A_{t}]-\alpha)^{2}]. Then

P(|1Tt=1Terrtα|ϵ)2exp(Tϵ2/8)+2exp(T(1η)ϵ28(1+η)σB2+20Bϵ)P(|\frac{1}{T}\sum_{t=1}^{T}\mathrm{err}_{t}-\alpha|\geq\epsilon)\leq 2\exp(-T\epsilon^{2}/8)+2\exp\left(\frac{-T(1-\eta)\epsilon^{2}}{8(1+\eta)\sigma_{B}^{2}+20B\epsilon}\right)

showing a finite-sample bound deviation bound for the empirical coverage.

Proof.

The structure of this proof follows that of Theorem 4.1 in [14], which depends on their Lemma A.2 and Theorem A.1. We outline the necessary modifications to their lemma (written below) to hold in our setting.

Lemma A.3.

For any τ\tau\in\mathbb{R} and t𝒩t\in\mathcal{N}

𝔼[s=1texp(τ(errs𝔼[errs|As]))]exp(τ2/2)𝔼[s=1t1exp(τ(errs𝔼[errs|As]))]\displaystyle\mathbbm{E}\left[\prod_{s=1}^{t}\exp(\tau(\mathrm{err}_{s}-\mathbbm{E}[\mathrm{err}_{s}|A_{s}]))\right]\leq\exp(\tau^{2}/2)\mathbbm{E}\left[\prod_{s=1}^{t-1}\exp(\tau(\mathrm{err}_{s}-\mathbbm{E}[\mathrm{err}_{s}|A_{s}]))\right]

To show this, we condition on (a1,b1,c1,d1)(a_{1},b_{1},c_{1},d_{1}) alongside (α1,A1,,At)(\alpha_{1},A_{1},\dots,A_{t}). As in [14], we require that

f(s=1t1errs)=s=1t1exp(τ(errs𝔼[errs|As]))f(\sum_{s=1}^{t-1}\mathrm{err}_{s})=\prod_{s=1}^{t-1}\exp(\tau(\mathrm{err}_{s}-\mathbbm{E}[\mathrm{err}_{s}|A_{s}]))

is non-decreasing and

g(s=1t1errs)=𝔼[exp(τ(errt𝔼[errt|At]))|err1,,errt1]=At(ytC^at,bt,ct,dt,αt(wt))exp(τ𝔼[errt|At])\displaystyle g(\sum_{s=1}^{t-1}\mathrm{err}_{s})=\mathbbm{E}\left[\exp(\tau(\mathrm{err}_{t}-\mathbbm{E}[\mathrm{err}_{t}|A_{t}]))|\mathrm{err}_{1},\dots,\mathrm{err}_{t-1}\right]=\mathbb{P}_{A_{t}}(y_{t}\in\hat{C}_{a_{t},b_{t},c_{t},d_{t},\alpha_{t}}(w_{t}))\exp(-\tau\mathbbm{E}[\mathrm{err}_{t}|A_{t}])
+(1At(ytC^at,bt,ct,dt,αt(wt)))exp(τ(1𝔼[errt|At]))\displaystyle\quad+(1-\mathbb{P}_{A_{t}}(y_{t}\in\hat{C}_{a_{t},b_{t},c_{t},d_{t},\alpha_{t}}(w_{t})))\exp(\tau(1-\mathbbm{E}[\mathrm{err}_{t}|A_{t}]))

is non-increasing in s=1t1errs\sum_{s=1}^{t-1}\mathrm{err}_{s} for τ0\tau\geq 0. This still holds because:

  • αt\alpha_{t} remains a non-increasing function of s=1t1errs\sum_{s=1}^{t-1}\mathrm{err}_{s} under the same update rule.

  • a,bt,ct,dta_{,}b_{t},c_{t},d_{t} are deterministic functions of the conditioning variables (including err1,,errt1\mathrm{err}_{1},\dots,\mathrm{err}_{t-1}), as the candidate set selection and parameter updates are fixed given the history.

  • As αt\alpha_{t} decreases, the validation set Λval\Lambda_{\mathrm{val}} shrinks (selecting only higher-width pairs), so the interval width is non-decreasing, thus At(ytC^at,bt,ct,dt,αt(wt))\mathbb{P}_{A_{t}}(y_{t}\in\hat{C}_{a_{t},b_{t},c_{t},d_{t},\alpha_{t}}(w_{t})) is non-decreasing, meaning that gg is overall non-increasing.

These observations ensure that the monotonicity conditions in Lemma A.2 still hold under our conditioning. The remainder of the proof, including the application of the Markov chain concentration inequality (Theorem A.1 of [14]), proceeds unchanged. ∎

Appendix B Conceptual information

B.1 Discussion on quantiles of residual components.

We observe that when quantiles of the residual components are taken then summed, they do not necessarily serve as upper-bounds on the corresponding quantiles of the full residuals, which we use in our implementation. Specifically, the inequality

Q1α({ΔR1})+Q1α({R2})<Q1α({R})Q_{1-\alpha}(\{\Delta R_{1}\})+Q_{1-\alpha}(\{R_{2}\})<Q_{1-\alpha}(\{R\})

may not hold. This observation is significant because while the triangle inequality holds for each individual point, i.e., RΔR1+R2R\leq\Delta R_{1}+R_{2}, such a relationship does not extend to the quantiles. This discrepancy is intentional, as it allows for a tighter combination of residuals compared to the traditional residual. Moreover, in practice, this situation occurs rarely.

B.2 Notation Table

We provide a table of notation in Table˜2.

Table 2: Summary of Notation
Symbol Meaning
μ^1\hat{\mu}_{1} First stage hypothesis
μ^2\hat{\mu}_{2} Second stage hypothesis
RR Residual |yμ^2(μ^1(w))||y-\hat{\mu}_{2}(\hat{\mu}_{1}(w))| for two-stage model
ΔR1\Delta R_{1} First stage residual component ||yμ^2(x)||yμ^2(x^)||\left|\,|y-\hat{\mu}_{2}(x)|-|y-\hat{\mu}_{2}(\hat{x})|\,\right|
R2R_{2} Second stage residual component |yμ^2(x)||y-\hat{\mu}_{2}(x)|
α\alpha Desired Miscoverage Rate
aa Scaled weight for ΔR1\Delta R_{1}
bb Scaled weight for R2R_{2}
cc Quantile level to take of ΔR1\Delta R_{1}
dd Quantile level to take of R2R_{2}
δ\delta Rejection stringency to control FWER
γ\gamma Step-size for updating αt\alpha_{t} with γ=0.01\gamma=0.01 empirically performing well
η\eta Step-size for updating ct,dtc_{t},d_{t} with η=0.01\eta=0.01 performing well
τ\tau Tolerance for calculating pp-values for hypothesis testing
kk Window length for adaptive method, with stabilized performance above 100 points

B.3 Choice of Components

Certainly, one could instead take quantiles or rescalings of the total sum ΔR1+R2\Delta R_{1}+R_{2} instead of separating it into components, however we choose to separate it multiple reasons. First, we note that by splitting the quantiles/weights, it becomes clear in the decomposition which stage contributes more to the error. Second under co-monotonicity of ΔR1,R2\Delta R_{1},R_{2}, the quantiles of ΔR1,R2\Delta R_{1},R_{2} match and thus both approaches result in the same interval Q1α({ΔR1+R2})Q_{1-\alpha}(\{\Delta R_{1}+R_{2}\}) which satisfies coverage guarantees as it is an upper-bound on the black box width. Furthermore, for smaller miscoverage levels α1,α2\alpha_{1},\alpha_{2} then Q1α1({ΔR1})+Q1α2({R2})Q_{1-\alpha_{1}}(\{\Delta R_{1}\})+Q_{1-\alpha_{2}}(\{R_{2}\}) is an upperbound on Q1α(R)Q_{1-\alpha}(R), so our use of Q1α({ΔR1})+Q1α({R2})Q_{1-\alpha}(\{\Delta R_{1}\})+Q_{1-\alpha}(\{R_{2}\}) can be an upper-bound. Lastly, we also desire the ability to be tighter than Q1α(R)Q_{1-\alpha}(R) at times, unlike Q1α({ΔR1+R2})Q_{1-\alpha}(\{\Delta R_{1}+R_{2}\}), to be able to quickly adapt to easy settings.

B.4 FWER-controlling algorithms for Scaling Selection

To construct the set of valid scaling parameter pairs Λval\Lambda_{\text{val}}, we apply multiple hypothesis testing procedures that control the family-wise error rate (FWER). This ensures that, with high probability, none of the accepted scaling pairs exhibit a true miscoverage rate above the nominal level α\alpha.

Definition 5 (FWER-Controlling Algorithm).

Let Λ={λ1,,λu}\Lambda=\{\lambda_{1},\dots,\lambda_{u}\} be a candidate set of scaling pairs, and let p1,,pu[0,1]p_{1},\dots,p_{u}\in[0,1] denote the associated pp-values testing the null hypothesis that the miscoverage of each λi\lambda_{i} exceeds α\alpha. A selection procedure 𝒜:[0,1]u2{1,,u}\mathcal{A}:[0,1]^{u}\to 2^{\{1,\dots,u\}} is said to control the family-wise error rate (FWER) at level δ\delta if:

(𝒜(p1,,pu)J)1δ,\mathbb{P}\left(\mathcal{A}(p_{1},\dots,p_{u})\subseteq J\right)\geq 1-\delta,

where J{1,,u}J\subseteq\{1,\dots,u\} is the set of true null hypotheses.

Bonferroni Correction.

A classical method for controlling FWER is the Bonferroni correction [10]. It accepts any λi\lambda_{i} for which piδ/up_{i}\leq\delta/u, where u=|Λ|u=|\Lambda|. By the union bound, this ensures that the probability of falsely including any λ\lambda with miscoverage above α\alpha is at most δ\delta. However, Bonferroni can be overly conservative when the number of tests uu is large.

Fixed Sequence Testing.

To improve power in large-scale settings, we adopt fixed sequence testing [9]. This procedure assumes a pre-specified ordering of hypotheses p(1),p(2),,p(u)p_{(1)},p_{(2)},\dots,p_{(u)}, based on prior knowledge or heuristics (e.g., larger values of a+ba+b are expected to have lower risk). Hypotheses are tested sequentially: each p(i)p_{(i)} is compared to δ\delta, and the first failure (i.e., p(i)>δp_{(i)}>\delta) terminates the process. All previous hypotheses are accepted, and the rest are rejected. This method is more powerful than Bonferroni in the presence of a meaningful ordering.

Practical Note.

In our algorithm, we use fixed sequence testing to select valid scaling parameters. If no candidate passes the test, i.e., Λval=\Lambda_{\text{val}}=\emptyset, we abstain from prediction. This outcome indicates a lack of statistical evidence and may signal that additional calibration data or model refinement is necessary.

B.5 Adaptive method pseudocode

Algorithm 1 Adaptive Risk Controlling with ΔR1\Delta R_{1}, R2R_{2}
1:Window size kk, Step sizes γ\gamma, η\eta; target coverage level α\alpha
2:for time step tt do
3:  Obtain window of kk recent observations
4:  Obtain Parameters αt\alpha_{t}, ctc_{t}, dtd_{t}, previous coverage covt1,covt1ΔR1,covt1R2,at1,bt1\text{cov}_{t-1},\text{cov}^{\Delta R_{1}}_{t-1},\text{cov}^{R_{2}}_{t-1},a_{t-1},b_{t-1}
5:  Split window into SconfS^{\text{conf}} and ScalS^{\text{cal}} and compute {ΔR1},{R2}\{\Delta R_{1}\},\{R_{2}\}
6:  Compute quantiles:
7:    ΔR1Qct({ΔR1})\Delta R_{1}\leftarrow Q_{c_{t}}(\{\Delta R_{1}\}), R2Qdt({R2})R_{2}\leftarrow Q_{d_{t}}(\{R_{2}\})
8:  Determine valid λ\lambda set:
9:    ΛvaltFixedSequenceTesting(Scal,αt,ΔR1,R2)\Lambda^{t}_{\text{val}}\leftarrow\textsc{FixedSequenceTesting}(S^{\text{cal}},\alpha_{t},\Delta R_{1},R_{2})
10:  Obtain at,bt,Δct,ΔdtSelectLambda(Λvalt,covt1,covt1ΔR1,covt1R2,at1,bt1)a_{t},b_{t},\Delta c_{t},\Delta d_{t}\leftarrow\textsc{SelectLambda}(\Lambda^{t}_{\text{val}},\text{cov}_{t-1},\text{cov}^{\Delta R_{1}}_{t-1},\text{cov}^{R_{2}}_{t-1},a_{t-1},b_{t-1})
11:  Calculate interval C^α,at,bt,ct,dt(wt)\hat{C}_{\alpha,a_{t},b_{t},c_{t},d_{t}}(w_{t})
12:  Check coverage at tt: covt,covtΔR1,covtR2\text{cov}_{t},\text{cov}^{\Delta R_{1}}_{t},\text{cov}^{R_{2}}_{t}
13:  Update: αt+1αt+γ(αcovt)\alpha_{t+1}\leftarrow\alpha_{t}+\gamma(\alpha-\text{cov}_{t})
14:     ct+1ct+ηΔctc_{t+1}\leftarrow c_{t}+\eta\Delta c_{t},  dt+1dt+ηΔdtd_{t+1}\leftarrow d_{t}+\eta\Delta d_{t}
15:end for

We define the coverage indicator at time tt as covt=𝟙ytC^α,a,b,c,d(wt).\text{cov}_{t}=\mathbbm{1}_{y_{t}\in\hat{C}_{\alpha,a,b,c,d}(w_{t})}. To capture deviations in the score components, we define covtΔR1\text{cov}^{\Delta R_{1}}_{t} to be 0 if the true value (ΔR1)t(\Delta R_{1})_{t} lies within the estimated quantile interval for ΔR1\Delta R_{1}, and equal to the excess (ΔR1)tΔR1(\Delta R_{1})_{t}-\Delta R_{1} otherwise. This can be succinctly written as covtΔR1=ReLU((ΔR1)tΔR1),\text{cov}^{\Delta R_{1}}_{t}=\text{ReLU}((\Delta R_{1})_{t}-\Delta R_{1}), and similarly for covtR2\text{cov}^{R_{2}}_{t}.

Appendix C Extensions

C.1 Additional Heuristic

One could alternatively define the residual components and combination as signed residual components to obtain tighter, asymmetric intervals, at the cost of coverage. We define new residual components

R~2\displaystyle\tilde{R}_{2} =yμ^2(x)\displaystyle=y-\hat{\mu}_{2}(x)
ΔR~1\displaystyle\Delta\tilde{R}_{1} =μ^2(μ^1(w))μ^2(x)\displaystyle=\hat{\mu}_{2}(\hat{\mu}_{1}(w))-\hat{\mu}_{2}(x)

such that the sum of the components is exactly the full error rather than an absolute value upper-bound. Because the components are now signed, one should construct intervals differently, which also results in asymmetric intervals. Similarly to before, we can scale each component to produce intervals of the form

[μ^2(μ^1(wtest)+aQc/2({ΔR~1})+bQd/2({R~2}),μ^2(μ^1(wtest)+aQ1c/2({ΔR~1})+bQ1d/2({R~2})]\displaystyle[\hat{\mu}_{2}(\hat{\mu}_{1}(w_{\text{test}})+aQ_{c/2}(\{\Delta\tilde{R}_{1}\})+bQ_{d/2}(\{\tilde{R}_{2}\}),\hat{\mu}_{2}(\hat{\mu}_{1}(w_{\text{test}})+aQ_{1-c/2}(\{\Delta\tilde{R}_{1}\})+bQ_{1-d/2}(\{\tilde{R}_{2}\})]

Experimentally in IID settings, it produces slightly tighter intervals without additional tolerance, however this heuristic is more prone to producing empty Λval\Lambda_{\text{val}} sets. For example, we report the performance of the method under IID settings with varying levels of nominal error rate α\alpha in Table˜3 exemplifying this behaviour.

Table 3: Signed residual heuristic performance over varying α\alpha
Metric α=0.13\alpha=0.13 α=0.01\alpha=0.01 α=0.07\alpha=0.07 α=0.05\alpha=0.05
Coverage 0.9013±\pm0.02 0.9196±\pm0.01 NA NA
Width 2.0265±\pm0.13 2.1523±\pm0.12 NA NA

Thus we have chosen to use the heuristic given in the main body of the text, but we provide this as an option for many possible ways of constructing intervals using the residual decomposition we introduce. Furthermore, because the decomposition itself is flexible, one could incorporate any recent conformal prediction advancements into the ways quantiles are taken, such as weighted data, or in the adaptive case how the step-size and αt\alpha_{t} etc. change.

C.2 Stationary Φ\Phi-mixing

We describe how the FWER control method outlined in Section˜5 can be extended beyond the IID setting to accommodate stationary ϕ\phi-mixing processes. Note that for dependent data, the concentration properties of empirical quantities degrade, resulting in looser tail bounds. This, in turn, inflates the values of the empirical error estimates pλp_{\lambda}, thereby requiring larger thresholds δ\delta to ensure that the selected set Λval\Lambda_{\text{val}} remains nonempty.

Our goal is to identify values of λ\lambda for which the expected coverage error is below a target level α\alpha. For each λ\lambda, we consider the null hypothesis

H0λ:(ytestC^λ(wtest))>α,H_{0}^{\lambda}:\quad\mathbb{P}\left(y_{\text{test}}\notin\hat{C}_{\lambda}(w_{\text{test}})\right)>\alpha,

and compute an empirical estimate of the coverage error using a held-out calibration dataset. To test this hypothesis, we apply a concentration inequality that bounds the deviation of the empirical error from its expectation under dependence. Specifically, we employ a result for bounded functions of ϕ\phi-mixing sequences from [16], with refinements from [18], which provides suitable control over the error rates in the non-IID case.

Theorem 4 (Mixing Concentration Inequality).

Let Z1,,ZnZ_{1},...,Z_{n} be random variables distributed according to a ϕ\phi-mixing distribution. Let f:ZnRf:Z^{n}\to R be a measurable function that is cc-Lipschitz with respect to the Hamming metric for some c>0c>0. Then, for any ϵ>0\epsilon>0, the following inequality holds:

P(|f(Z1,,Zn)𝔼[f(Z1,,Zn)]|ϵ)2exp(2ϵ2nc2Δn2)\displaystyle P(|f(Z_{1},\dots,Z_{n})-\mathbbm{E}[f(Z_{1},\dots,Z_{n})]|\geq\epsilon)\leq 2\exp{\left(\frac{-2\epsilon^{2}}{nc^{2}\left\lVert\Delta_{n}\right\rVert_{\infty}^{2}}\right)}

where Δn=1+i=1nϕ(i)\left\lVert\Delta_{n}\right\rVert_{\infty}=1+\sum_{i=1}^{n}\phi(i)

We apply this theorem with ^(λ)=f((w1,x1,y1),,(wl,xl,yl))=1li=1l1{yiC^λ(wi)}\hat{\mathcal{R}}(\lambda)=f((w_{1},x_{1},y_{1}),\dots,(w_{l},x_{l},y_{l}))=\frac{1}{l}\sum_{i=1}^{l}1_{\{y_{i}\notin\hat{C}_{\lambda}(w_{i})\}}. Clearly this is 1l\frac{1}{l}-Lipschitz with respect to Hamming metric. Furthermore,

𝔼[f(Scal)]=1l𝔼[i=1l1{yiC^λ(wi)}]=1li=1l𝔼[1{yiC^λ(wi)}]=α\mathbbm{E}[f(S^{\text{cal}})]=\frac{1}{l}\mathbbm{E}[\sum_{i=1}^{l}1_{\{y_{i}\notin\hat{C}_{\lambda}(w_{i})\}}]=\frac{1}{l}\sum_{i=1}^{l}\mathbbm{E}[1_{\{y_{i}\notin\hat{C}_{\lambda}(w_{i})\}}]=\alpha

due to stationarity and the null hypothesis. Then

P(𝔼[^(λ)]^(λ)ϵ)\displaystyle P(\mathbbm{E}[\hat{\mathcal{R}}(\lambda)]-\hat{\mathcal{R}}(\lambda)\geq\epsilon) 2exp(2ϵ2lΔl2)\displaystyle\leq 2\exp{\left(\frac{-2\epsilon^{2}}{l\left\lVert\Delta_{l}\right\rVert_{\infty}^{2}}\right)}

where we substitute our realized α^(λ)\alpha-\hat{\mathcal{R}}(\lambda) for ϵ\epsilon and plug it into the RHS to obtain a pp-value for λ\lambda.

Then we have the following:

Corollary 5.

Let ΛvalΛ\Lambda_{\text{val}}\subseteq\Lambda be the set selected by a FWER-controlling algorithm at level δ\delta, based on pp-values computed over ScalS^{\text{cal}} using the above method, which is assumed to be stationary ϕ\phi-mixing. Then, for any λ^Λval\hat{\lambda}\in\Lambda_{\text{val}}, we have

((ytestC^λ^(wtest)|Scal)α)1δ,\mathbb{P}\left(\mathbb{P}\left(y_{\text{test}}\notin\hat{C}_{\hat{\lambda}}(w_{\text{test}})\,\big|\,S^{\text{cal}}\right)\leq\alpha\right)\geq 1-\delta,

where the outer probability is over the randomness of ScalS^{\text{cal}}, and the inner probability is over the test point (wtest,xtest,ytest)(w_{\text{test}},x_{\text{test}},y_{\text{test}}) which is assumed to be drawn from the same stationary distribution.

While the statement of the result is similar to before, note that pλp_{\lambda} is an conservative upper-bound on the true pp-value, thus the threshold to accept a λ\lambda is more stringent. We show that the pλp_{\lambda} are still super-uniform. Let PP be the random variable representing empirical risk when the error rate is α\alpha and QQ by the empirical risk when the error rate is α\alpha^{\prime}. Under the null hypothesis, α>α\alpha^{\prime}>\alpha, which implies that QQ stochastically dominates PP (P(Qt)P(Pt)P(Q\leq t)\leq P(P\leq t)). Let f(q)=exp(2(αq)2lΔn2)f(q)=\exp{\left(\frac{-2(\alpha-q)^{2}}{l\left\lVert\Delta_{n}\right\rVert_{\infty}^{2}}\right)} and FP(t)F_{P}(t) and FQ(t)F_{Q}(t) represent CDFs of PP and QQ respectively. Then

pλ\displaystyle p_{\lambda} =f(Q)\displaystyle=f(Q)
P(αPαQ)\displaystyle\geq P(\alpha-P\geq\alpha-Q)
=FP(Q)\displaystyle=F_{P}(Q)

Furthermore, letting O=pλO=p_{\lambda}, by the stochastic dominance assumption,

P(Oo)\displaystyle P(O\leq o) =P(pλo)\displaystyle=P(p_{\lambda}\leq o)
P(FP(Q)o)\displaystyle\leq P(F_{P}(Q)\leq o)
P(FQ(Q)o)\displaystyle\leq P(F_{Q}(Q)\leq o)
=P(QFQ1(o))\displaystyle=P(Q\leq F_{Q}^{-1}(o))
=FQ(FQ1(o))\displaystyle=F_{Q}(F_{Q}^{-1}(o))
=o\displaystyle=o

Then we apply a FWER-controlling algorithm such as one from Section˜B.4 to set the acceptance threshold to reject the null hypothesis such that the total type-1 error rate is less than δ\delta, then use Theorem 1 of [2] (Theorem˜3).

We note that it is possible to adapt this result to exchangeable sequences using specific concentration inequalities such as ones found in [23] or more recently [6].

C.3 Decomposition for auxiliary data

We can also define ΔR1\Delta R_{1} and R2R_{2} when considering auxiliary data.

Definition 6.

If we have sample SS where data is of the form (w,x,x,y)(w,x,x^{\prime},y) where xx^{\prime} represents auxiliary data for the second stage such that we have learned upstream and downstream models μ^1:𝒲𝒳\hat{\mu}_{1}:\mathcal{W}\to\mathcal{X}, μ^2:𝒳×𝒳𝒴\hat{\mu}_{2}:\mathcal{X}\times\mathcal{X}^{\prime}\to\mathcal{Y}. Then we can decompose the residual into components as the following for some new point (w,x,x,y)(w,x,x^{\prime},y)

ΔR1(w,x,x,y)\displaystyle\Delta R_{1}(w,x,x^{\prime},y) =||yμ^2(x,x)||yμ^2(μ^1(w),x)||\displaystyle=||y-\hat{\mu}_{2}(x,x^{\prime})|-|y-\hat{\mu}_{2}(\hat{\mu}_{1}(w),x^{\prime})||
R2(w,x,x,y)\displaystyle R_{2}(w,x,x^{\prime},y) =|yμ^2(x,x)|.\displaystyle=|y-\hat{\mu}_{2}(x,x^{\prime})|.

We can use these residual components to create prediction intervals analogously to the original setting without auxiliary data.

C.4 Decomposition for multi-stage model

Now we consider extending the definition of ΔR1,R2\Delta R_{1},R_{2} to more than two stages. We provide a straightforward extension, but a cleverer one may exist.

Definition 7.

Suppose we have sample SS where data is of the form w1,w2,,wN+1w_{1},w_{2},\dots,w_{N+1} with NN stages and corresponding learned hypotheses μ^i:𝒲i𝒲i+1\hat{\mu}_{i}:\mathcal{W}_{i}\to\mathcal{W}_{i+1}. Then for a point (w1,,wN+1)(w_{1},\dots,w_{N+1}), define

ΔR1(w1,,wN+1)\displaystyle\Delta R_{1}(w_{1},\dots,w_{N+1}) =|||wN+1μ^N(μ^N1(μ^2(w2)))||wN+1μ^N(μ^N1(μ^1(w1)))||\displaystyle=|||w_{N+1}-\hat{\mu}_{N}(\hat{\mu}_{N-1}(\dots\hat{\mu}_{2}(w_{2})))|-|w_{N+1}-\hat{\mu}_{N}(\hat{\mu}_{N-1}(\dots\hat{\mu}_{1}(w_{1})))||
\displaystyle\dots
ΔRi(w1,,wN+1)\displaystyle\Delta R_{i}(w_{1},\dots,w_{N+1}) =|||wN+1μ^N(μ^N1(μ^i+1(wi+1)))||wN+1μ^N(μ^N1(μ^i(wi)))||\displaystyle=|||w_{N+1}-\hat{\mu}_{N}(\hat{\mu}_{N-1}(\dots\hat{\mu}_{i+1}(w_{i+1})))|-|w_{N+1}-\hat{\mu}_{N}(\hat{\mu}_{N-1}(\dots\hat{\mu}_{i}(w_{i})))||
\displaystyle\dots
Rn(w1,,wN+1)\displaystyle R_{n}(w_{1},\dots,w_{N+1}) =|wN+1μ^n(wN)|\displaystyle=|w_{N+1}-\hat{\mu}_{n}(w_{N})|

Observe that the sum of these terms is still an upperbound on the "full" residual R=|wN+1μ^N(μ^N1(μ^1(w1)))|R=|w_{N+1}-\hat{\mu}_{N}(\hat{\mu}_{N-1}(\dots\hat{\mu}_{1}(w_{1})))| by the triangle inequality and thus we have a residual decomposition.

An issue with this definition is that if we try a hypothesis testing method and FWER algorithm to weight each of these components, the set Λ\Lambda becomes prohibitively large as it grows exponentially with each stage. A FWER algorithm that is able to intelligently search the possibilities will become necessary, whereas for two-stages, it was feasible to search across the entire grid. One way to mitigate this is to restrict the proposed set Λ\Lambda at each step by focusing only on the most important residual components, fixing the weight of the remaining ones to remain; i.e. select the top kk residual components with the highest average width in the window of recent observations and define Λ\Lambda only for those components, freezing the other component weights. In fact, one can sort the residual components by magnitude then sequentially search through Λ\Lambda, freezing the others outside of the top kk, resulting in only kdkd tests, if dd is the number of sub-divisions for each component.

C.5 Decomposition for multiple upstream models

We consider an augmented two-stage model, with multiple upstream models each producing a component of xx, the intermediate value, which is now vector valued. Suppose there are NN upstream models, μ^11,,μ^N1\hat{\mu}^{1}_{1},\dots,\hat{\mu}^{1}_{N} that each map upstream features 𝐰=w1,,wM\mathbf{w}=w_{1},\dots,w_{M} to their corresponding intermediate feature in 𝐱=x1,,xN\mathbf{x}=x_{1},\dots,x_{N}. Denote the output of the ii-th upstream model to be x^i\hat{x}_{i}. Then we can define ΔR1,R2\Delta R_{1},R_{2} components for each downstream feature and take a weighted sum. Thus, for the ii-th downstream feature, we have

(R2)i\displaystyle(R_{2})_{i} =|yμ^2(x^1,,xi,x^N)|\displaystyle=|y-\hat{\mu}_{2}(\hat{x}_{1},\dots,x_{i}\dots,\hat{x}_{N})|
(ΔR1)i\displaystyle(\Delta R_{1})_{i} =||yμ^2(x^1,,x^N)||yμ^2(x^1,,xi,x^N)||.\displaystyle=||y-\hat{\mu}_{2}(\hat{x}_{1},\dots,\hat{x}_{N})|-|y-\hat{\mu}_{2}(\hat{x}_{1},\dots,x_{i}\dots,\hat{x}_{N})||.

Then see that a convex combination i=1nNθi((R2)i+(ΔR1)i),i=1Nθi=1,θi[0,1]\sum_{i=1}^{nN}\theta_{i}((R_{2})_{i}+(\Delta R_{1})_{i}),\sum_{i=1}^{N}\theta_{i}=1,\theta_{i}\in[0,1] forms an upperbound on R=|yμ^2(x^1,,x^N)|R=|y-\hat{\mu}_{2}(\hat{x}_{1},\dots,\hat{x}_{N})|. This approach has a similar issue as the multi-stage model, as the search-space Λ\Lambda scales with the number of upstream models.

Appendix D Additional experiments

We provide additional experiments on the non-adaptive methods and adaptive methods, providing further insight into the interpretability of the split residual method and the coverage guarantees. We include hyperparameter ablation studies as well as an additional real-world dataset. We include experiments on covariate shift to demonstrate that the split residual method is robust multiple types of shift and is not limited to concept shift, which was extensively explored in the main text. All experiments were conducted on a laptop with an Intel i7 processor and 8GB of RAM. No specialized hardware (e.g., GPUs or TPUs) was used for training or evaluation.

D.1 Experimental hyperparameter details

We describe the hyperparameters we use for the other methods in our experimens: SC, WSC, ACI, DtACI, PID, OCID.

  • SC: There is no other parameter used for split conformal prediction intervals other than the nominal desired miscoverage level α=0.1\alpha=0.1 which we keep the same across all methods

  • WSC uses the same α\alpha and also includes a weighting parameter 0.99 which exponentially weights points as they go further back in time, resulting in recent points having more weight.

  • ACI uses the same initial α=0.1\alpha=0.1 and updates αt\alpha_{t} using the same step size γ=0.01\gamma=0.01 as well as the same window size k=100k=100 as our method.

  • DtACI uses the same α=0.1\alpha=0.1 and also uses default multiple candidate values for step-size γ=[0.001,0.002,0.004,0.008,0.0160,0.032,0.064,0.128]\gamma=[0.001,0.002,0.004,0.008,0.0160,0.032,0.064,0.128].

  • PID uses the same initial α=0.1\alpha=0.1, step size γ=0.01\gamma=0.01, and window size k=100k=100. Two additional parameters KI and CsatC_{\text{sat}} are estimated with appropriate hypothesized bound BB depending on the data, using the default heuristic given in Appendix B of [1].

  • OCID uses the same α=0.1\alpha=0.1 and decay parameter 0.10.1 with the decay weight function given by [4]

Furthermore, note that ScalS^{\text{cal}} for our method is a subset of the held-out data which the other methods have full access to, so we do not use additional data relative to the other methods but simply partition it into smaller sets.

D.2 Non-adaptive experiments

We discuss results for the non-adaptive intervals using our residual decomposition as well as hyperparameter ablation studies.

D.2.1 Coverage under easy settings

We implement both Definition˜4 and Definition˜3 separately and test them on simple IID data with linear relationships between w,x,yw,x,y. We also include a results for an IID synthetic data featuring a nonlinear relationship in which the upstream task is a binary classification task of diagnosing a patient for a disease given their health (using logistic regression), and the downstream model predicts the cost of treatment given the upstream probabilistic prediction (using XGBoost). Additionally, we implement these prediction interval methods in a ϕ\phi-mixing setting which we provide guarantees for in Section˜C.2 where the data is generated by a stationary AR(1) process with bounded uniform noise (known to be a ϕ\phi-mixing process [5]) and linear relationships between w,x,yw,x,y. We report the width and coverage of the methods in Table˜4.

Table 4: Average coverage and width with standard deviation across 50 runs for the scaling components method with fixed c=d=α=0.1c=d=\alpha=0.1 and FWER as in Definition˜4, which we denote as SRa,b\text{SR}_{a,b}, the separate quantiles method with c=d=0.05c=d=0.05 as in Definition˜3, denoted as SRc,d\text{SR}_{c,d}, SC, and WSC. We examine these under linear and non-linear relationships and non-IID data
Shift Metric SRa,b\textbf{SR}_{a,b} SRc,d\textbf{SR}_{c,d} SC WSC
Linear Coverage 0.9439±0.02\pm 0.02 0.9492±0.01\pm 0.01 0.8971±0.01\pm 0.01 0.9062±0.02\pm 0.02
Width 2.3555±0.13\pm 0.13 2.3930±0.07\pm 0.07 1.9924±0.06\pm 0.06 2.0584±0.13\pm 0.13
Non-linear Coverage 0.9416±0.01\pm 0.01 0.9794±0.01\pm 0.01 0.8969±0.01\pm 0.01 0.9300±0.02\pm 0.02
Width 4.6343±0.14\pm 0.14 5.2007±0.06\pm 0.06 4.2179±0.08\pm 0.08 4.5089±0.16\pm 0.16
Non-IID (linear) Coverage 0.9561±0.02\pm 0.02 0.9886±0.01\pm 0.01 0.8997±0.02\pm 0.02 0.9067±0.03\pm 0.03
Width 3.1299±0.20\pm 0.20 3.6106±0.12\pm 0.12 2.6749±0.08\pm 0.08 2.7304±0.15\pm 0.15

We observe that both SRa,b,SRc,d\text{SR}_{a,b},\text{SR}_{c,d} are conservative for these settings, producing wider intervals but attaining the desired α\alpha coverage guarantee. Thus, our methods tend to over-cover in easier settings such as IID data and stationary ϕ\phi-mixing conditions because FWER hypothesis testing requires an error rate below α\alpha whereas for other methods simply being close is sufficient. However, as we see in the ablation studies below, this can be easily remedied if desired by including a tolerance for the hypothesis test itself.

D.2.2 Hyperparameter Study

We vary the tolerance parameter τ\tau which alters the hypothesis test, obtaining Bin(n,α+τ)\text{Bin}(n,\alpha+\tau) to demonstrate that when comparing our method with some tolerance, the width and coverage even under the IID settings are comparable and do not over-cover.

Table 5: Coverage and average width (±\pm std) across 50 runs for varying τ{0.01,0.03,0.05}\tau\in\{0.01,0.03,0.05\}. Only SRa,b\text{SR}_{a,b} depends on τ\tau, thus the other methods do not report values for each τ\tau
Method Metric τ=0.01\tau=0.01 τ=0.03\tau=0.03 τ=0.05\tau=0.05
SRa,b\text{SR}_{a,b} Coverage 0.9240±\pm0.01 0.9150±\pm0.02 0.9019±\pm0.02
Width 2.1750±\pm0.09 2.0922±\pm0.08 2.0192±\pm0.11
SRc,d\text{SR}_{c,d} Coverage 0.9516±\pm0.01
Width 2.4032±\pm0.05
SC Coverage 0.8981±\pm0.01
Width 2.0051±\pm0.04
WSC Coverage 0.9015±\pm0.02
Width 2.0320±\pm0.13

We observe that even low values of the tolerance parameter, such as τ=0.01\tau=0.01, are sufficient to make SRa,b\text{SR}_{a,b} competitive with the black-box methods in the IID setting. This suggests that SRa,b\text{SR}_{a,b} is not simply over-covering, but rather enforcing stricter criteria for valid prediction intervals than the black-box approaches. These stricter requirements likely contribute to its improved performance under distributional shifts. Therefore, if one seeks tighter but potentially less robust intervals, adjusting the tolerance parameter τ\tau provides a principled way to trade-off between coverage robustness and interval sharpness, however experimentally, unless stated otherwise, we use τ=0\tau=0.

We also vary the δ\delta parameter, which controls the FWER rejection threshold, which for smaller values encourages more sensitivity of the method overall to abstain from producing an interval; in a sense it has a similar effect to c,dc,d parameters but collectively for both stages at once rather than a stage-wise sensitivity.

Table 6: Coverage and average width (±\pm std) across 50 runs for varying δ{0.1,0.3,0.5,0.7,0.9}\delta\in\{0.1,0.3,0.5,0.7,0.9\}. Only SRa,b\text{SR}_{a,b} depends on δ\delta
Method Metric δ=0.1\delta=0.1 δ=0.3\delta=0.3 δ=0.5\delta=0.5 δ=0.7\delta=0.7 δ=0.9\delta=0.9
SRa,b\text{SR}_{a,b} Coverage NA NA 0.8969±\pm0.08 0.9294±\pm0.02 0.9268±\pm0.02
Width NA NA 2.2231±\pm0.09 2.2195±\pm0.10 2.1867±\pm0.11
SRc,d\text{SR}_{c,d} Coverage 0.9504±\pm0.01
Width 2.3919±\pm0.05
SC Coverage 0.8992±\pm0.01
Width 1.9987±\pm0.04
WSC Coverage 0.9082±\pm0.02
Width 2.0603±\pm0.14

D.3 Adaptive

We include additional experiments using additional synthetic experiments with covariate shift and real-world stock data to demonstrate the flexibility of our method. Furthermore, we provide visualizations of the scaling parameters to demonstrate the interpretability of the decomposition. In addition, we perform hyperparameter sweeps over the synthetic dataset described in Section˜7. We also include comparisons to DtACI, although DtACI performs poorly due to being designed for one step ahead updates, rather than with a forecasting horizon.

D.3.1 Stocks Dataset

We consider a dataset consisting of daily closing values of three S&\&P 500 stocks, AAPL, AMZN, MSFT with the goal of forecasting the AAPL closing value. We do so via a two-stage approach in which the first-stage is an N-BEATS forecasting model that forecasts the three stocks at a 6-day horizon, then the second-stage is ridge regression that predicts AAPL given the forecasted values to correct for error from the first stage. We rescale the data and compare against ACI, Conformal PID, and OCID methods to showcase the robustness of our method. The coverage of these methods is presented in Figure˜4.

Refer to caption
(a) SR prediction intervals from 2017-2018
Refer to caption
(b) ACI prediction intervals from 2017-2018
Refer to caption
(c) DtACI prediction intervals from 2017-2018
Refer to caption
(d) PID prediction intervals from 2017-2018
Refer to caption
(e) OCID prediction intervals from 2017-2018
Figure 4: Coverage of prediction intervals for α=0.1\alpha=0.1, δ=0.3\delta=0.3, γ=0.01\gamma=0.01, η=0.01\eta=0.01, k=24k=24

We observe that our adaptive method provides better coverage than the other three methods, particularly when the forecasts perform worse in July and October, seeing a much less dramatic drop in performance while not being overly tight. Our method uses almost the same average width (0.1347) as ACI (0.1344), while the other two methods have tighter intervals but significantly worse coverage. We list these results in the table below Table˜7 and also the sliding window coverage plots Figure˜5(b). We see that there are two main shocks, and that in comparison to ACI, our method is able to capture the second shock, while OCID is also able to but has volatile intervals during the stable period prior to the shocks.

Table 7: Average and standard deviation of width over time, and minimum average coverage over a sliding window of the last 100 observations for each method on the stocks dataset
Method Avg. Width Std. Width Min Coverage
SR(adaptive) 0.1348 0.0216 0.89
ACI 0.1345 0.0193 0.86
PID 0.1183 0.0226 0.85
OCID 0.1240 0.0348 0.83
Refer to caption
(a) Width of SR, ACI, PID, OCID
Refer to caption
(b) Sliding window coverage of SR, ACI, PID, OCID
Figure 5: Performance of various intervals on the stocks dataset for α=0.1\alpha=0.1, δ=0.3\delta=0.3, γ=0.01\gamma=0.01, η=0.01\eta=0.01, k=24k=24.

D.3.2 DtACI

We also include results on the DtACI method in comparison to the others for the same synthetic dataset with two shocks used in Section˜7.2, with the plots shown in Figure˜6(a).

Refer to caption
(a) Widths of the intervals for α=0.1\alpha=0.1, δ=0.1\delta=0.1, γ=0.01\gamma=0.01, η=0.01\eta=0.01, k=100k=100
Refer to caption
(b) Coverage over the last 200 points. Our method remains more robust, particularly to the upstream shock
Figure 6: Comparison of interval width and coverage robustness under synthetic distribution shifts with DtACI

We observe that DtACI outputs values of αt>1\alpha_{t}>1 which we clip to 1, resulting in very conservative intervals. This is most likely due to the fact that it is not designed for forecasting with a horizon and suffers in performance because of that. We see similar behavior for the automobile indicators dataset (Figure˜7), in which DtACI achieves good coverage but is wider ( compared to our method, ) because it often outputs αt=1\alpha_{t}=1. Interestingly, this does not hold for the stocks dataset, where DtACI is produces the tightest intervals but has the worst coverage.

Refer to caption
Figure 7: DtACI interval for automobile indicators dataset. It captures the jump in 2020 at the cost of wide intervals throughout

D.3.3 Automobile Visualizations

We also include additional visualizations of width and coverage for the automobile dataset Figure˜8. We note that our average width is strong in efficiency (21.085), only behind PID (18.469) and ACI (15.696) which have poor coverage, whereas OCID and DtACI are all wider (31.058, 31.356), respectively.

Refer to caption
(a) Widths of intervals for α=0.2\alpha=0.2, δ=0.3\delta=0.3, γ=0.01\gamma=0.01, η=0.01\eta=0.01, k=40k=40
Refer to caption
(b) Coverage averaged over sliding window of 40 months
Figure 8: Comparison of interval width and coverage for automobile indicator dataset

D.3.4 Hyperparameter Study

First, we visualize the effect of various γ\gamma values, which affects the rate at which αt\alpha_{t} changes. We find that our method excels at lower values of γ=0.01\gamma=0.01 and performance decays at higher values such as γ=0.1\gamma=0.1, particularly because this causes αt\alpha_{t} to drop to extremely low values, making the algorithm abstain from producing intervals under large shifts. We display the coverage effects for synthetic data in Figure˜9 and for Manheim data in Figure˜10, where we see a similar decay as γ\gamma increases.

Refer to caption
(a) Coverage with γ=0.01\gamma=0.01
Refer to caption
(b) Coverage with γ=0.05\gamma=0.05
Refer to caption
(c) Coverage with γ=0.1\gamma=0.1
Figure 9: We use α=0.1\alpha=0.1, δ=0.1\delta=0.1, η=0.01\eta=0.01, k=100k=100

ACI and PID improve slightly at higher γ\gamma on the synthetic dataset. While one might argue that this implies that one can just use higher values of γ\gamma without the need for the split residuals method, using a lower γ\gamma is “safer" as it does not adjust to noisy signals as significantly. Accordingly, the widths of the other methods become extremely volatile with increased γ\gamma. Thus, our method is able to still catch significant jumps while remaining at a safer, consistent step-size due to always considering conservative scaling parameters in Λval\Lambda_{\text{val}}. Furthermore, even when considering multiple values of γ\gamma, on the Manheim dataset, our method with γ=0.01\gamma=0.01 still outperforms the other methods in coverage at all values of γ\gamma, and the performance of the other methods actually decay in coverage as γ\gamma increases. OCID does not use same step-size parameter so it remains unaffected.

Refer to caption
(a) Coverage with γ=0.01\gamma=0.01
Refer to caption
(b) Coverage with γ=0.05\gamma=0.05
Refer to caption
(c) Coverage with γ=0.1\gamma=0.1
Figure 10: We use α=0.2\alpha=0.2, δ=0.3\delta=0.3, η=0.01\eta=0.01, k=40k=40
Table 8: Coverage, width, and maximum improvements in coverage vs. baselines (ACI, PID, OCID) across 50 runs for varying γ{0.01,0.03,0.05,0.1}\gamma\in\{0.01,0.03,0.05,0.1\}
Metric γ=0.01\gamma=0.01 γ=0.03\gamma=0.03 γ=0.05\gamma=0.05 γ=0.1\gamma=0.1
Width 22.92±\pm0.13 24.75±\pm0.10 21.26±\pm0.12 21.93±\pm0.12
Coverage 0.835±\pm0.02 0.815±\pm0.02 0.595±\pm0.02 0.715±\pm0.01
Max Δ\DeltaCoverage vs ACI 0.075±\pm0.02 0.030±\pm0.01 0.025±\pm0.01 0.020±\pm0.01
Max Δ\DeltaCoverage vs PID 0.075±\pm0.02 0.045±\pm0.01 0.045±\pm0.01 0.025±\pm0.01
Max Δ\DeltaCoverage vs OCID 0.130±\pm0.03 0.090±\pm0.02 0.110±\pm0.02 0.100±\pm0.03
Table 9: Coverage, width, and maximum improvements in coverage vs. baselines (ACI, PID, OCID) across 50 runs for varying η{0.01,0.03,0.05,0.1}\eta\in\{0.01,0.03,0.05,0.1\}.
Metric η=0.01\eta=0.01 η=0.03\eta=0.03 η=0.05\eta=0.05 η=0.1\eta=0.1
Width 23.50±\pm0.12 24.22±\pm0.13 24.16±\pm0.14 23.96±\pm0.10
Coverage 0.835±\pm0.03 0.835±\pm0.02 0.855±\pm0.01 0.860±\pm0.02
Max Δ\DeltaCoverage vs ACI 0.080±\pm0.01 0.065±\pm0.02 0.075±\pm0.01 0.075±\pm0.01
Max Δ\DeltaCoverage vs PID 0.105±\pm0.03 0.090±\pm0.01 0.090±\pm0.02 0.075±\pm0.01
Max Δ\DeltaCoverage vs OCID 0.125±\pm0.02 0.135±\pm0.01 0.115±\pm0.03 0.115±\pm0.02
Table 10: Coverage, width, and maximum improvements in coverage vs. baselines (ACI, PID, OCID) across 50 runs for varying δ{0.05,0.1,0.3,0.5,0.7}\delta\in\{0.05,0.1,0.3,0.5,0.7\}
Metric δ=0.05\delta=0.05 δ=0.1\delta=0.1 δ=0.3\delta=0.3 δ=0.5\delta=0.5 δ=0.7\delta=0.7
Width NA 23.30±\pm0.12 22.70±\pm0.13 24.11±\pm0.09 23.20±\pm0.11
Coverage NA 0.865±\pm0.03 0.840±\pm0.02 0.835±\pm0.01 0.845±\pm0.02
Max Δ\DeltaCoverage vs ACI NA 0.075±\pm0.01 0.065±\pm0.02 0.055±\pm0.01 0.065±\pm0.02
Max Δ\DeltaCoverage vs PID NA 0.080±\pm0.01 0.100±\pm0.02 0.080±\pm0.01 0.070±\pm0.01
Max Δ\DeltaCoverage vs OCID NA 0.115±\pm0.02 0.130±\pm0.02 0.110±\pm0.03 0.105±\pm0.03
Table 11: Coverage, width, and maximum improvements in coverage vs. baselines (ACI, PID, OCID) across 50 runs for varying window length k{10,50,100,150}k\in\{10,50,100,150\}
Metric k=10k=10 k=50k=50 k=100k=100 k=150k=150
Width NA 22.21±\pm0.08 23.24±\pm0.07 24.32±\pm0.13
Coverage NA 0.855±\pm0.02 0.856±\pm0.03 0.835±\pm0.01
Max Δ\DeltaCoverage vs ACI NA 0.050±\pm0.01 0.070±\pm0.02 0.070±\pm0.01
Max Δ\DeltaCoverage vs PID NA 0.080±\pm0.01 0.090±\pm0.02 0.080±\pm0.01
Max Δ\DeltaCoverage vs OCID NA 0.115±\pm0.01 0.100±\pm0.01 0.120±\pm0.01

We also provide tables for other hyperparameters γ,η,δ,k\gamma,\eta,\delta,k in Tables˜10, 9, 10 and 11, fixing them at default values of 0.01,0.01,0.3,1000.01,0.01,0.3,100 respectively when not varying them. We see that increasing γ\gamma has a rough impact of decreasing width, while increasing η\eta has a weaker widening effect, δ\delta has no particular effect on width, but low values encourage abstention, and lastly small window lengths also encourage abstention, and generally longer windows results in slightly larger intervals.

D.3.5 Visualization of residual components

Here, we provide extra visualizations of the ΔR1\Delta R_{1} and R2R_{2} parameters along with how the weights a,ba,b change over time, for the same synthetic setup as Section˜7.2 below in Figures˜11 and 12(a).

Refer to caption
(a) The value of the ΔR1\Delta R_{1} residual component
Refer to caption
(b) The value of the R2R_{2} residual component
Figure 11: A visualization of how ΔR1,R2\Delta R_{1},R_{2} change as shocks occur in the data. Note that they are taken at 1ct,1dt1-c_{t},1-d_{t} quantiles respectively. We use α=0.1\alpha=0.1, δ=0.1\delta=0.1, γ=0.01\gamma=0.01, η=0.01\eta=0.01, k=100k=100.

We observe that when the first shock occurs in the upstream, ΔR1\Delta R_{1} captures this change, then after it returns, ΔR1\Delta R_{1} correspondingly shrinks back to the baseline level. R2R_{2} correspondingly captures the second downstream shock.

Refer to caption
(a) Visualization of scaling parameters a,ba,b over time.
Refer to caption
(b) Ratio of RR to ΔR1\Delta R_{1} and R2R_{2} respectively
Figure 12: Plots of residual component relations over time on the synthetic dataset. We use α=0.1\alpha=0.1, δ=0.1\delta=0.1, γ=0.01\gamma=0.01, η=0.01\eta=0.01, k=100k=100

We also see this relationship reflected in the a,ba,b scaling parameters selected by our method as we see them adapt to the upstream shock, the return, then the downstream shock in a similar fashion (Figure˜12(a)). Lastly, we plot the ratios R/(ΔR1)R/(\Delta R_{1}) and R/R2R/R_{2} over time in Figure˜12(b).

We see that the shifts at t=100,500,900t=100,500,900 are visibly identifiable: initially R/R2R/R_{2} is high during the upstream shift as the majority of error comes from R1R_{1}. When the shock returns, we observe that the upstream remains the main source of error, and finally in the downstream shift, R2R_{2} becomes the dominant source of error. Concretely, ΔR1\Delta R_{1} dominates between t=100t=100 to t=500t=500 when the upstream undergoes a shock and R2R_{2} dominates when the downstream undergoes a shock from t=900t=900 onward. Thus we can observe which stage is responsible for the error with extreme clarity allowing diagnostic retraining. In particular, if we retrain at t=200t=200 after the first upstream shift, we can choose to either retrain just μ^1\hat{\mu}_{1} or both μ^1\hat{\mu}_{1} and μ^2\hat{\mu}_{2}. We observe that just retraining μ^1\hat{\mu}_{1} results in a width decrease of 24.643924.6439 for the SR method, while retraining both models, results in a decrease of 25.3451, demonstrating that it is sufficient for our framework to only retrain the affected upstream.

We similarly plot the used scaling parameters (a,b)(a,b) over time for the Manheim used vehicles dataset, and observe that there is in fact an asymmetry in the scaling factors in relation to the shifts occurring during 2018-2020. Likewise, we also observe asymmetry in the stocks dataset, with more emphasis on the downstream. This gives us an idea of which part of the model fails and how the shift is affecting it (Figure˜13).

Refer to caption
(a) Values of (a,b)(a,b) over time for the automobile dataset
Refer to caption
(b) Values of (a,b)(a,b) over time for the stocks dataset
Figure 13: A visualization of how (a,b)(a,b) change as shocks occur in real-world data

D.3.6 Without quantile level adjustments

In this experiment, we verify that the ct,dtc_{t},d_{t} heuristic update step of the algorithm is necessary for improved performance. Using the experimental setup with upstream and downstream shocks as Section˜7.2, we consider the case in which we fix ct,dtc_{t},d_{t} to initial values 0.01 t\forall t, in Figure˜14.

Refer to caption
(a) Width for fixed c,dc,d
Refer to caption
(b) Coverage for fixed c,dc,d
Figure 14: A visualization of how freezing ct,dtc_{t},d_{t} to fixed values with no updates affects performance. We use α=0.1\alpha=0.1, δ=0.1\delta=0.1, γ=0.01\gamma=0.01, η=0.00\eta=0.00, k=100k=100. The first shock leads to abstention

We observe that fixing c,dc,d makes the algorithm unable to adapt properly to shifts as it produces Λval\Lambda_{\text{val}} sets frequently after the initial upstream shock, which lowers average coverage (considering that to be a miscoverage). Thus, we see that adjusting ct,dtc_{t},d_{t} allows us to mitigate abstention, however, if desiring more sensitivity to shifts, keeping them fixed at larger values could serve as a diagnostic for retraining, similar to the non-adaptive case (Table˜1).

D.3.7 Covariate shift

We consider covariate shift, in which the distribution of ww suddenly changes, at t=100,500,900t=100,500,900 which affects the entire pipeline. This shift does not have a specific upstream/downstream affect, but our method is still able to perform comparably well to the ACI, PID, and OCID methods Figure˜15. However, all the methods perform similarly and it is difficult to distinguish them.

Refer to caption
(a) width for covariate shift settings
Refer to caption
(b) Coverage for covariate shift settings
Figure 15: The distribution of ww shifts from 𝒩(0,1)\mathcal{N}(0,1) to 𝒩(3,2)\mathcal{N}(3,2), then back, then 𝒩(3,2)\mathcal{N}(-3,2). We use α=0.1\alpha=0.1, δ=0.1\delta=0.1, γ=0.01\gamma=0.01, η=0.01\eta=0.01, k=100k=100

Thus, while our method provides clear interpretability for upstream/downstream shifts, it is still able to adapt to other types of shifts without overcompensating width in comparison to other methods. However, it does not obtain any explicit advantage compared to the other methods and is still wider on average; thus we see that our method excels under asymmetric shifts at each stage, and furthermore performs best under upstream shifts.