Conformal Prediction for Causal Effects of Continuous Treatments

Thank you for your constructive feedback. We will add the corresponding improvements labeled as “Action” to our revised manuscript.

 

Answer to weaknesses

 

Empirical evaluation:

We are happy to discuss our empirical evaluation and give more background on the TCGA and MIMIC experiments below.

(i) TCGA: We synthetically model a continuous treatment based on the sum of the 10 covariates with the highest variance and assign a treatment effect that is constant in the sum of the covariates. Specifically, we model the treatment to follow a normal distribution centered at 100*sum of the 10 covariate values, and the outcome to follow a normal distribution centered at the sum of the treatment and the covariate sum times 100.

We now also obtained prediction intervals from the other baselines (MC dropout and the ensemble method) and report the coverage results below. The results are in line with our findings on datasets 1 & 2.

(ii) MIMIC: Note that numerically evaluating causal inference methods on real-world data is not possible due to the fundamental problem of causal inference, i.e., only one potential outcome (for continuous treatments out of infinitely many) is observed for each individual. Therefore, we chose to provide “insights plots” (Fig. 7) for this dataset, as (i) assessing overall coverage across the test dataset is impossible, and (ii) aggregating prediction intervals is uninformative as an evaluation method.

As the reviewer correctly noted, merely comparing the intervals based on their width is difficult, as width does not directly map to coverage. Therefore, we interpret Fig. 7 as suggestive of the coverage (see Section 5.3) and perform the semi-synthetic experiment on the TCGA dataset to be able to properly evaluate our intervals on a medical dataset.

Action: We will address the challenges of real-world evaluation and clarify the meaning of Fig. 7 in the camera-ready version of our paper.

 

Section 4:

We are happy to provide a more detailed explanation of our notation and the theoretical results, acknowledging that our paper may be notationally dense for readers unfamiliar with the related work. Specifically, we will discuss Section 4 based on (i) the reformulation as quantile regression in the high-level outline, (ii) the reformulation of Lemma 4.1 in terms of the dual optimization problem in Theorem 4.2 (and similarly the formulation of Eq. $P_S$ in terms of Theorem 4.5), and (iii) the challenges when employing our approach.

(i) Quantile regression Eq. 34:

Eq. 3 & 4 represent the general form of a quantile regression as a minimization problem based on the so-called pinnball-loss $l_{\alpha}$, a standard result in the statistical regression literature (e.g. [1],[2]). Observe that $l_{\alpha}$ can equivalently be represented as $l_{\alpha}(\theta, S) = (\alpha - \mathbb{1}_\{[\theta-S < 0]\})(\theta-S)$. This asymmetric linear loss penalizes under-prediction $(\theta-S) > 0$ proportionally to $\alpha$ and over-prediction $(\theta-S) < 0$ proportionally to $1-\alpha$, which results in the fitted function (see Eq. 3) targeting the $\alpha$-th conditional quantile.

(ii) Dual optimization Theorem 4.2 with primal optimization problem in Lemma 4.1:

Note that directly solving the problem given in Lemma 4.1 is not feasible, as it would require solving Eq. 5 for all possible imputed values $S \in \mathbb{R}$ of $S_{n+1}$. It is unclear how to feasibly find the imputation of the “true” $S_{n+1}$ needed for the interval construction in Eq. 6. We thus exploit the efficient computation through dual optimization.

We can rewrite the primal problem in Lemma 4.1 as

$$\min_{\theta > 0} \quad \sum_{i=m+1}^{n+1} (1-\alpha)u_i + \alpha v_i\\ \qquad \text{ s.t. } S_i - \theta \frac{\pi(a_i + \Delta_A \mid x_i)}{\pi(a_i \mid x_i)}) -u_i + v_i = 0, \quad \forall i=m+1\ldots,n+1, S_{n+1} = S;\\ \quad u_i, v_i \geq 0, \quad \forall i=m+1\ldots,n+1.$$

Note that this is precisely the formulation of the respective primal problem stated in Eq. P_S as well as Theorem 4.5, where $u_i$ and $v_i$ represent the primal optimization variables.

The Lagrangian of the primal problem with Lagrangian multipliers $\eta, \gamma_1, \gamma_2$ states

$$\mathcal{L} = \sum_{i=m+1}^{n+1} (1-\alpha) u_i + \alpha v_i + \sum_{i=m+1}^{n+1} \eta_i \left( S_i - \theta \frac{\pi(a_i + \Delta_A \mid x_i)}{\pi(a_i \mid x_i)} - u_i+v_i \right) - \sum_{i=m+1}^{n+1} (\gamma_{1_i}u_i + \gamma_{2_i}v_i).$$

leading to the dual problem

$$\max_{\eta_i, i=m+1,\ldots,n+1} \min_{\theta > 0} \quad \sum_{i=m+1}^{n} \eta_i \left(S_i - \theta \frac{\pi(a_i + \Delta_A \mid x_i)}{\pi(a_i \mid x_i)} \right) + \eta_{n+1} \left( S - \theta \frac{\pi(a^{\ast} \mid x_{n+1})}{\pi(a_{n+1} \mid x_{n+1})} \right) \qquad \text{ s.t.} -\alpha \leq \eta_i \leq 1-\alpha, \quad \forall i=m+1,\ldots,n+1,$$

with dual optimization variables $\eta_{m+1},...,\eta_{n+1}$.

For more background on dual (convex) optimization theory, we refer to, e.g., [3], [4]. Specific details on our derivations can be found in our proofs of Theorems 4.2 and 4.5 in Appendix D. For more details on the coverage guarantees/reformulation from Eq. 6 to Eq. 9, we refer to the respective part of our proof of Thm. 4.2 in L 663ff in our manuscript.

(iii) Challenges:

The main challenge of our approach is the need for sufficient prior knowledge to set a bound $M$ on the propensity estimation error (see “answer to limitations” below). However, this bound is far less restrictive than the assumption of a known propensity score, and is common in the literature (see Section 2). In practice, one aims for valid but also narrow intervals. If the calibration dataset is not representative of the full population, the resulting intervals are potentially wider than necessary for test samples in segments of the output space with few calibration data. For transparency, we discuss all limitations in Section 6 of our manuscript.

Action: We will provide more background and explain the intuition for our derivations, as well as the challenges of our approach, in more detail in the camera-ready version.

 

Answer to questions: Numerical stability

We are happy to discuss the numerical stability of our proposed optimization approach. In Scenario 1, the only source of potential instability can be a very low propensity score in low-overlap regions of the covariate space. This, however, is only a problem if $\pi(a|x) <<\pi(a+\Delta|x)$. We can thus consider this unlikely in practice. For example, consider a patient who would be treated with a dosage of 10mg of some medication, which is prescribed in the range from 0 to 50mg. A practitioner is likely to be interested in the effect of an increase of the dosage to 15mg (i.e., $\Delta=5$) in contrast to an increase to 50mg. Furthermore, it is reasonable to assume the propensity function to be locally smooth. Therefore, the instability of $\pi(a|x) <<\pi(a+\Delta|x)$ is unlikely to occur.

In Scenario 2, we could additionally face an underflow of $\exp(-\frac{(a_i-a)^2}{2\sigma^2})$ or a blow-up of the pre-factor $(\sigma\hat{\pi}(x_i))^{-1}$. As a remedy, we can reparameterize the problem to work in the log domain and only exponentiate after subtracting a stable (maximum) constant across all datapoints. This additionally prevents the Jacobian/Hessian of the constraints from becoming nearly singular or wildly varying, causing Newton‐type steps to blow up or stall.

Action: We will add a discussion on the stability of the optimization approaches to our manuscript.

 

Answer to limitations: Considerations for practical application

We are happy to discuss the requirements of applying our method in practice. First, our method is designed for single continuous treatments with a univariate outcome, which is common in dosaging (e.g., determining the dosage of chemotherapy/ insulin/ hypertension drugs, etc).

• In randomized controlled trials (when the propensity score is known), we can directly apply Theorem 4.2 on top of the trained outcome prediction model $\phi$ to construct our CP intervals at target coverage level $\alpha$.
• In observational studies, when the propensity score is unknown, we require the practitioner to have sufficient prior knowledge to set a bound $M$ on the propensity estimation error. The practitioner should choose this bound with care and rather in a conservative way to prevent undercoverage. Then, the CP intervals can be derived based on Theorem 4.5.

Our intervals do not suffer from undercoverage due to limited data. Importantly, CP guarantees are valid for any sample size and our method is kept fully model agnostic, enabling the practitioner to make reliable judgements based on limited data and an ML model of choice.

Action: We will discuss practical considerations of applying our method in medicine.

 

[1] Koenker, et al. Regression quantiles. Econometrica (1978): 33-50.

[2] Koenker. Quantile regression. Vol. 38. Cambridge University Press, 2005.

[3] Boyd, et al. Convex optimization. Cambridge University Press, 2004.

[4] Luenberger, et al. Linear and nonlinear programming. Vol. 2. Reading, MA: Addison-Wesley, 1984.

Thank you for your constructive feedback and the positive evaluation of our manuscript. Below, we provide answers to all questions. We will add the corresponding improvements labeled as “Action” to our revised manuscript.

 

Answers to weaknesses:

 

Causal identifiability:

Throughout our work, we make the three standard assumptions in causal inference: positivity, consistency, and unconfoundedness (see Section 3). These assumptions are standard in the literature (e.g., [1],[2],[3],[4]). Under these assumptions, the conditional average treatment effect $\tau(x) = \mathbb{E}[Y(a)-Y(a^{‘})| X=x]$ is identified as $\tau(x) = \mathbb{E}[Y|A=a, X=x] -\mathbb{E}[Y|A=a^{‘},X=x]$. The potential outcome estimation thus simplifies to estimating the conditional expectation; in our case, estimated through the ML model $\phi$ (see Section 3). In the main paper, we provide conformal intervals for potential outcomes. We derive CP intervals for the treatment effect in Appendix A. Given that these assumptions are standard and generally satisfied by design (e.g., in clinical trials) (e.g., [4],[5],[6]), we do not share the concern raised by the reviewer. However, we are happy to clarify further aspects if needed.

 

Causal UQ baselines:

To foster the evaluation of our method, we additionally provide coverage results on two further causal UQ baselines: (i) causal Gaussian processes, and (ii) the causal CP method for binary treatments by Lei & Candés.

(i) The GP is only able to capture the true potential outcome in the prediction intervals for small distribution shifts (Δ = 1) on dataset 2. However, the empirical coverage is extremely low: $\alpha = 0.05: 0.125; \alpha = 0.1: 0.125; \alpha = 0.2: 0.0833$. This result aligns with our expectations, as the aleatoric uncertainty in our experiments is low. Therefore, the GP intervals will be small (average width of 0.1293) and barely provide valid intervals for the potential outcomes after the intervention.

(ii) We compare our method and the method by Lei and Candés on the semi-synthetic TCGA dataset. However, when comparing the interval width, we see that our method is superior: Our intervals have an average width of 0.6255 (sd = 0.1714). In contrast, the intervals on the binarized treatment via the Lei-Candés method have an average width of 3.2876 (sd = 0.5587). Hence, the intervals by our method are – by far – more informative.

 

Dependence on propensity score estimation, model misspecification, or unmeasured confounding:

Propensity estimation: Our intervals depend on the quality of propensity estimation through the parameter $M$. Larger $M$ represents more uncertainty in the propensity estimation and thus wider intervals. We conducted an ablation study on the TCGA dataset: the results confirm this pattern, as expected.

Model misspecification: Our CP intervals are exactly designed to account for the bias and uncertainty stemming from model misspecification. Misspecified models lead, on average, to higher non-conformity scores on the calibration dataset and, thus, to wider intervals.

Unmeasured confounding: We highlight again that we make the three standard assumptions in causal inference, including unconfoundedness. In the presence of unmeasured confounders, the potential outcomes are not (point-)identified. Our method is not designed for this setting. For future work, extending our method to partial identification of potential outcomes based on sensitivity analysis could be of interest.

Action: We will add the above discussion to the camera-ready version of our paper, including the results of the new ablation study.

 

Answers to questions:

 

(1) Causal UQ baselines and related work:

For your question on causal UQ baselines, please see our answer above.

There exists a variety of work for uncertainty quantification (UQ) in causal effect estimation for continuous treatments. Existing methods for UQ of causal quantities are often based on Bayesian methods (e.g., [1],[2],[7]). However, Bayesian methods require the specification of a prior distribution based on domain knowledge and are thus neither robust to model misspecification nor generalizable to model-agnostic machine learning models. Other methods only provide asymptotic guarantees (e.g. [3],[8]). The strength of conformal prediction, however, is to provide finite-sample uncertainty guarantees. Overall, none of the methods tackles the problem of distribution-free uncertainty quantification for causal effect estimation of continuous treatments in finite sample settings.

For an extensive discussion of related work on CP for causal quantities, we refer the reviewer to our extensive extended literature review in Appendix E.

 

(2) Exchangeability violation:

Coverage guarantees of existing CP intervals essentially rely on the exchangeability of the non-conformity scores. Exchangeability assures that the nonconformity score of the test point $n+1$ is equally likely to fall anywhere among the calibration scores, its rank is uniform, and that uniformity is exactly what yields the distribution‑free coverage guarantee. Without exchangeability, the rank is not guaranteed to be uniform, and the marginal coverage bound can fail.

However, intervening on treatment $A$ shifts the propensity function and, therefore, induces a shift in the covariates between calibration and test data, specifically in the treatment $A$. Therefore, exchangeability is not fulfilled, and the coverage guarantees might fail (at least for existing methods).

As a remedy, we present a novel and powerful remedy in our work: The overall distribution of the confounders $X$ is assumed to stay constant between train, calibration, and test data (as standard in ML problems). This is completely orthogonal to constructing intervals for different (e.g., young or old) patients. Note that CP intervals are constructed for only one sample/patient at a time. This means that different intervals are constructed for patients with different features $X$. In other words, the intervals are constructed conditionally on $X$, but the coverage guarantee is marginal across the complete population of $X$. Overall, the shift from one patient to another does not pose any challenges for CP methods.

Action: We will include a discussion on the role of the exchangeability assumption and the limitations it poses for existing causal CP methods in the camera-ready version of our paper. We will clarify how our proposed method circumvents these limitations, offering a more robust and practical approach under realistic causal settings.

 

(3) Aleatoric and epistemic uncertainty:

This is an interesting question! First, we like to highlight that CP guarantees marginal (total) coverage, agnostic to the split between aleatoric and epistemic error.

Aleatoric uncertainty denotes the irreducible noise inherent in the data–generation process, such as measurement error or intra‑subject variability. Epistemic uncertainty is present due to the limited data available and due to model misspecification. This uncertainty is, in principle, reducible by gathering more data or improving the model.

 

[1] A. Alaa and M. van der Schaar. Bayesian inference of individualized treatment effects using multi-task Gaussian processes. In Conference on Neural Information Processing Systems (NeurIPS), 2017

[2] K. Hess et. al. Bayesian neural controlled differential equations for treatment effect estimation. In International Conference on Learning Representations (ICLR), 2024.

[3] J. Jonkers et.al. Conformal Monte Carlo meta-learners for predictive inference of individual treatment effects. arXiv preprint, 2024.

[4] P. R. Rosenbaum and D. B. Rubin.The central role of the propensity score in observational studies for causal effects. Biometrika, 70(1):41–55, 1983

[5] S. Feuerriegel et. al. Causal machine learning for predicting treatment outcomes. Nature Medicine, 30(4):958–968, 2024.

[6] P. Sanchez et al. Causal machine learning for healthcare and precision medicine. Royal Society Open Science 9(8), 2022

[7] A. Jesson et.al. Identifying causal-effect inference failure with uncertainty-aware models. In Conference on Neural Information Processing Systems (NeurIPS), 2020.

[8] Y. Jin et.al. Sensitivity analysis of individual treatment effects: A robust conformal inference approach. Proceedings of the National Academy of Sciences of the United States of America, 120(6), 2023.

Thank you for your constructive feedback and the positive evaluation of our manuscript. Below, we provide answers to all questions. We will add the corresponding improvements labeled as “Action” to our revised manuscript.

 

Answers to weaknesses:

 

(1) Uncertainty quantification of propensity estimation:

Thank you for raising this question. Accounting for the propensity estimation error is exactly one of the strengths of our method (see Section 4.2) and differentiates our framework from former work in this field (see Section 2). Under Assumption 1 (by stating a user-specified upper bound on the estimation error), our intervals are valid for all misspecified propensity estimations with estimation error smaller than or equal to the bound.

 

(2) CP Baselines:

We followed your suggestion and added a comparison to the naive vanilla CP (V-CP), i.e., CP, which does not account for the distribution shift. Specifically, we now report the performance of (i) the naive vanilla CP (V-CP), i.e., CP, which does not account for the distribution shift, and (ii) the method proposed by Lei and Candés for binary treatments (for which we binarize our treatment).

(i) We observe that V-CP does not achieve any valid prediction interval, across all distribution shifts and confidence levels. This can be explained by the good prediction performance of the underlying model (see Supplement G). Thus, V-CP intervals are extremely small (average width of 0.0003) and can never cover the true potential outcome after the intervention. Overall, the results confirm the importance of accounting for the distribution shift induced by the intervention.

(ii) We compare our method and the method by Lei and Candés on the semi-synthetic TCGA dataset, where both methods fulfill the coverage guarantees. However, when comparing the interval width, we see that *our method is superior: Our intervals have an average width of 0.6255 (sd = 0.1714). In contrast, the intervals on the binarized treatment via the Lei-Candés method have an average width of 3.2876 (sd = 0.5587). Hence, the intervals by our method are – by far – more informative.

Action: We will include our performance results of the V-CP and the binary method by Lei and Candés in the camera-ready version of our manuscript.

 

(3) Functional form of $\delta$:

We are happy to explain our choice of representation of the Dirac delta function. Note that it is a common way to represent the Dirac delta as the limit of a family of probability densities that become infinitely concentrated at zero. One such family is the Gaussian density. Other options are, e.g., the representation through Lorentz spikes, rectangular pulses, or sinc kernels. For background on the representation, we refer the reviewer to, e.g., [1]. Please let us know if there are other functional forms in our paper that should be explained in more detail.

 

Answers to questions:

 

(1) Motivation for UQ:

The reviewer is correct that heterogeneous effects, such as CATE, are designed to capture variability in the treatment effect estimate; however, this only refers to the variability based on the heterogeneity in effects for varying covariates. It does not include the variability in the effects due to noise in the data or modeling error (which is our objective). Therefore, uncertainty quantification (as in our paper) is necessary in causal inference settings in the same way as for other ML applications.

 

(2) Exchangeability assumption:

Coverage guarantees of existing CP intervals essentially rely on the exchangeability of the non-conformity scores. Exchangeability assures that the nonconformity score of the test point $n+1$ is equally likely to fall anywhere among the calibration scores, its rank is uniform, and that uniformity is exactly what yields the distribution‑free coverage guarantee. Without exchangeability, the rank is not guaranteed to be uniform, and the marginal coverage bound can fail.

However, intervening on treatment $A$ shifts the propensity function and, therefore, induces a shift in the covariates between calibration and test data, specifically in treatment $A$. Therefore, exchangeability is not fulfilled, and the coverage guarantees might fail.

As a remedy, we present a novel and powerful remedy in our work: our calibration step differs from the standard procedure in that we conditionally calibrate the non-conformity scores depending on the tilting function $f$ to achieve marginal coverage for the interventional and thus shifted data. Thereby, we ensure that the exchangability assumption holds.

 

(3) Quantile regression Eq. 3 & 4:

Equations 3 & 4 represent the general form of a quantile regression as a minimization problem based on the so-called pinnball-loss $l_{\alpha}$, a standard result in the statistical regression literature (e.g. [2],[3],[4],[5]).

Observe, that $l_{\alpha}$ can equivalently be represented as $l_{\alpha}(\theta, S) = (\alpha - \mathbb{1}_\{[\theta-S < 0]\})(\theta-S)$. This asymmetric linear loss penalizes under-prediction $(\theta-S) > 0$ proportionally to $\alpha$ and over-prediction $(\theta-S) < 0$ proportionally to $1-\alpha$, which results in the fitted function (see Eq. 3) targeting the $\alpha$-th conditional quantile.

Action: We will include more background information on the reformulation of quantile regression in the camera-ready version of our manuscript.

 

(4) Dimension of vector in line 208:

Thank you for noting the typo! Initially, we referred to the sizes of both the training and calibration datasets as $n$, resulting in the dimension stated in line 208; however, we later adapted our formulation to enhance the reader’s understanding. The correct dimension of $\eta$ is $n-m$.

Action: We will correct the typo in the camera-ready version of our paper.

 

(5), (6) & (7) Dual optimization:

We are happy to explain the necessity and the derivations of the dual optimization problem in Theorem 4.2 from the primal optimization problem in Lemma 4.1.

First, observe that the objective function and loss formulation in Eq. 3 & 4 denote a general quantile regression, which does not directly correspond to our CP problem. The precise (primal) optimization problem for our CP setting is given in Lemma 4.1 (Eq. 5).

Note that directly solving the problem given in Lemma 4.1 is not feasible, as it would require solving Eq. 5 for all possible imputed values $S \in \mathbb{R}$ of $S_{n+1}$. It is unclear how to feasibly find the imputation of the “true” $S_{n+1}$ needed for the interval construction in Eq. 6. We thus follow former results in the literature (e.g., [6]) and exploit the efficient computation through dual optimization.

We can rewrite the primal problem in Lemma 4.1 as

$\min_{\theta > 0} \quad \sum_{i=m+1}^{n+1} (1-\alpha)u_i + \alpha v_i\\ \qquad \text{ s.t. } S_i - \theta \frac{\pi(a_i + \Delta_A \mid x_i)}{\pi(a_i \mid x_i)}) -u_i + v_i = 0, \quad \forall i=m+1\ldots,n+1, S_{n+1} = S, \quad u_i, v_i \geq 0, \quad \forall i=m+1\ldots,n+1.$

Note that this is precisely the formulation of the respective primal problem stated in Eq. P_S as well as Theorem 4.5, where $u_i$ and $v_i$ represent the primal optimization variables.

The Lagrangian of the primal problem with Lagrangian multipliers $\eta, \gamma_1, \gamma_2$ states $\mathcal{L} = \sum_{i=m+1}^{n+1} (1-\alpha) u_i + \alpha v_i + \sum_{i=m+1}^{n+1} \eta_i \left( S_i - \theta \frac{\pi(a_i + \Delta_A \mid x_i)}{\pi(a_i \mid x_i)} - u_i+v_i \right) - \sum_{i=m+1}^{n+1} (\gamma_{1_i}u_i + \gamma_{2_i}v_i).$

leading to the dual problem

$\max_{\eta_i, i=m+1,\ldots,n+1} \min_{\theta > 0} \quad \sum_{i=m+1}^{n} \eta_i \left(S_i - \theta \frac{\pi(a_i + \Delta_A \mid x_i)}{\pi(a_i \mid x_i)} \right) + \eta_{n+1} \left( S - \theta \frac{\pi(a^{\ast} \mid x_{n+1})}{\pi(a_{n+1} \mid x_{n+1})} \right) \\ \qquad \text{ s.t.} -\alpha \leq \eta_i \leq 1-\alpha, \quad \forall i=m+1,\ldots,n+1,$

with dual optimization variables $\eta_{m+1},...,\eta{n+1}$.

For more background on dual (convex) optimization theory, we refer the reviewer to, e.g., [7], [8], due to space constraints in the rebuttal. Specific details on our derivations can be found in the proofs of Theorems 4.2 and 4.5 in Appendix D.

For more details on the coverage guarantees/reformulation from Eq. 6 to Eq. 9, we refer the reviewer to the respective part of the proof of Thm. 4.2 in lines 663 and following in our manuscript.

Action: We will provide more background and explain the intuition for the dual optimization step in more detail in the camera-ready version of our manuscript.

 

[1] Gel’fand, I. M. and Shilov, G. E. . Generalized Functions, Vol. 1: Properties and Operations (Academic Press, 1964)

[2] Koenker, R., and Bassett Jr., G. "Regression quantiles." Econometrica: journal of the Econometric Society (1978): 33-50.

[3] Koenker, R. Quantile regression. Vol. 38. Cambridge University Press, 2005.

[4] Yu, K., and Moyeed, R. A. "Bayesian quantile regression." Statistics & Probability Letters 54.4 (2001): 437-447.

[5] Belloni, A., and Chernozhukov, V. "ℓ 1-penalized quantile regression in high-dimensional sparse models." (2011): 82-130.

[6] Gibbs, I. et al. Conformal prediction with conditional guarantees. arXiv preprint, 2023.

[7] Boyd, S. P., and Vandenberghe. L.. Convex optimization. Cambridge University Press, 2004.

[8] Luenberger, D. G., and Ye, Y. Linear and nonlinear programming. Vol. 2. Reading, MA: Addison-Wesley, 1984.

Thank you for your feedback and the positive evaluation of our manuscript. Below, we provide answers to all the questions and concerns you raised. We will add the corresponding improvements labeled as “Action” to our revised manuscript.

 

Comparison to naive CP baseline:

We followed your suggestion and added a comparison to the naive vanilla CP (V-CP), i.e., CP, which does not account for the distribution shift. This may be of interest to the reader to understand the importance of the proposed approach. In our experiments, we observed that V-CP does not provide any valid prediction interval, across all distribution shifts and confidence levels. This can be explained by the good prediction performance of the underlying model (see Supplement G). Thus, V-CP intervals are extremely small (average width of 0.0003) and can never cover the true potential outcome after the intervention. Overall, the results confirm the importance of accounting for the distribution shift induced by the intervention.

To investigate whether simple weighted CP methods for causal tasks on binary treatments sufficiently address this challenge, we compare our method with the method by Lei and Candés. Of note, both methods fulfill the coverage guarantees. However, when comparing the interval width, we see that our method is superior: Our intervals have an average width of 0.6255 (sd = 0.1714). In contrast, the intervals on the binarized treatment via the Lei-Candés method have an average width of 3.2876 (sd = 0.5587). Hence, the intervals by our method are – by far – more informative.

Action: We will include our observations on the performance of both methods in the camera-ready version of our manuscript.

 

Computational complexity:

You are correct that the overall complexity is $O(\log(\frac{S_{up} - S_{low}}{\epsilon} \cdot \sigma_s(n_c))$ following from the argumentation in the derivation. Please excuse the typo.

We acknowledge that our method is computationally rather complex, which we also stated when discussing the limitations of our method. However, this is not unique to our method, but also inherent to the baselines. As a case in point, the complexity of deriving intervals through the baselines (i.e., MC-dropout or ensemble methods) depends on the number of MC samples or models, respectively, both of which must be rather large to offer similar computational performance. In particular, both baselines do not offer valid prediction intervals, which is the contribution of our method. On top of that, the latter baseline (ensemble methods) scales with the precision of the intervals, which might be difficult to control and is thus impractical in real-world applications. In our work, we use optimization as a tool to provide conformal prediction intervals. Future research should focus on developing more efficient optimization algorithms for this task.

Action: We will correct the typo in the camera-ready version of our manuscript.

 

Runtime on MIMIC dataset:

This is an interesting question. Due to the optimization problem that needs to be solved in our approach, the time needed for computing our CP intervals is naturally higher than the time needed to calculate the MC intervals. While the theoretical runtime of CP may exhibit more complex scaling behavior (e.g., cubic or non-linear), our empirical results demonstrate that CP scales well in practice. Specifically, computing CP intervals takes roughly 10 times longer than computing MC intervals. However, we emphasize that MC intervals are generally not faithful and therefore not directly comparable. To better understand the scaling behavior of our CP method in practice, we examine computation time relative to the average interval width. Under this lens, both methods exhibit a similar time-to-width ratio of approximately one. In other words, our experiments do not indicate that CP intervals scale poorly—in fact, we observe good scaling properties in practice.

Action: We will include a full runtime comparison on the MIMIC dataset in the camera-ready version of our paper.