Conformal prediction for causal effects of continuous treatments

## Abstract

Abstract Uncertainty quantification of causal effects is crucial for safety-critical applications such as personalized medicine. A powerful approach for this is conformal prediction, which has several practical benefits due to model-agnostic finite-sample guarantees. Yet, existing methods for conformal prediction of causal effects are limited to binary/discrete treatments and make highly restrictive assumptions such as known propensity scores. In this work, we provide a novel conformal prediction method for potential outcomes of continuous treatments. We account for the additional uncertainty introduced through propensity estimation so that our conformal prediction intervals are valid even if the propensity score is unknown. Our contributions are three-fold: (1) We derive finite-sample prediction intervals for potential outcomes of continuous treatments. (2) We provide an algorithm for calculating the derived intervals. (3) We demonstrate the effectiveness of the conformal prediction intervals in experiments on synthetic and real-world datasets. To the best of our knowledge, we are the first to propose conformal prediction for continuous treatments when the propensity score is unknown and must be estimated from data.

## 1 Introduction

Machine learning (ML) for estimating causal quantities such as causal effects and the potential outcomes of treatments is nowadays widely used in real-world applications such as personalized medicine (Feuerriegel et al., ). However, existing methods from causal ML typically focus on point estimates (e.g., Nie et al., ; Patrick Schwab et al., ), which means that the uncertainty in the predictions is neglected and hinders the use of causal ML in safety-critical applications (Feuerriegel et al., ; Kneib et al., ).
As the following example shows, uncertainty quantification (UQ) of causal quantities is crucial for reliable decision-making.

*Motivating example:* Let us consider a doctor who seeks to determine the dosage of chemotherapy in cancer care. This requires estimating the tumor size in response to the dosage for a specific patient profile. A point estimate will predict the *average* size of the tumor post-treatment, but it will neglect that chemotherapy is ineffective for some patients. In contrast, UQ will give a *range* of the tumor size that is to be expected post-treatment, so that doctors can assess the probability that the patients will actually benefit from treatment. This helps to understand the risk of a treatment being ineffective and can guide doctors to choose treatments that are effective *with large probability*.

A powerful method for UQ is *conformal prediction* (CP) (Lei & Wasserman, ; Papadopoulos, ; Vovk et al., ). CP provides model-agnostic and distribution-free, finite-sample predictive intervals for quantifying uncertainty. CP has been widely used for traditional, predictive ML (e.g., Angelopoulos et al., ; Barber et al., ; Gibbs et al., ), where it has been shown to yield reliable prediction intervals in finite-sample settings (see Fig. 1). Recently, there have been works that adapt CP for estimating causal quantities (see Fig. 2 for an overview). Yet, existing methods for CP focus on binary or discrete treatments (e.g., Alaa et al., ; Jonkers et al., ; Lei & Candès, ), but *not* continuous treatments, which is our novelty.

Figure: Figure 1: CP intervals on finite-sample data. UQ methods with asymptotic guarantees might suffer from under-coverage and are often *not* faithful. Thus, we aim at CP with finite-sample guarantees.
Refer to caption: x1.png

However, adapting CP to causal quantities with continuous treatments is *non-trivial* for two main reasons. Challenge  a: Intervening on the treatment induces a shift in the covariate distribution, specifically, in the propensity score. As a result, the so-called *exchangeability assumption*, which is inherent to CP (Vovk et al., ), is violated between the observational and interventional distribution, and because of this, standard CP intervals are not valid. Thus, we must later account for the distribution shift and derive *treatment-conditional guarantees*. Challenge  b: Assessing the aforementioned shift in the distribution requires information about the propensity score; yet, the propensity score is typically unknown. Hence, estimating the propensity score introduces *additional uncertainty*. However, incorporating the additional uncertainty in the overall CP intervals cannot be done in a simple plug-in manner, and it is highly non-trivial.

In this paper, we develop a novel CP method for causal quantities, such as causal effects and potential outcomes of continuous treatments. Our method is carefully designed to account for the additional uncertainty introduced during propensity estimation and is therefore applicable to settings where the propensity score is known or even unknown.

Our contributions:(^1^11Code and data are available at our public GitHub repository: [https://github.com/m-schroder/ContinuousCausalCP](https://github.com/m-schroder/ContinuousCausalCP)) (1) We propose a novel method for CP of causal quantities such as potential outcomes or treatment effects of continuous treatments. For this, we mathematically derive finite-sample prediction intervals for potential outcomes under known and unknown propensity functions. (2) We provide an algorithm for efficiently calculating the derived intervals. (3) We demonstrate the effectiveness of the derived CP intervals in experiments on multiple datasets.

## 2 Related Work

UQ for causal effects: Existing methods for UQ of causal quantities are often based on Bayesian methods (e.g., Alaa & van der Schaar, ; Hess et al., ; Hill, ; Jesson et al., ). 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. A common ad hoc method for computing uncertainty intervals is Monte Carlo (MC) dropout (Gal & Ghahramani, ). However, MC dropout yields approximations of the posterior distribution, which are *not* faithful (Le Folgoc et al., ).

Conformal prediction: CP (Lei & Wasserman, ; Papadopoulos, ; Vovk et al., ) has recently received large attention for finite-sample UQ. For a prediction model $\phi$ trained on dataset $D_{T}=(X_{i},Y_{i})_{i=1,\ldots,m}$ and a new test sample $X_{k}$, CP aims to construct a prediction interval $C(X_{k})$ such that $P(Y_{k}\in C(X_{k}))\geq 1-\alpha$ for some significance level $\alpha$. We refer to (Angelopoulos et al., ) for an in-depth overview. Due to its strong finite-sample mathematical guarantees, CP is widely used for traditional, predictive ML with widespread applications such as in medical settings (Zhan et al., ) or drug discovery (Alvarsson et al., ; Eklund et al., ).

Several extensions have been developed for CP. One literature stream focuses on CP with *marginal coverage* under distribution shifts between training and test data (e.g., Cauchois et al., ; Fannjiang et al., ; Gendler et al., ; Ghosh et al., ; Gibbs & Candès, ; Gibbs et al., ; Guan, ; Lei & Candès, ; Podkopaev & Ramdas, ; Tibshirani et al., ; Yang et al., ). Our setting later also involves a distribution shift due to the intervention on the treatment but differs from the latter in that the true distribution shift is unknown. Another literature stream constructs intervals *conditional* on the variables following the shifted distribution. Since, in general, exact conditional coverage has been proven impossible (Lei & Wasserman, ; Vovk, ), the works in this literature stream have two key limitations: (1) they only guarantee *approximate* conditional coverage (e.g., Barber et al., ; Cai et al., ; Lei & Wasserman, ; Romano et al., ); or (2) they are restricted to specific data structures such as binary variables (e.g., Lei & Wasserman, ; Vovk, ). Because of that, none of the existing methods for marginal and conditional coverage can be applied to derive *exact* finite-sample prediction intervals for causal quantities of continuous treatments.

Conformal prediction for causal quantities: Only a few works focus on CP for causal quantities (see Fig. 2). Examples are methods aimed at off-policy learning (Taufiq et al., ; Zhang et al., ), conformal sensitivity analysis (Yin et al., ), or meta-learners for the conditional average treatment effect (CATE) (Alaa et al., ; Jonkers et al., ; Wang et al., ). However, there are crucial differences to our setting: First, the existing works (a) assume that the propensity is *known* and thus achieve exact finite-sample coverage guarantees, or the existing works (b) focus on the easier task of giving *asymptotic* guarantees but then might suffer from under-coverage because of which the intervals are *not* faithful. On top of that, the existing CP methods are designed for *binary* or *discrete* treatments. Applying such methods to discretized continuous treatments leads to ill-defined causal estimands.
Therefore, none of the existing methods are applicable to our continuous treatment setting. We offer a detailed discussion in Supplement D.

Figure: Figure 2: Key works on CP in causal inference.
Refer to caption: x2.png

Research gap: To the best of our knowledge, no work has provided exact finite-sample prediction intervals for causal quantities of continuous treatments, which is our novelty.

## 3 Problem formulation

Notation:
We denote random variables by capital letters $X$ with realizations $x$. Let $P_{X}$ be the probability distribution over $X$. We omit the subscript whenever it is obvious from the context. For discrete $X$, we denote the probability mass function by $P(x)=P(X=x)$ and the conditional probability mass functions by $P(y\mid x)=P(Y=y\mid X=x)$ for a discrete random variable $Y$. For continuous $X$, $p(x)$ is the probability density function w.r.t. the Lebesgue measure.

Setting: Let the data $(X_{i},A_{i},Y_{i})_{i=1,\ldots,n}$ consisting of observed confounders $X\in\mathcal{X}$, a continuous treatment $A\in\mathcal{A}$, and an outcome $Y\in\mathcal{Y}$ be drawn *exchangeably* from the joint distribution $P$. Additionally, let a new sample of confounders $X_{n+1}$ be drawn independently from the marginal distribution $P_{X}$. Throughout our work, we split the dataset into a proper training dataset $D_{T}=(X_{i},A_{i},Y_{i})_{i=1,\ldots,m}$, and a calibration dataset $D_{C}=(X_{i},A_{i},Y_{i})_{i=m+1,\ldots,n}$. Furthermore, let $\pi(a\mid x)$ define the generalized propensity score for treatment $A=a$ given $X=x$.

Throughout this work, we build upon the potential outcomes framework (Rubin, ). We denote the potential outcomes after a hard intervention $a^{*}$ by $Y(a^{*})$ and after a soft intervention $A^{*}(x)\sim\tilde{\pi}(a\mid x)=P_{A^{*}\mid X=x}$ by $Y(A^{*}(x))$.(^2^22Interventions are characterized by two classes: hard (structural) and soft (parametric) interventions. Hard interventions directly affect the treatment by setting it to a specific value and removing the edge in the graph (as in the do-operator). Soft interventions do not change the structure of the graph but affect the conditional distribution of the treatment given the confounders, i.e., the propensity score.) We make three standard identifiability assumptions for causal effect estimation: positivity, consistency, and unconfoundedness (e.g., Alaa et al., ; Jonkers et al., ). Finally, we consider an arbitrary machine learning model $\phi$ to predict the potential outcomes. Hence, we define the outcome prediction function as $\phi:\mathcal{X}\times\mathcal{A}\rightarrow\mathbb{R}$, $\phi(X,A)\mapsto Y$.

Our objective: In this work, we aim to derive *conformal prediction intervals* $C(X_{n+1},\Diamond)$ for the prediction of a potential outcome $Y_{n+1}(\Diamond)$ of a new data point under either hard, $\Diamond=a^{*}$, or soft intervention, $\Diamond=A^{*}(X_{n+1})\sim\tilde{\pi}(a\mid X_{n+1})$. The derived intervals are called *valid* for any new exchangeable sample $X_{n+1}$ with non-exchangeable intervention $\Diamond$, i.e.,

$$ $\displaystyle P(Y_{n+1}(\Diamond)\in C(X_{n+1},\Diamond))\geq 1-\alpha,\quad \Diamond\in\{a^{*},A^{*}(X_{n+1})\}.$ (1) $$

for some significance level $\alpha\in(0,1)$. Of note, our CP method can later be used with an arbitrary machine learning model $\phi$ for predicting the potential outcomes.

In CP, the interval $C$ is constructed based on so-called *non-conformity scores* (Vovk et al., ), which capture the performance of the prediction model $\phi$. For example, a common choice for the non-conformity score is the residual of the fitted model $s(X,A,Y)=\left|Y-\phi(X,A)\right|$, which we will use throughout our work. For ease of notation, we define $S_{i}:=s(X_{i},A_{i},Y_{i})$.

Why is CP for causal quantities non-trivial? There are two main reasons. First, coverage guarantees of CP intervals essentially rely on the exchangeability of the non-conformity scores. However, intervening on treatment $A$ shifts the propensity function and, therefore, induces a shift in the covariates $(X,A)$ ($\rightarrow$ Challenge  a). Formally, we have a *propensity shift* in which the intervention $\Diamond$ shifts the propensity function $\pi(a\mid x)$ to either a Dirac-delta distribution of the hard intervention, $\delta_{a^{*}}(a)$, or to the distribution of the soft intervention, $\tilde{\pi}(a\mid x)$, without affecting the outcome function $\phi(x,a)$. As a result, the test data sample under $\Diamond$ does *not follow the same distribution* as the train and calibration data, i.e., the exchangeability assumption is violated.

Second, the propensity score $\pi$ is commonly *unknown* in observational data and, therefore, must be estimated, which introduces additional uncertainty that one must account for when constructing CP intervals ($\rightarrow$ Challenge  b). Crucially, existing coverage guarantees (e.g., Vovk, ; Tibshirani et al., ) do *not* hold in our setting. Instead, we must derive new intervals with valid *coverage under propensity shift*.

In the following section, we address the above propensity shift by performing a calibration conditional on the propensity shift induced by the intervention, which allows us then to yield valid prediction intervals with significance level $(1-\alpha)$ for potential outcomes of a specific hard or soft intervention. We derive the potential outcomes and emphasize that the extension to causal effects is straightforward. For the latter, one combines the intervals for each potential outcome under a certain treatment and without treatment, so that eventually arrives at CP intervals for the individual treatment effect (ITE). Details are Supplement A.

## 4 CP intervals for potential outcomes of continuous treatments

Figure: Figure 3: Use cases of the two scenarios: 1 The new assignment is a function of the original policy (i.e., soft intervention). 2 The policy in the dataset is unknown. The new assignment cannot be expressed as a function of the original policy (i.e., hard intervention).
Refer to caption: x3.png

Recall that intervening on test data samples breaks the exchangeability assumption necessary for the validity, i.e., the guaranteed coverage of at least $(1-\alpha)$, of CP intervals. Therefore, we now construct CP intervals where we account for a (potentially unknown) covariate shift in the test data induced by the intervention.

Scenarios: In our derivation, we distinguish two different scenarios(see Fig. 3):

- 1
Known propensity score (see Section 4.1): If the propensity score in the observational data is known, it means that the treatment policy is known. Then, we aim to update the policy by increasing/decreasing the treatment by a value $\Delta_{A}$, i.e., $A^{*}(X)=A+\Delta_{A}$.
*Example:*
Imagine a doctor is about to prescribe a medication to a new patient. Instead of prescribing the same dosage as he would have prescribed to a similar patient in the past, the doctor is interested in the potential health outcome of the patient when increasing (or decreasing) the original dosage by amount $\Delta_{A}$.

- 2
Unknown propensity score (see Section
4.2): In observational data, the propensity score is typically unknown. Therefore, we are interested in the effect of hard interventions, i.e., $a^{\ast}$. Here, we face additional uncertainty from the propensity score estimation ($\rightarrow$  Challenge  b).
*Example:*
In our running example, a patient comes to a new doctor who has never prescribed the respective medication and thus will base the decision on observational data (electronic health records) from other physicians, yet the observational data was collected under a different, unknown treatment policy. Therefore, the prescribed intervention (and thus dosage) cannot be expressed in terms of the policy in the observational data.

In our derivations, we make use of the following two mathematical tools. First, we define the *propensity shift*. Formally, it is the shift between the observational and interventional distributions, $P$ and $\tilde{P}$, in terms of the tilting of the propensity function by a non-negative function $f$. Hence, we have

$$ $\displaystyle\tilde{\pi}(a\mid x)=\frac{f(a,x)}{\mathbb{E}_{P}[f(A,X)]}\pi(a \mid x).$ (2) $$

for some $f$ with $\mathbb{E}_{P}[f(X,A)]>0$ and $a\in\mathcal{A},x\in\mathcal{X}$.

Second, our CP method will build upon ideas from so-called split conformal prediction (Papadopoulos, ; Vovk et al., ), yet with crucial differences. In our methods, the 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.

High-level outline: Our derivation in Sections 4.1 and 4.2 proceed as follows. Following (Gibbs et al., ; Romano et al., ), we reformulate split conformal prediction as an augmented quantile regression. Let $S_{i}$ represent the non-conformity score of the sample $(X_{i},A_{i},Y_{i})$ for $i=m+1,\ldots,n$ of the calibration dataset and $S_{n+1}=S$ an imputed value for the unknown score of the new sample. We define

$$ $\displaystyle\hat{\theta}_{S}:=\operatorname*{arg\,min}\limits_{\theta\in \mathbb{R}}\frac{1}{n-m}\left(\sum_{i=m+1}^{n}l_{\alpha}(\theta,S_{i})+l_{ \alpha}(\theta,S)\right),$ (3) $$

where

$$ $\displaystyle l_{\alpha}(\theta,S):=\begin{cases}(1-\alpha)(S-\theta),&\text{ if }S\geq\theta,\\ \alpha(\theta-S),&\text{ if }S<\theta.\end{cases}$ (4) $$

Of note, $\hat{\theta}_{S}$ is an estimator of the $(1-\alpha)$-quantile of the non-conformity scores (Koenker & Bassett, ; Steinwart & Christmann, ). Using $\hat{\theta}_{S}$, we then construct the CP interval with the desired coverage guarantee $(1-\alpha)$. However, the interval is only valid for exchangeable data. Quantile regression might yield non-unique solutions, so we later restrict the analysis to solvers invariant to the data ordering.

### 4.1 Scenario 1: Known propensity score

We first consider scenario  1 with known propensity scores. Here, existing CP intervals are not directly applicable due to the shift from old to new propensity ($\rightarrow$ Challenge  a). For our derivation, we need the following lemma building upon and generalizing the intuition presented above.

###### Lemma 1 ( (Gibbs et al., 2023 ) ) .

Let $\mathcal{F}$ define a finite-dimensional function class that includes the function $f$ characterizing the shift in the (potentially unknown) propensity function $\pi$ (see Eq. equation 2). Define the distribution-shift-calibrated $(1-\alpha)$-quantile of the non-conformity scores as

$$ $\displaystyle\hat{g}_{S}(X_{n+1}):=\operatorname*{arg\,min}\limits_{g\in \mathcal{F}}\frac{1}{n-m}\left(\sum_{i=m+1}^{n}l_{\alpha}(g(X_{i}),S_{i})+l_{ \alpha}(g(X_{n+1}),S)\right)$ (5) $$

for an imputed guess $S$ of the $(n+1)$-th non-conformity score $S_{n+1}$. The prediction interval

$$ $\displaystyle C(X_{n+1}):=\{y\mid S_{n+1}(y)\leq\hat{g}_{S_{n+1}(y)}(X_{n+1})\}$ (6) $$

for the true $S_{n+1}$ given a realization of $Y_{n+1}=y$ satisfies the desired coverage guarantee under all distribution shifts $f\in\mathcal{F}$, i.e.,

$$ $\displaystyle P_{f}(Y_{n+1}\in C(X_{n+1}))\geq 1-\alpha.$ (7) $$

Building upon Lemma 1, we derive our first main result in Theorem 1. We define the finite-dimensional function class of interest as $\mathcal{F}:=\{\theta\frac{\pi(a+\Delta_{A}\mid x)}{\pi(a\mid x)}\mid\theta\in
\mathbb{R}^{+}\}$. It is easy to verify that all $f\in\mathcal{F}$ represent the desired propensity shift to $\tilde{\pi}(a\mid x)=\pi(a+\Delta_{A}\mid x)$ as defined in Eq. equation 2.

###### Theorem 1 (Conformal prediction intervals for known baseline policy) .

Consider a new datapoint with $X_{n+1}=x_{n+1}$, $A_{n+1}=a_{n+1}$, and $A^{\ast}(X_{n+1})=a^{*}=a_{n+1}+\Delta_{A}$. Let $\eta^{S}=\{\eta_{m+1}^{S},\ldots,\eta_{n+1}^{S}\}\in\mathbb{R}^{n+1}$ be the optimal solution to

$$ $\displaystyle\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)$ (8) s.t. $\displaystyle\qquad\qquad\qquad\qquad-\alpha\leq\eta_{i}\leq 1-\alpha,\quad \forall i=m+1,\ldots,n+1,$ $$

for an imputed unknown $S_{n+1}=S$.
Furthermore, let $S^{\ast}$ be defined as the maximum $S$ s.t. $\eta_{n+1}^{S}<1-\alpha$. Then, the prediction interval

$$ $\displaystyle C(x_{n+1},a^{*}):=\{y\mid S_{n+1}(y)\leq S^{\ast}\}$ (9) $$

satisfies the desired coverage guarantee

$$ $\displaystyle P(Y(A^{\ast}(X_{n+1}))\in C(X_{n+1},A^{*}(X_{n+1})))\geq 1-\alpha.$ (10) $$

###### Proof.

We provide a full proof in Supplement C.2.
Here, we briefly outline the underlying idea of the proof. First, we show that the function class $\mathcal{F}$ indeed satisfies Eq. equation 2 for the intervention $A^{*}(X)=A+\Delta_{A}$, and we then rewrite Eq. equation 5 as a convex optimization problem. Next, we exploit the strong duality property, we optimize over the corresponding dual problem to receive a dual prediction set with equal coverage probability. Finally, we derive $S^{\ast}$ from the dual prediction set to construct $C_{n+1}$ and prove the overall coverage guarantee.
∎

###### Lemma 1 ( (Gibbs et al., 2023 ) ) .

###### Theorem 1 (Conformal prediction intervals for known baseline policy) .

###### Proof.

### 4.2 Scenario 2: Unknown treatment policy

If the underlying treatment policy is unknown, the only possible intervention is a hard intervention $a^{\ast}$. As described above, measuring the induced propensity shift is non-trivial due to two reasons: (i) The propensity model needs to be estimated, which introduces additional uncertainty affecting the validity of the intervals ($\rightarrow$ Challenge  b). (ii) The density function corresponding to a hard intervention is given by the Dirac delta function

$$ $\displaystyle\delta_{a^{\ast}}(a):=\begin{cases}0,&\quad\text{for}\quad a\neq a ^{\ast},\\ \infty,&\quad\text{for}\quad a=a^{\ast},\end{cases}$ (11) $$

which hinders a direct adaptation of Theorem 1 due to the inherent discontinuity of the improper function. Hence, to proceed, we make the following assumption on the propensity estimator.

Assumption 1. *The estimation error of the propensity function $\hat{\pi}(a\mid x)$ is bounded in the sense that, for all $i=1,\ldots,n+1$, there exists $M>0$ such that*

$$ $\displaystyle c_{a_{i}}:=\frac{\hat{\pi}(a_{i}\mid x_{i})}{\pi(a_{i}\mid x_{i} )}\in\left[\tfrac{1}{M},M\right].$ (12) $$

Under Assumption 1, the distribution shift induced by the intervention is then defined as

$$ $\displaystyle\tilde{\pi}(a\mid x)=\frac{\delta_{a^{\ast}}(a)}{\hat{\pi}(a\mid x )}\frac{\hat{\pi}(a\mid x)}{\pi(a\mid x)}\pi(a\mid x)=c_{a}\frac{\delta_{a^{ \ast}}(a)}{\hat{\pi}(a\mid x)}\pi(a\mid x)=\frac{f(a,x)}{\mathbb{E}_{P}[f(A,X) ]}\pi({a\mid x}),$ (13) $$

for a suitable function $f$. We further formulate $\delta_{a^{\ast}}(a)$ in terms of a Gaussian function as

$$ $\displaystyle\delta_{a^{\ast}}(a)=\lim\limits_{\sigma\rightarrow 0}\frac{1}{ \sqrt{2\pi}\sigma}\exp\left(-\frac{(a-a^{\ast})^{2}}{2\sigma^{2}}\right).$ (14) $$

This motivates the following lemma. Therein, we specify the class $\mathcal{F}$ of tilting functions $f$ that represent the distribution shift induced by the hard intervention $a^{\ast}$.

###### Lemma 2 .

For $\sigma>0$, we define

$$ $\displaystyle f(a,x):=\frac{c_{a}}{\sqrt{2\pi}\sigma}\frac{\exp(-\frac{(a-a^{ \ast})^{2}}{2\sigma^{2}})}{\hat{\pi}(a\mid x)}$ (15) $$

with $\mathbb{E}_{P}[f(A,X)]=1$.
Furthermore, we define the finite-dimensional function class $\mathcal{F}$

$$ $\displaystyle\mathcal{F}:=\Bigg{\{}\frac{c_{a}}{\sqrt{2\pi}\sigma}\frac{\exp(- \frac{(a-a^{\ast})^{2}}{2\sigma^{2}})}{\hat{\pi}(a\mid x)}\Bigm{|}0<\sigma, \frac{1}{M}\leq c_{a}\leq M\Bigg{\}}.$ (16) $$

Then, $f(a,x)\in\mathcal{F}$ for all $c_{a}\in[\frac{1}{M},M]$ and $\sigma\rightarrow 0$. As a result, the distribution shift

$$ $\displaystyle\tilde{\pi}(a\mid x)=\lim\limits_{\sigma\rightarrow 0}\frac{c_{a} }{\sqrt{2\pi}\sigma}\,\frac{\exp(-\frac{(a-a^{\ast})^{2}}{\sigma^{2}})}{\hat{ \pi}(a\mid x)}\pi({a\mid x})$ (17) $$

can be represented in terms of Eq. equation 2 through functions $f\in\mathcal{F}$.

###### Proof.

See Supplement C.1.
∎

Following the motivation in scenario  1, we thus aim to estimate the $(1-\alpha)$-quantile of the non-conformity scores under the distribution shift in Lemma 1. We can reformulate this problem as

$$ $\displaystyle\min\limits_{{\sigma>0,\frac{1}{M}\leq c_{a}\leq M}}\quad\sum_{i= m+1}^{n+1}(1-\alpha)u_{i}+\alpha v_{i}$ (P S ) s.t. $\displaystyle\qquad\qquad S_{i}-\frac{c_{a}}{\sqrt{2\pi}\sigma}\frac{\exp{(- \frac{(a_{i}-a^{\ast})^{2}}{2\sigma^{2}})}}{\hat{\pi}(a_{i}\mid x_{i})}-u_{i}+ v_{i}=0,\quad\forall i=m+1\ldots,n+1,$ $\displaystyle\qquad\qquad u_{i},v_{i}\geq 0,\quad\forall i=m+1\ldots,n+1$ $$

for the imputed score $S_{n+1}=S$. As the score is unknown, computing the CP interval would require solving equation PS for all $S\in\mathbb{R}$, yet which is computationally infeasible. As a remedy, we previously exploited properties of the dual optimization problem and the Lagrange multipliers of the convex problem in Theorem 1 to efficiently compute the CP intervals (see the proof of Theorem 1 in Supplement C). However, the present non-convex problem does not automatically allow for the same simplifications. Instead, we now present a remedy for efficient computation of the CP intervals in the following lemma. We prove Lemma 3 in Supplement C.

###### Lemma 3 .

The problem equation PS is Type-I invex and satisfies the linear independence constraint qualification (LICQ).

Lemma 3 allows us to derive properties of the present non-convex optimization problem in terms of the Karush-Kuhn-Tucker (KKT) conditions. For this, we note that the fulfillment of the LICQ serves as a sufficient regularity condition for the KKT to hold at any (local) optimum of equation PS. Combined with the Type-I invexity of the objective function and the constraints, the KKT conditions are not only necessary but also sufficient for a global optimum. As a result, we can employ the KKT conditions at the optimal values(^3^33We provide an interpretation of the optimal values in Supplement E. Therein, we further discuss the implications of the proposed kernel smoothing of $\delta_{a^{\ast}}(a)$.) $\sigma^{\ast}$ and $c_{a}^{\ast}$ to derive coverage guarantees of our CP interval in a similar fashion as in Theorem 1. We thus arrive at the following theorem to provide CP intervals for the scenario with unknown propensity scores.

###### Theorem 2 (Conformal prediction intervals for unknown propensity scores) .

Let $u^{S}=\{u_{m+1}^{S},\ldots,u_{n+1}^{S}\},v^{S}=\{v_{m+1}^{S},\ldots,v_{n+1}^{S
}\}\in\mathbb{R}^{n-m}$, $\sigma^{S},c_{a}^{S}\in\mathbb{R}$ be the optimal solution to

$$ $\displaystyle\min\limits_{{\sigma>0,\frac{1}{M}\leq c_{a}\leq M}}\quad\sum_{i= m+1}^{n+1}(1-\alpha)u_{i}+\alpha v_{i}$ (18) s.t. $\displaystyle\qquad\qquad S_{i}-\frac{c_{a}}{\sqrt{2\pi}\sigma}\frac{\exp{(- \frac{(a_{i}-a^{\ast})^{2}}{2\sigma^{2}})}}{\hat{\pi}(a_{i}\mid x_{i})}-u_{i}+ v_{i}=0,\quad\forall i=m+1\ldots,n+1,$ $\displaystyle\qquad\qquad u_{i},v_{i}\geq 0,\quad\forall i=m+1\ldots,n+1.$ $$

for an imputed unknown $S_{n+1}=S$.
Let $S^{\ast}$ be defined as the maximum $S$ s.t. $v_{n+1}^{S}>0$. Then, the interval

$$ $\displaystyle C(X_{n+1},a^{\ast}):=\{y\mid S_{n+1}(y)\leq S^{\ast}\}$ (19) $$

satisfies the desired coverage guarantee

$$ $\displaystyle P(Y(a^{\ast})\in C(X_{n+1},a^{\ast}))\geq 1-\alpha.$ (20) $$

###### Proof.

See Supplement C.3.
∎

In certain applications, it might be beneficial to fix $\sigma$ to a small value $\sigma_{0}$ to approximate $\delta_{a^{\ast}}(a)$ though a soft intervention and only construct the CP interval through optimizing over $c_{a}$. We present an alternative theorem for this case in Supplement A.

###### Lemma 2 .

###### Proof.

###### Lemma 3 .

###### Theorem 2 (Conformal prediction intervals for unknown propensity scores) .

###### Proof.

### 4.3 Algorithm

Figure: Algorithm 1 Algorithm for computing CP intervals of potential outcomes of continuous interventions under unknown propensities.

We now use our Theorem 2 to present an algorithm for computing CP intervals of potential outcomes from continuous treatment variables under unknown propensities in Alg. 1. We present a similar algorithm for scenario  1 with known propensities and discuss the computational complexity in Supplement B.

We make the following comments: In our algorithm, an optimization solver is used to calculate $v_{n+1}$ according to Theorem 2. The specific choice of the solver is left to the user. In our experiments in Section 5, we perform the optimization via interior point methods. Further, the overall goal of our algorithm is to find the optimal imputed non-conformity score $S^{\ast}$ such that the coverage guarantees hold. This search can be implemented through suitable iterative search algorithms.

## 5 Experiments

Baselines: As we have discussed above, there are *no* baselines that directly compute finite-sample prediction intervals for potential outcomes of continuous treatments. Therefore, we compare our method against MC dropout (Gal & Ghahramani, ) and deep ensemble methods (Lakshminarayanan et al., ). Yet, we again emphasize that MC dropout is an ad hoc method with poor approximations of the posterior, which is known to give intervals that are *not* faithful (Le Folgoc et al., ).

Implementation: All methods are implemented with $\phi$ as a multi-layer perception (MLP) and an MC dropout regularization of rate 0.1. Crucially, we use the *identical* MLP for both our CP method and MC dropout. Hence, all performance gains must be attributed to the coverage guarantees of our conformal method. In the MC dropout baseline, the uncertainty intervals are computed via Monte Carlo sampling. In scenario 2 with unknown propensity scores, we perform the conditional density estimation through conditional normalizing flows (Trippe & Turner, ). Details about our implementation and training are in Supplement E.

Performance metrics: We evaluate the methods in terms of whether the prediction intervals are *faithful* (e.g., Hess et al., ). That is, we compute whether the *empirical coverage* of the prediction intervals surpasses the threshold of $1-\alpha$ for different significance levels $\alpha\in\{0.05,0.1,0.2\}$.

### 5.1 Datasets

Synthetic datasets: We follow common practice in causal ML and evaluate our methods using synthetic datasets (e.g., Alaa et al., ; Jin et al., ). The reason is the fundamental problem of causal inference, because of which counterfactual outcomes are never observable in real-world datasets. Synthetic datasets enable us to access counterfactual outcomes and, thus, to benchmark methods in terms of whether the computed intervals are faithful.

We consider two synthetic datasets with different propensity scores and outcome functions. Dataset 1 uses a step-wise propensity function and a concave outcome function. Dataset 2 is more complex and uses a Gaussian propensity function and oscillating outcome functions. Both datasets contain a single discrete confounder, a continuous treatment, and a continuous outcome. By choosing low-dimensional datasets, we later render it possible to plot the treatment–response curves so that one can inspect the prediction intervals visually. (We later also show that our method scales to high-dimensional settings as part of the real-world dataset.) Details about the data-generating processes are in Supplement E.

Figure: Figure 4: Comparison of *faithfulness* on dataset 1 across 50 runs. Larger values are better. For each $\alpha$, the plots show how often the empirical intervals contain the true outcome. Intervals should ideally yield a coverage of $1-\alpha$ (red line).
Refer to caption: x4.png

Medical dataset: We demonstrate the applicability of our CP method to datasets from medicine by leveraging the MIMIC-III dataset (Johnson et al., ). MIMIC-III contains de-identified health records from patients admitted to critical care units at large tertiary care hospitals. Our goal is to predict patient outcomes in terms of blood pressure when treated with a different duration of mechanical ventilation. We use 8 confounders from medical practice (e.g., respiratory rate, hematocrit). Overall, we consider 14,719 patients, which we randomly split into train (60%), validation (10%), calibration (20%), and test (10%). Further details are in Supplement E.

### 5.2 Results for synthetic datasets

Figure: Figure 5: Comparison of *faithfulness* on dataset 2 across 50 runs. Larger values are better.
Refer to caption: x5.png

We now evaluate the faithfulness of our CP intervals through a series of experiments. On each dataset, we analyze the performance of the prediction intervals in the presence of various soft interventions $\Delta_{a}\in\{1,5,10\}$ and hard interventions $a^{\ast}\in\{5x,7x,10x\}$ for each $X=x$. We average the empirical coverage across 10 runs with random seeds. The results are in Fig. 4 (dataset 1) and Fig. 5 (dataset 2).

We make the following observations. First, the intervals of our CP method comply with the targeted significance level $\alpha$ and, therefore, are faithful. Second, both MC dropout and the deep ensemble method have, in contrast, a considerably lower coverage, implying that the intervals are *not* faithful. This is in line with the literature, where MC dropout is found to produce poor approximations of the posterior (Le Folgoc et al., ). In particular, the ensemble method is highly unfaithful. Thus, we will not consider this baseline in all of the following experiments. Third, our method has only a small variability in terms of empirical coverage, whereas the empirical coverage of MC dropout varies greatly. This corroborates the robustness of our method. Fourth, the results are consistent for both datasets.
In sum, this demonstrates the effectiveness of our proposed CP method.

**Table 1: Coverage of the intervals from our CP method / MC dropout for various hard interventions $a^{\ast}$ and significance levels $\alpha$. Intervals with coverage $\geq 1-\alpha$ are considered faithful (bold numbers).**
|  | Coverage on dataset 1 | Coverage on dataset 2 |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Intervention | $\alpha=0.05$ | $\alpha=0.1$ | $\alpha=0.2$ | $\alpha=0.05$ | $\alpha=0.1$ | $\alpha=0.2$ |
| $a^{\ast}=7x$ | 1.00 / 0.19 | 0.90 / 0.13 | 0.83 / 0.11 | 1.00 / 0.02 | 0.94 / 0.02 | 0.85 / 0.02 |
| $a^{\ast}=10x$ | 1.00 / 0.28 | 0.91 / 0.23 | 0.88 / 0.11 | 1.00 / 0.08 | 0.84 / 0.07 | 0.83 / 0.07 |

Table LABEL:tab:results_unknown presents the empirical coverage of the intervals from our CP method vs. MC dropout across different $\alpha$ and hard interventions $a^{\ast}$. We observe that our conformal prediction intervals are effective and achieve the intended coverage. In contrast, MC dropout does not provide faithful intervals. Our findings are again in line with the literature, where MC dropout is found to produce poor approximations of the posterior and thus might provide poor coverage (Le Folgoc et al., ). We present further results in Supplement F.

Figure: Figure 6: Prediction intervals for multiple significance levels $\alpha$ for synthetic dataset 1 with intervention $\Delta=5$.
Refer to caption: x6.png

Insights: We plot the prediction intervals across different significance levels $\alpha$ and different covariates $X$ (see Fig. 6). This allows us to inspect the intervals visually.
We observe that the intervals behave as expected; that is, the intervals become sharper with increasing significance level $\alpha$. We further see that our CP intervals are slightly wider (see details in Supplement F), yet this is intended because it ensures that the intervals are faithful. For example, our CP intervals (blue) generally include the true outcome. In contrast, the intervals from MC dropout (orange) often do *not* include the true outcome (e.g., see the bottom row Fig. 6) and are thus *not* faithful.

### 5.3 Results for the MIMIC dataset

In Fig. 7, we compare the CP intervals of each two male and female patients of differing ages when treated with increasing duration of mechanical ventilation. Our intervals show higher uncertainty in treatment regions rarely included in the training data (high medium to high treatments). The intervals given by MC-dropout are narrower, which also indicates lower cover-
age, thus confirming the effectiveness of our proposed method. This finding aligns with our observation from the synthetic datasets.

Figure: Figure 7: CP intervals for potential outcomes of increasing duration of mechanical ventilation for four exemplary patients.
Refer to caption: x7.png

## 6 Discussion

Limitations:
As with any other method, our UQ method has limitations that offer opportunities for future research. Our method relies on the quality of the propensity estimator. Although we incorporate estimation errors in the construction of our intervals, poorly estimated propensities could potentially lead to wide prediction intervals. We acknowledge that our prediction intervals are conservative for the point intervention and for segments of the output space with limited calibration data, implying that a representative calibration dataset is essential for the performance of our method. As for all CP methods, the use of sample splitting may reduce data efficiency.

Broader impact: Our method makes a significant step toward UQ for potential outcomes and, therefore, toward *reliable* decision-making. To this end, we provided strong theoretical and empirical evidence that our prediction intervals are valid. To this end, our method fills an important demand for using causal ML in medical practice Feuerriegel et al. () and other safety-critical applications with limited data.

Conclusion: We presented a novel conformal prediction method for potential outcomes of continuous treatments with finite-sample guarantees. Our method is applicable to potential outcomes and extends naturally to treatment effects. A key strength of our method is that the intervals are valid under distribution shifts introduced by the treatment assignment, even if the propensity score is unknown and has to be estimated.

## Acknowledgements

This paper is supported by the DAAD program Konrad Zuse Schools of Excellence in Artificial Intelligence, sponsored by the Federal Ministry of Education and Research.

## Appendix A Additional theoretical results

### A.1 Calculating prediction intervals for further causal quantities and differences

We presented a method for calculating conformal prediction intervals for potential outcomes of continuous treatments. In the following, we show how the intervals can be combined to yield valid prediction intervals for further causal quantities, such as the individual treatment effect (ITE) $\gamma_{i}$ of treatment $a$:

$$ $\displaystyle\gamma_{i}(a):=Y_{i}(a)-Y_{i}(0).$ (21) $$

###### Lemma 4 .

Let $S_{a}^{\ast}$ and $S_{0}^{\ast}$ denote the optimal imputed non-conformity scores $S_{n+1}$ for treatment $a$ and no treatment at significance level $1-\frac{\alpha}{2}$ for $\alpha\in(0,1)$, respectively. Furthermore, let

$$ $\displaystyle C^{+}:=\phi(x_{i},a)+S_{a}^{\ast}-\phi(x_{i},0)+S_{0}^{\ast},$ (22) $\displaystyle C^{-}:=\phi(x_{i},a)-S_{a}^{\ast}-\phi(x_{i},0)-S_{0}^{\ast}.$ (23) $$

Then the interval $C_{\gamma}(X_{i},a):=[C^{-},C^{+}]$ contains the ITE $\gamma_{i}$ with probability $1-\alpha$.

###### Proof.

Let $\varepsilon_{i}(a)$ be the estimation error of the potential outcome, i.e.

$$ $\displaystyle\varepsilon_{i}(a):=Y_{i}(a)-\phi(x_{i},a).$ (24) $$

We can rewrite the coverage guarantee of the original conformal prediction intervals for the potential outcome $Y(a)$ as

$$ $\displaystyle P(Y_{i}(a)\in C(X_{i},a))=P(\lvert\varepsilon_{i}(a)\rvert\leq S _{a}^{\ast})\geq 1-\frac{\alpha}{2}.$ (25) $$

Now observe that

$$ $\displaystyle P(\gamma_{i}(a)\in C_{\gamma}(x_{i},a))=P((Y_{i}(a)-Y_{i}(0))\in C _{\gamma}(x_{i},a))$ (26) $\displaystyle=$ $\displaystyle P((Y_{i}(a)\geq C^{-}+Y_{i}(0))\land(Y_{i}(a)\leq C^{+}+Y_{i}(0)))$ (27) $\displaystyle=$ $\displaystyle P((\varepsilon_{i}(a)\geq\varepsilon_{i}(0)-(S_{a}^{\ast}+S_{0}^ {\ast}))\land(\varepsilon_{i}(a)\leq\varepsilon_{i}(0)+(S_{a}^{\ast}+S_{0}^{ \ast})))$ (28) $\displaystyle=$ $\displaystyle P(\lvert\varepsilon_{i}(a)-\varepsilon_{i}(0)\rvert\leq S_{a}^{ \ast}+S_{0}^{\ast})$ (29) $\displaystyle\geq$ $\displaystyle P(\lvert\varepsilon_{i}(a)\rvert+\lvert\varepsilon_{i}(0)\rvert \leq S_{a}^{\ast}+S_{0}^{\ast}).$ (30) $$

Thus, it follows directly that

$$ $\displaystyle P(\gamma_{i}(a)\in C_{\gamma}(X_{i},a))\geq 1-\alpha.$ (31) $$

∎

###### Lemma 4 .

###### Proof.

### A.2 Alternative scenario 2: Fixing an approximation of δ a ∗ ⁢ ( a ) subscript 𝛿 superscript 𝑎 ∗ 𝑎 \delta_{a^{\ast}}(a) italic_δ start_POSTSUBSCRIPT italic_a start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT ( italic_a )

In Section 4.2, we formulated the unknown propensity shift in terms of

$$ $\displaystyle\delta_{a^{\ast}}(a)=\lim\limits_{\sigma\rightarrow 0}\frac{1}{ \sqrt{2\pi}\sigma}\exp\left(-\frac{(a-a^{\ast})^{2}}{\sigma^{2}}\right)$ (32) $$

and minimized over $\sigma$ and $c_{a}$ in Theorem 2 to construct the CP intervals. However, in certain applications, it might be beneficial to control the spread of the approximation of $\delta_{a^{\ast}}(a)$ through fixing $\sigma$ to a small value $\sigma_{0}$ and performing the soft intervention $\tilde{\pi}(a\mid x)=\frac{c_{a}}{\sqrt{2\pi}\sigma_{0}}\frac{\exp{(-\frac{(a-
a^{\ast})^{2}}{\sigma_{0}^{2}})}}{\hat{\pi}(a\mid x)}$. In this case, the resulting optimization problem is a convex problem similar to Theorem 1. We present the alternative optimization problem below.

###### Theorem 3 (Alternative for Theorem 2 : Conformal prediction intervals for unknown propensity scores) .

Let a new datapoint be given with $X_{n+1}=x_{n+1}$ and $A_{n+1}=a_{n+1}$. Let $\eta^{S}=\{\eta_{1}^{S},\ldots,\eta_{n+1}^{S}\}\in\mathbb{R}^{n+1}$ be the optimal solution to

$$ $\displaystyle\max_{\begin{subarray}{c}\eta_{i},\\ i=1,\ldots,n+1\end{subarray}}\min_{\frac{1}{M}\leq c_{a}\leq M}\sum_{i=1}^{n} \eta_{i}\,\left(S_{i}-\frac{c_{a}}{\sqrt{2\pi}\sigma_{0}}\frac{\exp{\left(- \frac{(a_{i}-a^{\ast})^{2}}{\sigma_{0}^{2}}\right)}}{\hat{\pi}(a_{i}\mid x_{i} )}\right)+\eta_{n+1}\left(S-\frac{c_{a}}{\sqrt{2\pi}\sigma_{0}}\frac{1}{\hat{ \pi}(a_{i}\mid x_{n+1})}\right)$ (33) $\displaystyle\text{ s.t.}\qquad\qquad\qquad\qquad-\alpha\leq\eta_{i}\leq 1- \alpha,\quad\forall i=1,\ldots,n+1,$ $$

for an imputed unknown $S_{n+1}=S$.
Furthermore, let $S^{\ast}$ be defined as the maximum $S$ s.t. $\eta_{n+1}^{S}<1-\alpha$. Then, the prediction interval

$$ $\displaystyle C(x_{n+1},a^{*}):=\{y\mid S_{n+1}(y)\leq S^{\ast}\}$ (34) $$

satisfies the desired coverage guarantee

$$ $\displaystyle P(Y(a^{\ast}\in C(X_{n+1},a^{*})\geq 1-\alpha,$ (35) $$

where with a slight abuse of notation $Y(a^{\ast})$ denotes the potential outcome under the soft intervention $\tilde{\pi}$ above.

###### Proof.

The statement follows from Theorem 1.
∎

###### Theorem 3 (Alternative for Theorem 2 : Conformal prediction intervals for unknown propensity scores) .

###### Proof.

## Appendix B Algorithm

In our main paper, we presented an algorithm for computing the prediction intervals if the propensity score is unknown. Below, we state a second algorithm that is applicable if the propensity score is known. In this case, a convex optimization solver can be used.

Figure: Algorithm 2 Algorithm for computing CP intervals of potential outcomes of continuous interventions under known propensities.

Computational complexity:
The complexity of running our algorithms depends heavily on the employed optimization solver with complexity $\sigma_{s}(n_{c})$ (e.g., polynomial complexity for suitable convex solvers) and the size of the calibration dataset $n_{c}$. This might become costly for large-scale calibration datasets in practice. The outer algorithm has a time complexity of at most $O(\log(\frac{S_{up}-S_{low}}{\varepsilon})+1)$. Overall, our algorithm has a fixed complexity of $O(\log(\frac{S_{up}-S_{low}}{\varepsilon})+\sigma_{s}(n_{c}))$. The complexity of deriving intervals through MC-dropout or ensemble methods depends, however, on the number of MC samples or models, respectively. The latter thus scales with the precision of the intervals, which might be difficult to control. 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.

## Appendix C Proofs

### C.1 Proofs of the supporting lemmas

In the following, we prove Lemma 2 and Lemma 3 from our main paper.

#### Proof of Lemma 2

Recall the definition of the hard intervention

$$ $\displaystyle\tilde{\pi}(a\mid x)=\delta_{a^{\ast}}(a)=\frac{\delta_{a^{\ast}} (a)}{\hat{\pi}(a\mid x)}\frac{\hat{\pi}(a\mid x)}{\pi(a\mid x)}\pi(a\mid x),$ (36) $$

where

$$ $\displaystyle\delta_{a^{\ast}}(a)=\lim\limits_{\sigma\rightarrow 0}\frac{1}{ \sqrt{2\pi}\sigma}\exp\left(-\frac{(a-a^{\ast})^{2}}{\sigma^{2}}\right).$ (37) $$

Under Assumption 1, we have

$$ $\displaystyle\frac{\hat{\pi}(a\mid x)}{\pi(a\mid x)}=:c_{a}\in[\frac{1}{M},M]$ (38) $$

for some $M>0$ and all $a,x$. Then

$$ $\displaystyle\lim\limits_{\sigma\rightarrow 0}\frac{1}{\sqrt{2\pi}\sigma}\frac {\exp\left(-\frac{(a-a^{\ast})^{2}}{\sigma^{2}}\right)}{\hat{\pi}(a\mid x)} \frac{1}{M}\leq\pi(a^{\ast}|x)\leq\lim\limits_{\sigma\rightarrow 0}\frac{1}{ \sqrt{2\pi}\sigma}\frac{\exp\left(-\frac{(a-a^{\ast})^{2}}{\sigma^{2}}\right)} {\hat{\pi}(a\mid x)}M.$ (39) $$

Therefore, the distribution shift induced by the hard intervention can be represented as

$$ $\displaystyle f(a,x)=\lim\limits_{\sigma\rightarrow 0}\frac{c_{a}}{\sqrt{2\pi} \sigma}\frac{\exp\left(-\frac{(a-a^{\ast})^{2}}{\sigma^{2}}\right)}{\hat{\pi}( a\mid x)}\in\mathcal{F}:=\left\{\frac{c_{a}}{\sqrt{2\pi}\sigma}\frac{\exp\left (-\frac{(a-a^{\ast})^{2}}{\sigma^{2}}\right)}{\hat{\pi}(a\mid x)}\;\middle|\;0 <\sigma,c_{a}\in[\frac{1}{M},M]\right\}.$ (40) $$

∎

#### Proof of Lemma 3

We first prove that the problem equation PS fulfills the linear independence constraint qualifications. For all $i=m+1,\ldots,n+1$, we denote the constraints of problem equation PS as

$$ $\displaystyle h_{i}(u,v,c_{a},\sigma):=S_{i}-u_{i}+v_{i}-\frac{c_{a}}{\sqrt{2 \pi}\sigma}\frac{\exp(-\frac{(a_{i}-a^{\ast})^{2}}{2\sigma})}{\hat{\pi}(a_{i}, x_{i})}.$ (41) $$

The gradient of $h_{i}$ is given by

$$ $\displaystyle\nabla h_{i}(u,v,c_{a},\sigma)=\left[\begin{array}[]{c}\dfrac{ \partial h_{i}}{\partial u_{m+1}}(\left.u,v,c_{a},\sigma\right)\\ \dfrac{\partial h_{i}}{\partial v_{m+1}}(\left.u,v,c_{a},\sigma\right)\\ \vdots\\ \dfrac{\partial h_{i}}{\partial u_{i}}(\left.u,v,c_{a},\sigma\right)\\ \dfrac{\partial h_{i}}{\partial v_{i}}(\left.u,v,c_{a},\sigma\right)\\ \vdots\\ \dfrac{\partial h_{i}}{\partial c_{a}}(\left.u,v,c_{a},\sigma\right)\\ \dfrac{\partial h_{i}}{\partial\sigma}(\left.u,v,c_{a},\sigma\right)\end{array }\right]=\left[\begin{array}[]{c}0\\ 0\\ \vdots\\ -1\\ 1\\ \vdots\\ -\frac{1}{\sqrt{2\pi}\sigma}\frac{\exp(-\frac{(a_{i}-a^{\ast})^{2}}{2\sigma})} {\hat{\pi}(a_{i},x_{i})}\\ \frac{\exp(-\frac{(a_{i}-a^{\ast})^{2}}{2\sigma})}{\hat{\pi}(a_{i},x_{i})}( \frac{c_{a}}{\sqrt{2\pi}\sigma^{2}}-\frac{c_{a}(a_{i}-a^{\ast})^{2}}{2\sqrt{2 \pi}\sigma^{4}})\end{array}\right].$ (58) $$

Therefore, with $\nabla h:=(\nabla h_{m+1},\ldots,\nabla h_{n+1})$ and $\lambda\in\mathbb{R}^{n+1}$, we obtain

$$ $\displaystyle\nabla h\cdot\lambda=0\iff\lambda=0\in\mathbb{R}^{n+1}.$ (59) $$

As a result, the constraints are linearly independent. This property suffices for the KKT conditions to hold at any (local) optimum of equation PS.To furthermore show that the KKT conditions are also sufficient for a global optimum, we show that equation PS is Type-I invex. An optimization problem with objective function $f(x)$ and constraints $g(x)<=0$ with $x\in\mathbb{R}^{n+1}$ is Type-I invex at $x_{0}$, if there exists $\nu(x,x_{0})\in\mathbb{R}^{n+1}$, such that

$$ $\displaystyle f(x)-f(x_{0})\geq\nu(x,x_{0})^{T}\nabla f(x_{0}),$ (60) $$

and

$$ $\displaystyle-g(x_{0})\geq\nu(x,x_{0})^{T}\nabla g(x_{0})$ (61) $$

(Hanson & Mond, ).
In problem equation PS, the gradients of the objective function and of each constraint $h_{i}$ for all $i,j=m+1,\ldots,n+1$ at $x_{0}$ are given by

$$ $\displaystyle\frac{\partial\text{obj}(u_{0},v_{0},c_{a_{0}},\sigma_{0})}{ \partial u_{i}}=1-\alpha,\quad\frac{\partial\text{obj}(u_{0},v_{0},c_{a_{0}}, \sigma_{0})}{\partial v_{i}}=\alpha,$ (62) $\displaystyle\frac{\partial\text{obj}(u_{0},v_{0},c_{a_{0}},\sigma_{0})}{ \partial c_{a}}=\frac{\partial\text{obj}(u_{0},v_{0},c_{a_{0}},\sigma_{0})}{ \partial\sigma}=0$ (63) $$

$\forall i=m+1,\ldots,n+1$ and

$$ $\displaystyle\frac{\partial h_{i}}{\partial u_{j}}\rvert_{u_{0},v_{0},c_{a_{0} },\sigma_{0}}=\begin{cases}-1,\qquad&\text{for }i=j,\\ 0,&\text{else},\end{cases}\qquad\frac{\partial h_{i}}{\partial v_{j}}\rvert_{u _{0},v_{0},c_{a_{0}},\sigma_{0}}=\begin{cases}1,\qquad&\text{for }i=j,\\ 0,&\text{else},\end{cases}$ (64) $\displaystyle\frac{\partial h_{i}}{\partial c_{a}}\rvert_{u_{0},v_{0},c_{a_{0} },\sigma_{0}}=-\frac{1}{\sqrt{2\pi}\sigma_{0}}\frac{\exp\left(-\frac{(a_{i}-a^ {\ast})^{2}}{2\sigma_{0}}\right)}{\hat{\pi}(a_{i},x_{i})},$ (65) $\displaystyle\frac{\partial h_{i}}{\partial\sigma}\rvert_{u_{0},v_{0},c_{a_{0} },\sigma_{0}}=\frac{\exp\left(-\frac{(a_{i}-a^{\ast})^{2}}{2\sigma_{0}}\right) }{\hat{\pi}(a_{i},x_{i})}\left(\frac{c_{a_{0}}}{\sqrt{2\pi}\sigma_{0}^{2}}- \frac{c_{a_{0}}(a_{i}-a^{\ast})^{2}}{2\sqrt{2\pi}\sigma_{0}^{4}}\right).$ (66) $$

For

$$ $\displaystyle\eta((u,v,c_{a},\sigma),(u_{0},v_{0},c_{a_{0}},\sigma_{0})):=(-u_ {0_{1}},\ldots,-u_{0_{n+1}},-v_{0_{1}},\ldots,-v_{0_{n+1}},-c_{a_{0}},0)^{T},$ (67) $$

the definition of Type-I invexity holds for equation PS. Thus, the KKT conditions are also sufficient for a global optimum.
∎

### C.2 Proof of Theorem 1

We prove Theorem 1 in three steps: (i) We show that function class $\mathcal{F}:=\{\theta\frac{\pi(a+\Delta_{A}\mid x)}{\pi(a\mid x)}\mid\theta\in
\mathbb{R}^{+}\}$ indeed satisfies Eq. equation 2 for the intervention $a^{\ast}=a+\Delta_{A}$ and rewrite Eq. equation 5 as a convex optimization problem. (ii) We retrieve the corresponding dual problem, derive a dual prediction set, and show the equality of the coverage guarantee of the dual and the primal prediction sets. (iii) We derive $S^{\ast}$ from the dual prediction set to construct $C_{n+1}$ and prove the overall coverage guarantee. For further theoretical background on the idea of the proof, we refer to Gibbs et al. ().

#### Justification of the distribution shift

Observe that $\mathbb{E}[f(A,X)]=\theta$ for all $f\in\mathcal{F}:=\{\theta\frac{\pi(a+\Delta_{A}\mid x)}{\pi(a\mid x)}\mid
\theta\in\mathbb{R}^{+}\}$. Therefore, Eq. equation 2 simplifies to

$$ $\displaystyle\tilde{\pi}(a,x)=\frac{\pi(a+\Delta_{A}\mid x)}{\pi(a\mid x)}\pi( a\mid x)=\pi(a+\Delta_{A}\mid x).$ (68) $$

Thus, $\mathcal{F}$ satisfies the propensity shift from Eq. equation 2 for the soft intervention $a^{\ast}=a+\Delta_{A}$. Following Lemma 1, we thus aim to find

$$ $\displaystyle\hat{q}_{S}:=\text{arg\,}\min\limits_{\theta>0}\frac{1}{n-m}\left (\sum_{i=m+1}^{n}l_{\alpha}(\theta\frac{\pi(a_{i}+\Delta_{A}\mid x_{i})}{\pi(a _{i}\mid x_{i})},S_{i})+l_{\alpha}(\theta\frac{\pi(a^{\ast}\mid x_{n+1})}{\pi( a_{n+1}\mid x_{n+1})},S)\right).$ (69) $$

#### Dual problem formulation

First, we rewrite the primal problem as

$$ $\displaystyle\min_{\theta>0}\quad$ $\displaystyle\sum_{i=m+1}^{n+1}(1-\alpha)u_{i}+\alpha v_{i}$ (70) s.t. $\displaystyle 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$ $\displaystyle u_{i},v_{i}\geq 0,\quad\forall i=m+1\ldots,n+1.$ $$

For a reference, see (Gibbs et al., ).
The Lagrangian of the primal problem states

$$ $\displaystyle\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}).$ (71) $$

Setting derivative of $\mathcal{L}$ w.r.t. $u_{i}$ and $v_{i}$ to 0 0 results in

$$ $\displaystyle\frac{\partial\mathcal{L}}{\partial u_{i}}=(1-\alpha)-\eta_{i}- \gamma_{1_{i}}\overset{!}{=}0,\quad\forall i=m+1\ldots,n+1$ (72) $\displaystyle\frac{\partial\mathcal{L}}{\partial v_{i}}=(1-\alpha)-\eta_{i}- \gamma_{2_{i}}\overset{!}{=}0,\quad\forall i=m+1\ldots,n+1.$ (73) $$

Since $\gamma_{1_{i}},\gamma_{2_{i}}\geq 0$ $\forall i$, it follows for all $i=m+1,\ldots,n+1$ that

$$ $\displaystyle(1-\alpha)-\eta_{i}\geq 0\quad\text{and}\quad\alpha-\eta_{i}\geq 0 \quad\Rightarrow-\alpha\leq\eta_{i}\leq 1-\alpha.$ (74) $$

Therefore, the dual problem is formulated as

$$ $\displaystyle\max_{\eta_{i},i=m+1,\ldots,n+1}\min_{\theta>0}\quad$ $\displaystyle\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)$ (75) s.t. $\displaystyle-\alpha\leq\eta_{i}\leq 1-\alpha,\quad\forall i=m+1,\ldots,n+1.$ $$

#### Coverage guarantee

Recall from Lemma 1 that, for

$$ $\displaystyle\hat{q}_{S_{n+1}}(y)=\text{arg\,}\min\limits_{\theta>0}\frac{1}{n -m}(\sum_{i=m+1}^{n}l_{\alpha}(\theta\frac{\pi(a_{i}+\Delta_{A}\mid x_{i})}{ \pi(a_{i}\mid x_{i})},S_{i})+l_{\alpha}(\theta\frac{\pi(a^{\ast}\mid x_{n+1})} {\pi(a_{n+1}\mid x_{n+1})},S_{n+1}(y))),$ (76) $$

we can construct $C_{n+1}=\{y\mid S_{n+1}(y)\leq\hat{q}_{S_{n+1}}(y)\}$ to achieve the desired coverage guarantee

$$ $\displaystyle P_{f}(Y(a^{\ast})\in C(X_{n+1},A^{\ast}(X_{n+1})))\geq 1-\alpha.$ (77) $$

It is infeasible to calculate $\hat{q}_{S_{n+1}}(y)$ directly. Therefore, we optimize the dual problem to receive

$$ $\displaystyle C(X_{n+1},a^{\ast}):=\{y\mid S_{n+1}(y)\leq S^{\ast}\}$ (78) $$

with $S^{\ast}$ the maximum $S$, s.t. for $\eta_{n+1}^{S}$ maximizing the dual problem, $\eta_{n+1}^{S}<1-\alpha$. Hence, it is left to show that replacing $\hat{q}_{S_{n+1}}(y)$ by $S^{\ast}$ in $C$ does not change to coverage guarantee.

To do so, we fix some $\theta>0$ to obtain a specific $f(a,x):=\theta\frac{\pi(a^{\ast}\mid x)}{\pi(a\mid x)}$. Let $\hat{g}(a,x)\in\mathcal{F}$ denote the primal optimal solution. Recall the Lagrangian

$$ $\displaystyle\mathcal{L}=\sum_{i=m+1}^{n+1}(1-\alpha)u_{i}+\alpha v_{i}+\sum_{ i=m+1}^{n+1}\eta_{i}(S_{i}-f(a_{i},x_{i})-u_{i}+v_{i})-\sum_{i=m+1}^{n+1}( \gamma_{1_{i}}u_{i}+\gamma_{2_{i}}v_{i}).$ (79) $$

Deriving wrt. $f$ yields the stationarity condition

$$ $\displaystyle 0$ $\displaystyle\overset{!}{=}-\sum_{i=m+1}^{n+1}\eta_{i}^{S}f(a_{i},x_{i})$ (80) $\displaystyle=-\sum_{S_{i}<\hat{g}(a_{i},x_{i})}\eta_{i}^{S}f(a_{i},x_{i})- \sum_{S_{i}>\hat{g}(a_{i},x_{i})}\eta_{i}^{S}f(a_{i},x_{i})-\sum_{S_{i}=\hat{g }(a_{i},x_{i})}\eta_{i}^{S}f(a_{i},x_{i}).$ (81) $$

The complementary slackness Karush-Kuhn-Tucker conditions yield

$$ $\displaystyle\eta_{i}^{S}\in\begin{cases}-\alpha,&\text{if }S_{i}<\hat{g}(a_{i },x_{i}),\\ [-\alpha,1-\alpha],&\text{if }S_{i}=\hat{g}(a_{i},x_{i}),\\ 1-\alpha,&\text{if }S_{i}>\hat{g}(a_{i},x_{i}).\end{cases}$ (82) $$

Therefore, we can rewrite the equation from above as

$$ $\displaystyle 0$ $\displaystyle=\sum_{S_{i}<\hat{g}(a_{i},x_{i})}\alpha f(a_{i},x_{i})-\sum_{S_{ i}>\hat{g}(a_{i},x_{i})}(1-\alpha)f(a_{i},x_{i})-\sum_{S_{i}=\hat{g}(a_{i},x_{ i})}\eta_{i}^{S}f(a_{i},x_{i})$ (83) $\displaystyle=\sum_{\eta_{i}^{S}<1-\alpha}\alpha f(a_{i},x_{i})-\sum_{\eta_{i} ^{S}=1-\alpha}(1-\alpha)f(a_{i},x_{i})-\sum_{\begin{subarray}{c}\eta_{i}^{S}<1 -\alpha,\\ S_{i}=\hat{g}(a_{i},x_{i})\end{subarray}}(\alpha+\eta_{i}^{S})f(a_{i},x_{i})$ (84) $\displaystyle=\sum_{i=m+1}^{n+1}(\alpha-\mathbbm{1}_{[\eta_{i}^{S}=1-\alpha]}) f(a_{i},x_{i})-\sum_{\begin{subarray}{c}\eta_{i}^{S}<1-\alpha,\\ S_{i}=\hat{g}(a_{i},x_{i})\end{subarray}}(\alpha+\eta_{i}^{S})f(a_{i},x_{i}).$ (85) $$

Before deriving the coverage guarantee from the stationarity condition, we state the following lemma to underline the definition of $S^{\ast}$.

###### Lemma 5 ( Gibbs et al. ( 2023 ) ) .

The mapping $S\mapsto\eta_{n+1}^{S}$ is non-decreasing in S for all $\eta_{n+1}^{S}$ maximizing

$$ $\displaystyle\max_{\eta_{i},i=m+1,\ldots,n+1}\min_{g\in\mathcal{F}}\quad$ $\displaystyle\sum_{i=1}^{n}\eta_{i}(S_{i}-g(a_{i},x_{i}))+\eta_{n+1}(S-g(a_{n+ 1},x_{n+1}))$ (86) s.t. $\displaystyle-\alpha\leq\eta_{i}\leq 1-\alpha,\quad\forall i=m+1,\ldots,n+1$ $$

for non-negative function classes $\mathcal{F}$.

To prove the final coverage guarantee, we observe that

$$ $\displaystyle\mathbb{E}[f(a_{n+1},x_{n+1})(\mathbbm{1}_{[Y(a^{\ast})\in C(X_{n +1},a^{\ast})]}-(1-\alpha))]$ (87) $\displaystyle=$ $\displaystyle\mathbb{E}[f(a_{n+1},x_{n+1})(\alpha-\mathbbm{1}_{[Y(a^{\ast}) \notin C(X_{n+1},a^{\ast})]})]$ (88) $\displaystyle=$ $\displaystyle\mathbb{E}[f(a_{n+1},x_{n+1})(\alpha-\mathbbm{1}_{[S(y)>S^{\ast}] })].$ (89) $$

With the definition of $S^{\ast}$ as the maximum optimizer $\eta_{n+1}^{S}$ with $\eta_{n+1}^{S}<1-\alpha$ and Lemma 5, it follows that

$$ $\displaystyle\mathbb{E}[f(a_{n+1},x_{n+1})(\alpha-\mathbbm{1}_{[S(y)>S^{\ast}] })]=\mathbb{E}[(\alpha-\mathbbm{1}_{[\eta_{n+1}^{S}=1-\alpha]})f(a_{n+1},x_{n+ 1})]$ (90) $$

and, by exchangeability of $(f(a_{i},x_{i}),\hat{q}_{S}(a_{i}.x_{i}),S_{i})$, that

$$ $\displaystyle\mathbb{E}[(\alpha-\mathbbm{1}_{[\eta_{i}^{S}=1-\alpha]})f(a_{i}, x_{i})]$ $\displaystyle=\mathbb{E}\left[\frac{1}{n-m}\sum_{i=m+1}^{n+1}(\alpha-\mathbbm{ 1}_{[\eta_{i}^{S}=1-\alpha]})f(a_{i},x_{i})\right]$ (91) $\displaystyle=\mathbb{E}\left[\frac{1}{n-m}\sum_{\begin{subarray}{c}\eta_{i}^{ S}<1-\alpha,\\ S_{i}=\hat{g}(a_{i},x_{i})\end{subarray}}(\alpha+\eta_{i}^{S})f(a_{i},x_{i}) \right].$ (92) $$

Since $f$ is positive and $\eta_{i}\in[-\alpha,1-\alpha]$, it follows

$$ $\displaystyle\mathbb{E}[f(a_{n+1},x_{n+1})(\mathbbm{1}_{[Y(a^{\ast})\in C(X_{n +1},a^{\ast})]}-(1-\alpha))]\geq 0$ (93) $$

and thus

$$ $\displaystyle P_{f}(Y(a^{\ast})\in C_{n+1}C(X_{n+1},A^{\ast}(X_{n+1})))\geq 1-\alpha.$ (94) $$

∎

###### Lemma 5 ( Gibbs et al. ( 2023 ) ) .

### C.3 Proof of Theorem 2

We follow the same outline as in the proof of Theorem 1 in Section C.2. In Lemma 2, we motivated the functional class of distribution shifts. Therefore, it is left to prove the coverage guarantee of $C_{n+1}$.

Key to our proof is the following lemma.

###### Lemma 6 .

The mapping $S\mapsto v_{n+1}^{S}$ is non-increasing in S for all $g^{S}(x,a)$ minimizing

$$ $\displaystyle\min_{g\in\mathcal{F}}\quad$ $\displaystyle\sum_{i=m+1}^{n+1}(1-\alpha)u_{i}+\alpha v_{i}$ (95) s.t. $\displaystyle S_{i}-g(x_{i},a_{i})-u_{i}+v_{i}=0,\quad\forall i=m+1,\ldots,n+1$ $$

for non-negative function classes $\mathcal{F}$ and imputed $S_{n+1}=S$ stemming from a non-negative non-conformity score function (e.g., the residual of the prediction).

###### Proof.

Assume for contradiction that there exists $\tilde{S}>S$ such that $v^{\tilde{S}}_{n+1}>v^{S}_{n+1}$. Then

$$ $\displaystyle(\tilde{S}-S)(v^{\tilde{S}}_{n+1}-v^{S}_{n+1})>0.$ (96) $$

We observe that

$$ $\displaystyle\tilde{S}(S-g^{S}(x_{n+1},a_{n+1})-u^{S}_{n+1}+v^{S}_{n+1})=S( \tilde{S}-g^{\tilde{S}}(x_{n+1},a_{n+1})-u^{\tilde{S}}_{n+1}+v^{\tilde{S}}_{n+ 1})=0.$ (97) $$

Reformulating the equation above yields

$$ $\displaystyle(\tilde{S}-S)(v^{\tilde{S}}_{n+1}-v^{S}_{n+1})$ (98) $\displaystyle=\quad$ $\displaystyle\tilde{S}u^{S}_{n+1}-Su^{\tilde{S}}_{n+1}+\tilde{S}g^{S}(x_{n+1}, a_{n+1})-Sg^{\tilde{S}}(x_{n+1},a_{n+1})+\tilde{S}v^{\tilde{S}}_{n+1}-Sv^{S}_{ n+1}$ (99) $\displaystyle<\quad$ $\displaystyle S(u^{S}_{n+1}-u^{\tilde{S}}_{n+1}+g^{S}(x_{n+1},a_{n+1})-g^{ \tilde{S}}(x_{n+1},a_{n+1})-(v^{S}_{n+1}-v^{\tilde{S}}_{n+1}))$ (100) $\displaystyle=\quad$ $\displaystyle S(S-\tilde{S}).$ (101) $$

This is equivalent to

$$ $\displaystyle(S-\tilde{S})(v^{S}_{n+1}-v^{\tilde{S}}_{n+1})<S(S-\tilde{S})$ (102) $\displaystyle\iff\quad$ $\displaystyle v^{S}_{n+1}-v^{\tilde{S}}_{n+1}>S\geq 0,$ (103) $$

which contradicts the assumption that $v^{\tilde{S}}_{n+1}>v^{S}_{n+1}$.
∎

Coverage guarantees. As in C.2, we fix some $\sigma>0$ and $c_{a}\in[\frac{1}{M},M]$ to obtain a specific $f(a,x):=\frac{c_{a}}{\sqrt{2\pi}\sigma}\frac{\exp{\left(-\frac{(a_{i}-a^{\ast}
)^{2}}{2\sigma^{2}}\right)}}{\hat{\pi}(a_{i}\mid x_{i})}$. We further denote $\hat{g}(a,x)\in\mathcal{F}$ the optimal solution given by the optimal values $\hat{\sigma}$ and $\hat{c_{a}}$.

With the definition of $S^{\ast}$ as the minimum $S$ such that $v_{n+1}^{S^{\ast}}=0$ and Lemma 6, we now can state

$$ $\displaystyle\mathbb{E}[f(a_{n+1},x_{n+1})(\mathbbm{1}_{[Y(a^{\ast})\in C_{n+1 }]}-(1-\alpha))]$ $\displaystyle=\mathbb{E}[f(a_{n+1},x_{n+1})(\mathbbm{1}_{[v_{n+1}^{S}>0]}-(1- \alpha))]$ (104) $\displaystyle=\mathbb{E}[f(a_{n+1},x_{n+1})(\alpha-\mathbbm{1}_{[v^{S}_{n+1}=0 ]})]$ (105) $$

and, by exchangeability of $(f(a_{i},x_{i}),\hat{q}_{S}(a_{i}.x_{i}),S_{i})$, that

$$ $\displaystyle\mathbb{E}[f(a_{n+1},x_{n+1})(\alpha-\mathbbm{1}_{[v^{S}_{n+1}=0] })]$ $\displaystyle=\mathbb{E}\left[\frac{1}{n-m}\sum_{i=m+1}^{n+1}f(a_{n+1},x_{n+1} )(\alpha-\mathbbm{1}_{[v^{S}_{n+1}=0]})\right]$ (106) $\displaystyle=\frac{1}{n-m}\mathbb{E}\left[\sum_{v^{S}_{i}>0}\alpha f(a_{i},x_ {i})-\sum_{v^{S}_{i}=0})(1-\alpha)f(a_{i},x_{i})\right]$ (107) $\displaystyle=\frac{1}{n-m}\mathbb{E}\left[\sum_{S_{i}<\hat{g}(a_{i},x_{i})} \alpha f(a_{i},x_{i})-\sum_{S_{i}\geq\hat{g}(a_{i},x_{i})})(1-\alpha)f(a_{i},x _{i})\right].$ (108) $$

Deriving the Lagrangian above wrt. $f$ yields the stationarity condition

$$ $\displaystyle 0$ $\displaystyle\overset{!}{=}\sum_{i=m+1}^{n+1}\eta_{i}^{S}f(a_{i},x_{i})$ (109) $\displaystyle=\sum_{S_{i}<\hat{g}(a_{i},x_{i})}\eta_{i}^{S}f(a_{i},x_{i})+\sum _{S_{i}>\hat{g}(a_{i},x_{i})}\eta_{i}^{S}f(a_{i},x_{i})+\sum_{S_{i}=\hat{g}(a_ {i},x_{i})}\eta_{i}^{S}f(a_{i},x_{i}).$ (110) $$

The complementary slackness Karush-Kuhn-Tucker conditions yield

$$ $\displaystyle\eta_{i}^{S}\in\begin{cases}-\alpha,&\text{if }S_{i}<\hat{g}(a_{i },x_{i}),\\ [-\alpha,1-\alpha],&\text{if }S_{i}=\hat{g}(a_{i},x_{i}),\\ 1-\alpha,&\text{if }S_{i}>\hat{g}(a_{i},x_{i}).\end{cases}$ (111) $$

Therefore, we receive

$$ $\displaystyle\mathbb{E}[f(a_{n+1},x_{n+1})(\alpha-\mathbbm{1}_{[v^{S}_{n+1}=0] })]$ $\displaystyle=\frac{1}{n-m}\mathbb{E}\left[\sum_{\eta^{S}_{i}<1-\alpha}\alpha f (a_{i},x_{i})-\sum_{\eta^{S}_{i}=1-\alpha})(1-\alpha)f(a_{i},x_{i})\right]$ (112) $\displaystyle=\frac{1}{n-m}\mathbb{E}\left[\sum_{\begin{subarray}{c}\eta_{i}^{ S}<1-\alpha,\\ S_{i}=\hat{g}(a_{i},x_{i})\end{subarray}}(\alpha+\eta_{i}^{S})f(a_{i},x_{i}) \right].$ (113) $$

Since $f$ is positive and $\eta_{i}\in[-\alpha,1-\alpha]$, it follows

$$ $\displaystyle\mathbb{E}[f(a_{n+1},x_{n+1})(\mathbbm{1}_{[Y(a^{\ast})\in C(X_{n +1},a^{\ast})]}-(1-\alpha))]\geq 0$ (114) $$

and thus

$$ $\displaystyle P_{f}(Y(a^{\ast})\in C(X_{n+1},a^{\ast}))\geq 1-\alpha.$ (115) $$

∎

###### Lemma 6 .

###### Proof.

## Appendix D Additional background

### D.1 Extended literature review

In the following, we present related work on CP for causal quantities in more detail.

Recently, Alaa et al. () provided predictive intervals for CATE meta-learners under the assumption of full knowledge of the propensity score. As an extension, Jonkers et al. () proposed a Monte-Carlo sampling approach to receive less conservative intervals. Chen et al. () provide prediction intervals for counterfactual outcomes. However, the proposed method requires access to additional interventional data and is thus not applicable to real-world applications on observational data. All methods are restricted to binary treatments.

Other works focus on prediction intervals for off-policy prediction (Taufiq et al., ; Zhang et al., ) and conformal sensitivity analysis (Yin et al., ), thus neglecting estimation errors arising from propensity or weight estimation or for randomized control trials (Kivaranovic et al., ; Lei & Candès, ). Wang et al. () constructed intervals with treatment-conditional coverage of discrete treatments.

Aiming for group-conditional coverage, Wang et al. () adapted CP to cluster randomized trials. Nevertheless, the method only applies to a finite number of treatments and thus is not applicable to continuous treatments.

Overall, no method can provide exact intervals for continuous treatments. Especially, no method considers the error arising from propensity estimation in the analysis.

### D.2 A note on challenges and difficulties in CP for causal effects of continuous treatments

Existing works on conformal prediction for binary or (low-dimensional) discrete treatment are commonly based on (a) weighted conformal prediction (Tibshirani et al., ) or (b) conformal prediction local coverage guarantees (Lei & Wasserman, ). The first approach provides marginal coverage under a distribution shift through reweighting. It requires computing the weights based on the probability of treatment $A=a$. However, for continuous treatments, this is always zero. Although applicable to binary or low-dimensional discrete treatments (e.g., Lei & Candès, ), this weighting approach cannot be extended similarly to continuous treatments. Furthermore, the propensity of a continuous treatment given by the Dirac delta function $\delta_{a}$ would require us to restrict the calibration to data samples of the specific treatment, which are extremely rare or even *might be missing*. Therefore, the calibration step cannot be employed in our setting. The second approach provides treatment group-conditional coverage. Although again possible for binary or low-dimensional treatments, this approach *does not apply to continuous treatment*s as no treatment groups can be defined. Instead, we propose a novel method for conformal predictions that circumvents the above problems and is carefully tailored to continuous treatments.

### D.3 Causal effects of continuous treatments & kernel smoothing

Causal inference becomes challenging with continuous treatments primarily due to the infinite number of potential outcomes per sample, from which only one outcome is observed. Continuous treatments thus result in causal effects that are generally represented by curves (called dose-response curves) (Kennedy.2018). This is unlike binary treatment, where the causal effects are represented by a single discrete value.

For continuous treatments, the dose-response curves are typically assumed to fulfill some smoothness criterion (e.g., Patrick Schwab et al., ; Schweisthal et al., ). Hence, when estimating treatment effects, interpolation and kernel smoothing of the outcome function are commonly employed (e.g., Kennedy.2018; Nagalapatti et al., ).

Underlying causal estimation with continuous treatments is the generalized propensity score (Imbens, ). It is defined as the conditional probability of receiving treatment $a^{\ast}$ given the covariates $X$ under the following regularity conditions: (i) For each $i$, $Y_{i}(a),x_{i},A_{i}$ are defined on a common probability space; (ii) $A_{i}$ is continuously distributed with respect to the Lebesgue measure; and (iii) $Y_{i}=Y_{i}(A_{i})$ is a well-defined random variable.

Approximating the density $\delta_{a^{\ast}}(a)$ of the hard intervention $a^{\ast}$ through a Gaussian kernel follows directly from the definition of $\delta_{a^{\ast}}(a)$ as the limit of such kernel. This is also common in the literature (e.g., Kallus & Zhou, ). Importantly, we note that we do not directly approximate the potential outcome $Y(a^{\ast})$ (but only the propensity scores). Thus, we do not have a bias-variance trade-off of the estimated outcome. Due to the smoothness of the dose-response curve, it is now valid to employ observed samples within a treatment region of $a^{\ast}$ defined by $\sigma$ to construct the intervals. We note that the importance of the samples is weighted by the inverse distance of the sample to $a^{\ast}$ in treatment space. We give further intuition on the relationship between $\sigma$, the importance of observational samples, and the prediction interval width in the following.

### D.4 Interpretation of optimal parameters

To obtain CP intervals under an unknown distribution shift, we approximate the Dirac-delta distribution representing the hard intervention by a Gaussian function as

$$ $\displaystyle\delta_{a^{\ast}}(a)=\lim\limits_{\sigma\rightarrow 0}\frac{1}{ \sqrt{2\pi}\sigma}\exp\left(-\frac{(a-a^{\ast})^{2}}{2\sigma^{2}}\right).$ (116) $$

In Theorem 2, we thus optimize over $\sigma>0$ and $c_{a}\in[\frac{1}{M},M]$ to obtain the $(1-\alpha)$-quantile of the distribution shift-calibrated non-conformity scores. The optimal parameter $\sigma^{\ast}$ represents a trade-off between the uncertainty in the prediction and the uncertainty in the interval construction: A small $\sigma^{\ast}$ resembles the propensity of the hard intervention best. Thus, with sufficient or even infinite data close to $a^{\ast}$, we could construct the narrowest CP interval. However, the smaller $\sigma$, the less data close to $a^{\ast}$ will be available to calculate the prediction interval in practice. As a result, many calibration data samples will be strongly perturbed during the calculation, which increases the uncertainty and, thus, the interval size.

## Appendix E Experimentation details

### E.1 Synthetic dataset generation

We consider two different propensity and outcome functions. In each setting, we assign two types of interventions: a known propensity shift of $\Delta=1,5,10$, i.e., three soft interventions $a^{\ast}=a+\Delta$, and the point interventions $a^{\ast}\in\{1x,5x,10,\}$ given the confounder $X=x$.

We generate synthetic datasets for each setting. Specifically, we draw each 2000 train, 1000 calibration, and each 1000 test samples per intervention from the following structural equations. Dataset 1 is given by

$$ $\displaystyle X\sim\mathrm{Uniform}[1,4]\quad\text{(integer)}$ $\displaystyle A\sim p\cdot\mathrm{Uniform}[0,5X)+(1-p)\mathrm{Uniform}[5X,40], \quad p\sim\mathrm{Bernoulli}(0.3)$ $\displaystyle Y\sim\sin\left(\frac{\pi}{6}(0.1A-0.5X)\right)+\mathrm{Normal}(0 ,0.1),$ $$

and dataset 2 by

$$ $\displaystyle X\sim\mathrm{Uniform}[1,4]\quad\text{(integer)}$ $\displaystyle A\sim\mathrm{Normal}(5X,10)$ $\displaystyle Y\sim\sin\left(\frac{\pi}{2}(0.1A-0.1X)\right)+\mathrm{Normal}(0 ,0.1),$ $$

### E.2 Medical dataset

We use the MIMIC-III dataset (Johnson et al., ), which includes electronic health records (EHRs) from patients admitted to intensive care units. From this dataset, we extract 8 confounders (heart rate, sodium, blood pressure, glucose, hematocrit, respiratory rate, age, gender) and a continuous treatment (mechanical ventilation) using an open-source preprocessing pipeline (Wang et al., ). From each patient trajectory in the EHRs, we sample random time points and average the value of each variable over the ten hours before the sampled time point. We define the variable blood pressure after treatment as the outcome, for which we additionally apply a transformation to be more dependent on the treatment and less on the blood pressure before treatment. We remove all patients (samples) with missing values and outliers from the dataset. Outliers are defined as samples with values smaller than the 0.1th percentile or larger than the 99.9th percentile of the corresponding variable. The final dataset contains 14719 samples, which we split into train (60%), val (10%), calibration (20%), and test (10%) sets.

### E.3 Implementation details

Our experiments are implemented in PyTorch Lightning. We provide our code in our [GitHub](https://github.com/m-schroder/ContinuousCausalCP) repository.

We limited the experiments to standard multi-layer perception (MLP) regression models, consisting of three layers of width 16 with ReLu activation function and MC dropout at a rate of 0.1, optimized via Adam. We did not perform hyperparameter optimization, as our method aimed to provide an agnostic prediction interval applicable to any prediction model. All models were trained for 300 epochs with batch size 32.

Our algorithm requires solving (non-convex) optimization problems through mathematical optimization. We chose to employ two interior-point solvers in our experiments: For the experiments with soft interventions that pose convex optimization problems, we use the solver MOSEK. For the hard interventions, which included non-convex problems, we used the solver IPOPT. Both solvers were run with default parameters.

## Appendix F Further results

We present further results from our experiments in Section 5. Specifically, we state the prediction performance of the underlying models $\phi$, discuss the scalability of our approach, and show the prediction intervals per covariate for various significance levels $\alpha$ and soft interventions $\Delta$ of our synthetic experiments on dataset 1 and dataset 2.

Performance:
We first report the performance of the underlying prediction models $\phi$ for the synthetic datasets across 10 runs in Table 2. The prediction model on the real-world dataset achieved a mean squared error loss of 1.2373.

**Table 2: Mean and standard deviation of MSE loss of prediction models $\phi$ across 10 runs.**
|  | Dataset 1 | Dataset 2 | MIMIC |
| --- | --- | --- | --- |
| $\phi$ | 0.0216 (0.0056) | 0.9029 (0.3908) | 0.0141 (0.0057) |
| Ens. | 0.0094 (2.1169$e^{-5}$) | 0.0130 (0.0003) | - |

Figure: Figure 8: Prediction intervals for multiple significance levels $\alpha$ for the synthetic dataset 1 with intervention $\Delta=1$.
Refer to caption: x8.png

Scalability:
Calculating the prediction intervals requires an iterative search for an optimal value $S^{\ast}$. Therefore, the underlying optimization problem must be fitted multiple times throughout the algorithm, potentially posing scalability problems. In our empirical studies, however, we did not encounter scalability issues. Importantly, we found that the average runtime of our algorithm on a standard desktop CPU is only 16.43 seconds.

Prediction intervals:
In Figures 8, 10, and 10, we present the prediction bands given by our method and MC dropout on dataset 1. In particular, for confounder $X=1$, our method shows a large increase in the uncertainty in the potential outcomes of treatments affected by the intervention. MC dropout does not capture this uncertainty.

Figure: Figure 9: Prediction intervals for multiple significance levels $\alpha$ for the synthetic dataset 1 with intervention $\Delta=5$
Refer to caption: x9.png

In Figures 12 and 12, we present the prediction bands given by our method and MC dropout on dataset 2 for the soft interventions $\Delta=1$ and $\Delta=10$ (the results for $\Delta=5$ were presented in the main paper). We observe that the prediction intervals for $\Delta=10$ become extremely wide for high treatments. This aligns with our expectation, as data for high treatments in combination with low confounders is rare or even absent in the dataset. Thus, the expected uncertainty is very high.

Figure: Figure 11: Prediction intervals for multiple significance levels $\alpha$ for the synthetic dataset 2 with intervention $\Delta=1$
Refer to caption: x11.png

## Peer Reviews

confidence: 2
contribution: 3
initial_rating: 6
presentation: 3
questions: See weakness.
soundness: 3
strengths: Clear writing & solid theoretical results
summary: This paper proposes a conformal prediction method for continuous treatments, addressing the challenge of unknown propensity scores. It derives finite-sample prediction intervals, provides an algorithm for calculating them, and demonstrates effectiveness on synthetic and real datasets. This is the first method for continuous treatments with unknown propensity scores.
weaknesses: The reviewer is not an expert in conformal prediction for causal inference but would assume there are working targeting the studies setting: uncertainty quantification of causal effects for continuous treatment.

Equation 14 seems to use the kernel function to approximate the indicator function. Does this imply some smoothness assumptions?

The optimal parameter $\sigma^*$ represents a trade-off between the uncertainty in the prediction and the uncertainty in the interval construction. It is similar to the bandwidth approximation of the indicator function using a kernel function. Are there any empirical values ​​here, such as \sigma = cn^{-2/5} to  trade-off between the prediction and interval.

The method is demonstrated on specific datasets and scenarios. While the authors claim that their method scales to high-dimensional settings, further validation on a broader range of datasets and treatment types would strengthen the generalizability of the findings.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: 1. Major question is in **Weakness** part 
2. I think the real data in the paper can not justify the superiority of authors' method well as I can see at range [0.6, 1], it yields not very reasonable prediction (blood pressure can be less than 0). 
Also, the interval estimate is too wide to be meaningful.
soundness: 3
strengths: 1. They study a very important question in the causal inference and provide novel methods. 
2. They provide the finite-sample theorem for their methods. 
3. They clearly and accessibly present the research question and methods, using figures and concrete examples.
summary: In the paper, the authors propose a Conformal-prediction (CP) based method to obtain the interval estimation of the potential outcomes under continuous treatments. 
To achieve this, the authors follows (Gibbs et al., 2023; Romano et al., 2019) and reformulated split CP as an augmented quantile regression.
The authors derive finite-sample prediction intervals for cases where the propensity score function is known and unknown.
Moreover, they provide finite-sample theoretical guarantees.
weaknesses: The authors do not have a good evaluation of their method in terms of numerical studies. 

Firstly, the authors only focus on the empirical coverage, and claim that it should be the higher the better. 
In my opinion, it seems not very reasonable. 
If we simply choose the [-inf, inf] as the interval, it is always has the 100% empirical coverage. 
A more balanced evaluation should consider both the interval length and empirical coverage. Additionally, the empirical coverage should closely match the nominal coverage for optimal results.


When checking the numerical results, I do believe the method in the implementation is not very stable under some cases. 
For example, in Fig 4 ($\alpha=0.2$, 5), with only 50 repetitions, I can see many outliers in the boxplot, indicating there may be some issues there. Similar issue can be found in Fig 5 ($\alpha=0.1$, 1) and  ($\alpha=0.2$, 5). 
Anyway, I know such unstableness may be inevitable for a complicated method, I think the authors should have some discussions on this issue.

---

confidence: 4
contribution: 3
initial_rating: 5
presentation: 2
questions: * Consider building semi-synthetic experiments for stronger results.
 * How can $M$ of Assumption 1 be calibrated in practice, for tightness?
 * Do you assume some kind of continuity in the dose response?
soundness: 2
strengths: The problem of conformal prediction of dose responses is clearly significant. The authors proposed an original formulation that uses Gaussian smoothing over the treatments and finds an optimal bandwidth for each problem instance. It is also important to consider soft interventions (by the authors' definition) in addition to hard interventions, since it is common in medicine and other settings to study incremental changes.
summary: Conformal prediction for causal estimands is useful because it provides finite-sample coverage guarantees. Predicting potential outcomes is challenging due to the potentially unknown distribution shift between the training set and the test set, which is subject to intervention. This work extends conformal prediction under domain shift to address the challenges that are introduced when the treatment variable is continuous, and the causal estimand is a dose response.
weaknesses: * I do not understand some core elements of the motivation. What is the intention behind the separation of challenge (a) and challenge (b)? Estimating the propensity score is more difficult with continuous treatments, but it is still a challenge with discrete treatments. By my understanding, the main challenge for bringing conformal prediction to continuous treatments is that you rarely see data points with the same treatment value, so you have to do some sort of smoothing. The motivating challenges do not mention this (unless I am misreading one of them).
 * I also found the roles of soft and hard interventions confusing. Under "Scenarios" in Section 4, the authors explain that soft interventions are relevant when the propensity is known, while hard interventions are reserved for when the propensity is unknown. This storyline is convenient for the mathematical derivations. However, I am not convinced: even if the propensity function is unknown, it can still be of interest to estimate the effect of increasing the treatment by a small amount on top of the existing policy. These cases have been studied before:
   * Many papers by Chernozhukov et al. that mention "average causal derivatives"
   * Marmarelis et al. "Policy Learning for Localized Interventions from Observational Data" in AISTATS 2024.
 * When the propensity / policy is unknown, the validity of the conformal prediction relies on a constant $M$ that bounds the error of the estimated propensity. This resembles a sensitivity analysis. It appears difficult to know which $M$ to use.
 * The results from the medical dataset do not seem convincing. It is good to see that intervals are wider in treatment regions that are less common in the data. However, the fact that intervals are wider than MC-dropout is not necessarily a good thing.

Minor notes.
 * "Quantile regression might yield non-unique solutions, so we later restrict the analysis to solvers invariant to the data ordering". Please expand on this sentence. How does ordering invariance help?
 * The algorithm text is small and quite dense.
 * An interpretation of the optimized $\sigma$ is provided in the supplement, but further discussion would be helpful. Please clarify its role in the optimization problem.

---

confidence: 3
contribution: 3
initial_rating: 3
presentation: 3
questions: - What is your added contribution to the work of Gibbs et al. (2023), which you leverage in this work?
- Can you explain the variability in empirical coverage across different significance levels, particularly for α = 0.1? Do you plan to increase the number of runs or experiments to solidify your claims?
- Could you indicate how long it takes to run these experiments?
- I think there is a typo on line 476 that should be Table 1, right?
- In your synthetic experiments, why are the hard interventions related to x and not just hard interventions?
- Related to Assumption 1 (line 295-299), does this assumption mean that the estimation error needs to be the same for each X given an intervention $a$, otherwise should it not be $c_a(x)$ instead of $c_a$.
- Related to Theorem 2, what is the meaning of resulting $\sigma$ and $c_a$ of the optimization problem.
- Related to Lemma 4 (lines 712-736):
    - I think you need to clarify that this is only relevant for absolute residual non-conformity score.
    - Next, I think the proof is complete. You do not prove the connection between the Bonferonni corrected significance level and the targeted significance level.
    - In my opinion, the proof for this lemma is even elementary since it is just an application of the Bonferonni correction.
soundness: 2
strengths: - **Clear Motivation:** The paper's objective is well-placed within causal machine learning, particularly in high-stakes areas like medicine, where uncertainty quantification is critical for decision-making.
- **Novel Contribution:** The authors propose a novel method for deriving conformal prediction intervals in the presence of continuous treatments and unknown propensity scores. This is a significant contribution, as most prior work focuses on binary or discrete treatments.
summary: The paper applies the framework introduced by Gibbs et al. (2023) to produce prediction intervals, with a finite-sample coverage guarantee, for potential outcomes of continuous treatments. Gibbs et al. (2023) can provide certain types of conditional coverage guarantees by reformulating the coverage as coverage over a class of covariate shifts. Gibbs et al. (2023) already showed that given an appropriate choice of $\mathcal{F}$, that their prediction set achieves coverage under covariate shift, so long that the covariate shift can be represented as, $P_X$ titled by $f$, and $f \in \mathcal{F}$.

This work postulates the causal effect of continuous treatments as soft and hard interventions.  A soft intervention can be seen as a value shift in the generalized propensity score. A hard intervention can be seen as shifting the generalized propensity score to a Dirac delta function to the desired treatment value. In the paper, they define an appropriate function class for the different interventions and deal with their specific covariate shifts. For the hard intervention setting, they also ensure that the function class can deal with an estimated generalized propensity score with a bounded estimation error of the propensity function. 

They provide some experiments on synthetic datasets (from their own) and medical ICU datasets, benchmarking their prediciton intervals against an MC dropout and Ensemble approach.

[1] I. Gibbs, J. J. Cherian, and E. J. Candès, “Conformal Prediction With Conditional Guarantees,” Dec. 20, 2023, arXiv: arXiv:2305.12616. doi: 10.48550/arXiv.2305.12616.
weaknesses: - **Terminology Issues:** Terms such as "finite-sample predictive intervals" should be revised to "finite-sample validity guarantees" to align with established terminology in the field. Additionally, (line 090) “strong finite-sample mathematical guarantees” should be “finite-sample validity guarantees”. Similarly, the frequent use of the term "exact" in relation to conditional coverage is misleading, as exact conditional coverage is known to be unattainable. I also assume you do not provide exact coverage in the context of exact and conditional coverage.
- **Relation to Lei and Candes (2021):** In Figure 2, and the appendix, the paper claims that the work of Lei and Candes can not deal with an unknown propensity, this is not true it defines and proves a coverage gap based on the error in the estimation of the likelihood ratio’s (= propensity in the case of treatment effect) in weighted CP.
- **Relation to Gibbs et al. (2023): T**he paper mentions the paper of Gibbs et al. (2023) in that it reformulates the split CP as an augmented quantile regression and the use of Lemma 1. But throughout the work, several notions from this work are taken, such as the tilting of the propensity function (lines 210-215), Theorem 1 (lines 262-276) and Theorem 2 (lines 352-371) which is heavily inspired by this work as it can even be seen in the borrowed notation of Section 4 of Gibbs et al. (2023), Algorithm 1 similarly with Algorithm 1 of Gibbs et al. (2023).
- **Under-coverage in Empirical Results**: The experimental results show some under-coverage, particularly for specific significance levels (e.g., α = 0.1). This variability raises concerns about the method's robustness, and the number of runs (10) in some experiments is too low to draw reliable conclusions. More runs would improve the statistical robustness of the findings.
- **Computational Complexity**: While the paper introduces an algorithm, the appendix discusses the method's computational efficiency. It points out that complexity can be high due to the optimization solver. This should be included in the limitations. Additionally, you should include your experiments' fitting and inference times in the paper, as this could indicate if computational complexity is an issue in practice.
- **Too simple simulation experiments: T**he proposed synthetic experiments are way too simplistic, with only one discrete confounder, continuous treatment, and continuous outcome. This is not in line with other treatment effect literature; we cannot adequately evaluate the method with these experiments.
- **Baseline**: The choice of baseline methods, particularly MC dropout, may not be ideal for comparison. Since these approaches only quantify epistemic uncertainty, which is probably minimal due to the simple experiment settings, consider a benchmark approach that quantifies both aleatoric and epistemic uncertainty, or at least aleatoric such as the Gaussian process or maximum likelihood estimator. The MLE can be easily combined with MC dropout to get both aleatoric and epistemic to get a reasonable benchmark. Also, you did not evaluate vanilla CP without considering covariate shift as a benchmark.

Meta-review: This paper proposed a conformal prediction method for continuous treatments, addressing the challenge of unknown propensity scores. The problem is well motivated and the extension of conformal inference to continuous treatments is important. However, some main concerns remain:

 - The novelty and contribution compared with Gibbs et al. (2023) is limited, and theoretical results largely follows this paper.

 - There are concerns raised about the method's robustness and the empirical performance. For example, the authors report zero coverage for vanilla CP and GP methods, which suggests the presence of purely epistemic uncertainty and extremely heteroscedastic uncertainty and seems inconsistent with the data-generating processes. Additionally, there is a critical issue with the reported average width of prediction intervals for V-CP (claimed to be 0.0003). This value is mathematically improbable and is several orders of magnitude smaller than these theoretical lower bounds, which raises serious questions about the experimental methodology or reporting.

 - About utilities. The time complexity of the proposed approach is not very practical. The reviewers also suggested to have semi-synthetic benchmarks in the main text.
 
Given the above reasons, after discussion with the reviewers, we agreed it is not quite ready for publication and we'd encourage the authors to take into consideration all the feedback provided by the reviewers to strengthen their manuscript for resubmission.

Decision: Reject