Quantifying Aleatoric Uncertainty of the Treatment Effect: A Novel Orthogonal Learner

Authors: valentyn  melnychuk
stefan  feuerriegel
mihaela  van der schaar
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-13 10:59:02

Abstract: Estimating causal quantities from observational data is crucial for understanding the safety and effectiveness of medical treatments. However, to make reliable inferences, medical practitioners require not only estimating averaged causal quantities, such as the conditional average treatment effect, but also understanding the randomness of the treatment effect as a random variable. This randomness is referred to as aleatoric uncertainty and is necessary for understanding the probability of benefit from treatment or quantiles of the treatment effect. Yet, the aleatoric uncertainty of the treatment effect has received surprisingly little attention in the causal machine learning community. To fill this gap, we aim to quantify the aleatoric uncertainty of the treatment effect at the covariate-conditional level, namely, the conditional distribution of the treatment effect (CDTE). Unlike average causal quantities, the CDTE is not point identifiable without strong additional assumptions. As a remedy, we employ partial identification to obtain sharp bounds on the CDTE and thereby quantify the aleatoric uncertainty of the treatment effect. We then develop a novel, orthogonal learner for the bounds on the CDTE, which we call AU-learner. We further show that our AU-learner has several strengths in that it satisfies Neyman-orthogonality and is doubly robust. Finally, we propose a fully-parametric deep learning instantiation of our AU-learner.

Main paper: 
 introduction: 
 Estimating causal quantities from observational data is crucial for decision-making in medicine [9, 12, 22, 30, 70] . For example, medical practitioners are interested in estimating the effect of chemotherapy vs. immunotherapy on patient survival from electronic health records to understand the best treatment strategies in cancer care. Here, common estimation targets are averaged causal quantities such as the average treatment effect (ATE) and the conditional average treatment effect (CATE), yet averaged causal quantities do not allow for understanding the variability of the treatment effect.What is needed for the reliability of causal quantities in medicine? To obtain reliable causal quantities, one often needs to "move beyond the mean" [44, 68] and consider the inherent randomness in the treatment effect as a random variable. This randomness is referred to as aleatoric uncertainty [17, 60, 110] . Quantifying the aleatoric uncertainty of the treatment effect is relevant in medical practice to understand the probability of benefit from treatment [26, 60] and the quantiles and variance of the treatment effect [5, 17, 26, 33, 59] . As an example, averaged quantities such as the CATE would simply suggest a positive effect for some patients, while the probability of benefit from treatment can inform patients about the odds of being negatively affected by the treatment. Hence, aleatoric uncertainty of the treatment effect promises additional, fine-grained insights beyond simple averages.Methods for quantifying the aleatoric uncertainty of the treatment effect have gained surprisingly little attention in the causal machine learning community. So far, machine learning for treatment Our main contribution relates to the estimation, shown in yellow . However, moving from CATE identification and estimation to our setting comes with important challenges: 1 CATE (shown in green) is point identifiable but the CDTE is not (shown in blue); 2 there is no closed-form expression of the target estimand in terms of nuisance functions and, because of that, CATE learners cannot be directly adapted for estimation; and 3 CATE is an unconstrained target estimand whereas Makarov bounds (shown in gray ) are monotonous and contained in the interval [0, 1].effect estimation was primarily focused on estimating averaged causal quantities [20, 57, 64, 69, 96, 98, 114, 125, 128] . Some research aims to quantify the epistemic uncertainty in treatment effect estimation [51] or the total uncertainty (but without distinguishing the types of uncertainty) [1, 67, 76] . Other works focused on the aleatoric uncertainty of the potential outcomes [13, 31, 65, 66, 94] 1 or on contrasts between distributions of potential outcomes (also known as distributional treatment effects) [16, 29, 62, 92, 100] . 2 However, to the best of our knowledge, there is no comprehensive meta-learning theory for the estimation of the aleatoric uncertainty in the treatment effect.In this paper, we aim to quantify the aleatoric uncertainty of the treatment effect at the covariateconditional level in the form of a conditional distribution of the treatment effect (CDTE). Knowing the CDTE would automatically allow one to compute the above-mentioned quantities of aleatoric uncertainty, namely, the probability of benefit from treatment and the quantiles and variance of the treatment effect, at both population and covariate-conditional levels.Yet, the identification and estimation of the CDTE in contrast to CATE come with three challenges as follows (see Fig. 1 ):Challenge 1 is that the CDTE does not allow for point identifiable, neither in the potential outcomes framework nor in randomized control trials due to the fundamental problem of causal inference as counterfactual outcomes can not be observed [26, 88] . We thus employ partial identification [48] to obtain bounds on the CDTE and thereby quantify the aleatoric uncertainty of the treatment effect. Specifically, we focus on Makarov bounds [87, 126, 129] that give sharp bounds for both the cumulative distribution function (CDF) and the quantiles of the CDTE.Challenge 2 is that there is no closed-form expression of the target estimand in terms of nuisance functions. Because of this, existing CATE learners cannot be directly adapted to our task of estimating Makarov bounds. For example, there are no orthogonal learners in the general setting, and existing approaches only use naïve plug-in estimators/learners. Furthermore, even the derivation of the orthogonal loss is non-trivial as there is no efficient influence function at hand for the Makarov bounds.Challenge3 is that CATE is an unconstrained target estimand whereas Makarov bounds are monotonous and contained in the interval [0, 1]. Notably, any constraints of the target estimand could be violated by orthogonal learners [71, 125] . Therefore, an orthogonal learner for Makarov bounds needs to be carefully adapted, especially to perform well in low-sample settings. Both upper and lower plots have the same total uncertainty but vastly different aleatoric and epistemic components. Yet, aleatoric uncertainty is nonidentifiable (see Challenge 1 ).In this paper, we develop a novel, orthogonal learner for estimating Makarov bounds which we call AU-learner, which allows to quantify the aleatoric uncertainty of the treatment effect. Our AU-learner addresses all of the above-mentioned challenges 1 -3 . Further, our AUlearner has several useful theoretical properties, such as satisfying Neyman-orthogonality and double robustness. Finally, we propose a flexible, fully-parametric deep learning instantiation of our AU-learner. For this, we make use of conditional normalizing flows and call our method AU-CNFs.To summarize, our contributions are as follows: 3 1. We derive a novel, orthogonal leaner called AU-learner to quantify the aleatoric uncertainty of the treatment effect. For this, we estimate Makarov bounds on the CDF/quantiles of the conditional distribution of the treatment effect (CDTE). 2. We prove several favorable theoretical properties of our AU-learner, such as Neyman-orthogonality and double robustness. 3. We propose a flexible deep learning instantiation of our AU-learner based on conditional normalizing flows, which we call AU-CNFs, and demonstrate its effectiveness over several benchmarks. 
related work: 
 In the following, we briefly summarize the existing works on uncertainty quantification in the potential outcomes framework; on the identification of the CDTE; and on the estimation of Makarov bounds. For a more detailed overview of literature, we refer to Appendix A.Uncertainty quantification in the potential outcomes framework. The (total) uncertainty of a predictive model in machine learning is generally split into (a) epistemic and (b) aleatoric uncertainty [41, 50] . 4 This split is important, as it informs a decision-maker about the source of uncertainty (see Fig. 2 ), especially in the context of the potential outcomes framework. (a) Epistemic uncertainty was studied for predictive models targeting at identifiable averaged causal quantities, such as conditional average potential outcomes (CAPOs) and CATE [51, 52] . (b) Aleatoric uncertainty, on the other hand, is only identifiable for potential outcomes [88] . Prominent methods focus on interventional (counterfactual) quantities such as: (i) CDF/quantiles estimation [6, 13, 31, 43, 86] ; (ii) density estimation [65, 66, 91, 94, 97, 102, 124] ; and (iii) distributional distances (also known as distributional treatment effects) estimation [16, 29, 62, 92, 100 ]. Yet, our work differs substantially from the above, as we aim at inferring the aleatoric uncertainty of the treatment effect, which is only partially identifiable. Identification of the distribution of the treatment effect. Point identification of the distribution of the treatment effect (or, equivalently, a joint distribution of potential outcomes) is only possible under additional assumptions on the data-generating mechanism. A common example is, e. g., invertibility of latent outcome noise [2, 3, 6] . Other works have rather focused on partial identification [48] . For example, [88] proposed sharp bounds for the distribution of the treatment effect under a monotonicity assumption. Later, assumption-free sharp bounds were proposed for both the joint CDF of potential outcomes [26] and for the variance of treatment effect [5] , both known as Fréchet-Hoeffding bounds [39, 49] . Finally, [26, 28, 33] proposed sharp bounds on the CDF/quantiles of the treatment effect without any additional assumptions, so-called Makarov bounds [87, 126] . Makarov bounds were further generalized [81, 129] and applied to other settings [17, 25, 27, 36, 59] but different from ours.Estimation of Makarov bounds. [26, 110] and for potential outcomes framework [17, 77, 110] . Crucially, these methods are not orthogonal and, thus, are sensitive to the misspecification of the nuisance functions. Nevertheless, we include the latter methods [17, 77] as baselines for our experiments as they use a highly flexible CDF estimator based on kernel density estimators. Some works also developed efficient estimators for Makarov bounds at the population level (analogous to the two-stage orthogonal learners at the covariate-conditional level) but only in highly restricted settings. In particular, [60] is restricted to binary outcomes, [112] assumed a known propensity score, and [54] made special optimization assumptions. 5 In addition, all three works [54, 60, 112] suggest fixing a value of δ/α, which the CDF/quantiles of the treatment effect are evaluated at; when our work suggests targeting at several values of δ/α at once. Therefore, the previous methods are not applicable to our general setting of estimating covariate-conditional level Makarov bounds.Research gap. To the best of our knowledge, we are the first to propose an orthogonal learner for estimating Makarov bounds on the CDF/quantiles of conditional distribution of the treatment effect.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes]Justification: We include a paragraph discussing the limitations of our method in Sec. 7.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
identification of distribution of treatment effect: 
 Notation. Let capital letters X, A, Y, ∆ denote random variables and small letters x, a, y, δ their realizations from domains X , A, Y, ∆. Let P(Z) denote a distribution of some random variable Z, and let P(Z = z) be the corresponding density or probability mass function. Furthermore, π(x) = P(A = 1 | X = x) is propensity score, µ a (x) = E(Y | X = x, A = a) are conditional expectations, and F a (y | x) = P(Y ≤ y | x, a) is a conditional outcome CDF. For other conditional quantities or distributions, we use short forms whenever possible; e. g.,E(Y | x) = E(Y | X = x). Further, P n {f (Z)} = 1 n n i=1 f (z i ) is a sample average of a random f (Z),where n is the sample size. We denote linear rectifier functions as [x] + = max(x, 0) and [x] -= min(x, 0), and sup/inf convolutions of two functions [116] f 1 (• | x), f 2 (• | x) as (f 1 * f 2 ) Y (δ | x) = sup y∈Y {f 1 (y | x) -f 2 (y -δ | x)} and (f 1 * f 2 ) Y (δ | x) = inf y∈Y {f 1 (y | x) -f 2 (y -δ | x)}.Problem setup. We consider the standard setting of the Neyman-Rubin potential outcomes framework [109] . That is, we have an observational dataset D with a binary treatment A ∈ A = {0, 1}, potentially high-dimensional covariates X ∈ X ⊆ R dx and a continuous outcome Y ∈ Y ⊆ R. For instance, a typical scenario is in cancer therapy, where the outcome is tumor growth, the treatment is whether chemotherapy is given, and the covariates include patient details like age and sex. We define a joint random variable Z = (X, A, Y ). D = {x i , a i , y i } n i=1 is sampled i.i.d. from the observational distribution P(Z) = P(X, Y, A), where n is the sample size. The potential outcomes framework then makes three (causal) assumptions, i. e., (1) consistency: if A = a, then Y [a] = Y ; (2) overlap: P(0 ≤ π(X) ≤ 1) = 1; and (3) exchangeability: A ⊥ ⊥ (Y [0], Y [1]) | X.Treatment effect distribution. In this paper, we refer to the treatment effect ∆ = Y [1] -Y [0] as a random variable. 6 The CATE is given by τ (x) = E(∆ | x), which is identifiable as µ 1 (x)µ 0 (x) under the causal assumptions (1)-( 3). We are interested in identifying a conditional distribution of the treatment effect (CDTE), specifically, its CDF or quantiles:EQUATIONThe CDF and the quantiles of the CDTE are point non-identifiable due to the fundamental problem of causal inference, i. e., that the counterfactual outcome, Y [1 -A], is never observed. This is illustrated in Fig. 3 , where both conditional potential outcome distributions are identifiable asP(Y [a] | x) = P(Y | x, a); but a conditional joint distribution, P(Y [0], Y [1] | x), and the CDTE,P(∆ | x), are not.Partial identification of the CDTE. Fan et al. [26] proposed pointwise sharp bounds on the CDF and the quantiles of the CDTE, so-called Makarov bounds [87, 126, 129] . Given that the outcome Y is continuous, the Makarov bounds for the CDF of the CDTE are given by linearly rectified sup/inf convolutions [116] of conditional CDFs of potential outcomes:EQUATIONwhere δ ∈ ∆, andF a (y | x) is the CDF of P(Y [a] | x).Similarly, Makarov bounds can be formulated for the quantiles of the CDTE:EQUATIONEQUATIONwhere α ∈ [0, 1], and F -1 a (u | x) are the quantiles of P(Y [a] | x). Then, under the causal assumptions (1)-( 3), conditional distributions of potential outcomes coincide with observed ones, i. e., P(Y [a]) = P(Y | x, a). Notably, Makarov bounds on the CDFs are CDFs themselves, but these CDFs do not correspond to the solution of the partial identification task (which implies pointwise sharpness). We refer to Appendix B for more illustrations about the inference of the Makarov bounds, the explanation of pointwise sharpness, and Makarov bounds for categorical/mixed-type outcomes.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [Yes]• The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [Yes] Justification: We include a paragraph discussing the broader impacts of our method in Sec. 7. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
an au-learner for estimating makarov bounds: 
 In the following, we develop a theory of orthogonal learning for Makarov bounds, which then gives rise to our AU-learner. For this, we first review the plug-in learner and its shortcomings. Motivated by this, we then derive two-stage learners. Here, we first present a novel CA-learner in an intermediate step and finally our AU-learner. Note that both are novel but we frame our contributions around the AU-learner because of favorable theoretical properties. For notation, we use over-and underlines as in, e. g., F(δ | x) to refer to the upper and/or lower bound.Plug-in learner. A naïve way to construct estimators for Makarov bounds [17, 26, 77, 110] is to estimate conditional outcome CDFs, F a (y | x), and plug them into Eq. (3):EQUATIONThe plug-in learner for the bounds of the quantiles, F PI -1 (α | x), can be obtained in a similar way but where one uses Eq. ( 4) and relies on the estimators of either the conditional outcome quantiles or inverse of the estimated conditional outcome CDFs, F -1 a (u | x). P(Y [1] -Y [0] ≤ 0 | x 7 ) ∈ [0, 0.242].Hence, Makarov bounds are informative almost everywhere (except δ = τ (x 7 )).Shortcomings. The plug-in learner suffers from two important shortcomings [96] . (a) The plug-in learner does not account for the selection bias, meaning that F 1 is estimated better for the treated population and F 0 for the untreated. Hence, it could be necessary to re-weight the loss wrt. to the propensity score. A remedy is to employ an inverse propensity of treatment weighted (IPTW) learner for both F 0 and F 1 [6, 43] (see Appendix C). (b) The plug-in learner does not target the Makarov bounds directly but rather the conditional outcome distributions. Therefore, it is unclear how to incorporate an inductive bias that the Makarov bounds are less heterogeneous than either of the conditional outcome CDFs (i. e., the Makarov bounds can depend on a subset of covariates X). The second shortcoming thus motivates our derivation of two-stage learners.In order to address the above shortcomings of the plug-in learners, a two-stage learning theory was proposed [15, 35, 72, 96 ]. Yet, the two-staged learning theory is primarily built for simple target estimands (e. g., CATE) and therefore requires non-trivial adaptations by us to extend to Makarov bounds, which we do in the following.Working model & target risk. Our two-stage learners seek to find the best approximation of the ground-truth Makarov bounds, functional target estimands, by a (parametric) working model G.Formally, the working model is givenG = {g(δ, x) | g : ∆ × X → [0, 1]; g(•, x)is non-decreasing}.The best approximation g * is then obtained by minimizing a (population) target risk via g * = arg min g∈G L(g). 7 In our setting, L is chosen as some distributional distance between the target estimands and the working model. Specifically, we use continuous ranked probability score (CRPS) [40, 117] as a target risk for learning Makarov bounds on the CDF viaEQUATIONor squared Wasserstein-2 distance as a target risk for learning Makarov bounds on the quantiles viaEQUATION)The target risks in Eq. ( 6) and ( 7) cannot be directly minimized as we do not observe the ground-truth Makarov bounds. Yet, due to the identifiability results in Sec. 3, the Makarov bounds depend on the nuisance functions, F a (y | x), which can be estimated from the observational data. We perform that in the following. From now on, we denote the ground-truth nuisance functions as η and their estimates as η. Also, we make the dependence on the target risks of the nuisance functions explicit; that is, we write the target risk as L(g, η) and the target estimands as F(δ | X; η) and F -1 (α | X; η).Covariate-adjusted learner. A straightforward way to estimate and then minimize the target risk is to plug-in the estimates of conditional outcome CDFs, F a (y | x), into Eq. ( 6) and (7) , respectively. This yields the so-called covariateadjusted (CA) learner, which aims at minimizing the following losses (empirical risks):EQUATIONEQUATIONwhere we call bothF PI (δ | x; η) = F PI (δ | x) and F PI -1 (α | x; η) = F PI -1 (α | x)pseudo-CDFs and pseudo-quantiles, respectively, and where both can be obtained from Eq. ( 5).The CA-learner addresses the shortcoming (b) of the plug-in learner from above in that loss minimization in the equations above targets directly at Makarov bounds. However, the shortcoming (a) of the selection bias still persists. Furthermore, a new shortcoming (c) now emerges: The losses can be highly sensitive to badly estimated nuisance functions so that L PI (g, η) and L PI (g, η) differ significantly. Next, we develop an orthogonal learner that addresses all of the shortcomings. One-step bias correction. In order to address the before-mentioned shortcomings of the CA-learner, we employ the concept of (Neyman-)orthogonal losses [15, 35] . Informally, orthogonal losses are first-order insensitive to the misspecification of the nuisance functions, which introduces many favorable properties such as double robustness [127] . The CA-learner losses in Eq. ( 8) and ( 9) can be made orthogonal by performing a one-step bias correction [10, 63] . The one-step bias correction requires the knowledge of an efficient influence function, which has not yet been derived for Makarov bounds (→ Challenge 2 ). Hence, the following theorem presents one of our main theoretical results. Theorem 1 (Efficient influence function for Makarov bounds). Let P denotes P(Z) = P(X, A, Y ),and let y * Y (• | x) and u * [α,1] (• | x) be argmax/argmin sets of the convolutions (F 1 * F 0 ) Y (• | x) and (F -1 1 * F -1 0 ) [α,1] (• -0 | x), respectively. Then, under mild conditions on the conditional outcome distributions and for almost all values of δ ∈ ∆ and for all values of α ∈ (0, 1) (see Appendix D), average Makarov bounds are pathwise differentiable. Further, the corresponding efficient influence functions, ϕ, are as follows:EQUATIONEQUATIONEQUATION)where In the above theorem, we use red color to show the nuisance functions of P that are influencing the target estimand, i. e., averaged Makarov bounds. Therein, we also provide a Corollary 1, where we derive efficient influence functions for the target risks from Eq. ( 6) and ( 7), namely ϕ(L(g); P).I(X; η) = 1 (F 1 * F 0 ) Y (δ | X) > 0 ; y * isNote on infinite argmax/argmin sets. When argmax/argmin sets of the sup/inf-convolutions are infinite, average Makarov bounds are pathwise non-differentiable, and, thus, one-step bias correction is not possible as statistical inference becomes non-regular [47] . This result also holds for other causal quantities that contain sup/inf operators (e. g., for the policy value of the optimal treatment strategy [84, 106] ). Although, there exist approaches to perform inference in the non-regular setting [84] , we focus solely on the regular setting where pathwise differentiability holds (see Appendix D for a discussion on the generality of such a setting).Orthogonal leaner (AU-learner). Given the derived efficient influence function for the target risks, we perform a γ-scaled one-step bias correction [63] of the CA-learner losses, namely, L PI (g, η) + γP n ϕ(L(g); P) . The latter then yields our novel orthogonal AU-learner (see Corollary 2 in Appendix D). Our AU-learner effectively resolves all the above-mentioned shortcomings (see a comparison in Fig. 4 ). Formally, it aims at minimizing one of the following losses:EQUATIONEQUATIONF DR (δ, Z; η, γ) = F PI (δ | X; η) + γ C(δ, Z; η) and F DR -1 (α, Z; η, γ) = F PI -1 (α | X; η) + γ C -1 (α, Z; η),where C(δ, Z; η) and C -1 (α, Z; η) are given by Eq. ( 11) and ( 12), respectively; and γ ∈ (0, 1] is a scaling hyperparameter. We present a meta-algorithm of our AU-learner (with the CRPS target risk) based on cross-fitting in Algorithm 1 (AU-learner with the W 2 2 target risk follows analogously). Scaling hyperparameter. The scaling hyperparameter γ is introduced to tackle Challenge 3 from above, namely, that the pseudo-CDF term, F DR (δ, Z; η, γ), is not guaranteed to be a valid CDF for γ > 0 (both monotonicity wrt. δ and [0, 1]-constraint can be violated). 8 The same happens with the pseudo-quantiles of the AU-learner, F DR -1 (α, Z; η, γ), which could be non-monotonous wrt. α. The main intuition behind scaling is that it interpolates between the full AU-learner (γ = 1), that has favorable theoretical properties; and the CA-learner (γ = 0), for which the pseudo-CDFs and pseudoquantiles are valid CDFs and quantiles, respectively (we refer to Appendix E with visual examples). Hence, scaling mimics a learning rate of a Newton-Raphson method (usually considered as an analogy to the one-step bias correction [34] ). We found fixed values for the scaling hyperparameter γ to work well in all of our experiments and to improve the low-sample performance of our AU-learner.Algorithm 1 AU-learner (CRPS) via cross-fitting1: Input. Training dataset D = {xi, ai, yi} n i=1 , scaling γ ∈ (0, 1], folds K ≥ 2, δ-grid {δj ∈ ∆} n δ j=1 2: Output. Estimator of Makarov bounds g(δ, x) 3: for k ∈ {1, . . . , K} do ▷ First stage 4:Use {xi, ai, yi} k-1̸ =(i mod K) to fit η = (π, F0, F1) 5:for i : k -1 = (i mod K) do 6:Use η to infer the pseudo-CDF for the δ-grid:EQUATION)) 7:end for 8: end for 9: Fit the working model based on δ-grid: ▷ Second stage g = arg min g∈G L AU, CRPS (g, η)In the following, we formulate our second main theoretical result. For the results to hold, the nuisance functions η = (π, F 0 , F 1 )/η = (π, F -1 0 , F -1 1 ) need to be estimated independently from the second stage model g/g -1 . This could be done by either assuming a not-too-flexible class of models (such as a Donsker class of estimators and fitting all the models on the same dataset D) or by using a generic approach of cross-fitting [15, 35] .Theorem 2 (Neyman-orthogonality of AU-learner (informal)). Under the assumptions of the Theorem 1, the following holds for AU-learner from Algorithm 1 with the scaling hyperparameter γ = 1: 1. Neyman-orthogonality. Losses in Eq. (13) and Eq. ( 14) are first-order insensitive wrt. to the misspecification of the nuisance functions. 2. Rate double robustness. The bias from the misspecification of the nuisance functions is of second order and vanishes at the same rate as the fastest estimated nuisance functions.We refer to the Appendix D for the detailed formulation of the theorem and the proof. Notably, the rate double robustness has two important implications. (i) If at least one of the nuisance functions is estimated consistently, the target estimands are also estimated consistently. (ii) It allows to achieve a so-called quasi-oracle property [98] . This means that our AU-learner with (sufficiently fast) estimated nuisance functions performs nearly identical to the AU-learner with the ground-truth nuisance functions. 
neural instantiation with au-cnfs: 
 We now introduce a flexible fully-parametric instantiation of our AU-learner, which we call AU-CNFs. Therein, we employ conditional normalizing flows (CNFs) [105, 121] as main backbone for our AU-learner. CNFs are a flexible neural probabilistic method with tractable conditional densities, CDFs, and quantiles. Importantly, all three attributes of the CNFs (densities, CDFs, and quantiles) can be used for both training via back-propagation and inference [115] , which makes them a perfect model for both stages of our AU-learner.Architecture. The architecture of our AU-CNFs is inspired by interventional normalizing flows (INFs) [94] (a two-stage model for efficient estimation of potential outcomes densities). Our AU-CNFs consist of several CNFs corresponding to the two stages of learning, namely a nuisance CNF, which fits the nuisance functions, (π, F 0 , F 1 ) or, equivalently, (π, F -1 0 , F -1 1 ); and two target CNFs, which implement second stage working models for upper and lower bounds, G and G, respectively.Training & implementation. At the first stage of AU-CNFs learning, the nuisance CNF aims at maximizing the conditional log-likelihood and minimizing a binary cross-entropy via a joint loss. Then, we generate the pseudo-CDFs and pseudo-pseudo quantiles, as described in Algorithm 1. Therein, we set the δ/α-grid size to n δ = n α = 50 and discretize the Y-space/[0, 1]-interval to infer the argmax/argmin values, ŷ * /û * . Then, we proceed with the second stage of AU-CNFs learning, where we set γ = 0.25 for the CRPS loss and γ = 0.01 for the W 2 2 loss. We found the fixed values of γ to work well in all of the synthetic and semi-synthetic experiments (except for the IHDP100 dataset, where the overlap assumption is violated). We use the same training data for two stages of learning, as (regularized) CNFs as neural networks belong to the Donsker class of estimators [123] . We refer to Appendix F for more details on our AU-CNFs. 
experiments: 
 We now evaluate our AU-learner. For this, we use (semi-)synthetic benchmarks with the ground-truth conditional CDFs/quantiles of potential outcomes, i. e., F a (y | x)/F -1 a (u | x). In this way, we can infer the ground-truth Makarov bounds and use them for evaluation.Evaluation metric. We use evaluation metrics based on the target risks (as introduced in Sec 4.2). Specifically, we report root continuous ranked probability score (rCRPS) and Wasserstein-2 distance (W 2 ) based on training data (in-sample) and test data (out-sample).Baselines. We compare the proposed hierarchy of learners from Sec. 4 with CNFs as backbones. These are the plug-in learner (Plug-in CNF), IPTW-learner (IPTW-CNF), CA-learners (CA-CNFs (CRPS / W 2 2 )), and AU-learners (AU-CNFs (CRPS / W 2 2 )). The only relevant baseline found in the literature is a plug-in learner based on kernel density estimation [17, 77, 110] . 9 For this, we used distributional kernel mean embeddings [97] (Plug-in DKME), a standard conditional kernel density estimation method. Details on the baselines are in Appendix G. Synthetic data. We adapt the synthetic data generator (d x = 2) from [61, 93] by creating three settings with different conditional outcome distributions: normal, multi-modal, and exponential (see data generation details in Appendix H). In the synthetic data, the ground-truth Makarov bounds are less heterogeneous than the potential outcomes, and, hence, two-stage learners are expected to perform the best. We sample n train ∈ {100; 250; 500; 750; 1000} training and n test = 1000 test datapoints. The out-sample results are in Fig. 5 . Here, our AU-CNFs perform the best wrt. rCRPS in the majority of settings and different sizes of training data. We also report the results wrt. W 2 2 in Appendix I. HC-MNIST dataset. HC-MNIST is a high-dimensional semi-synthetic dataset (d x = 785), built on top of the MNIST dataset [52] (see details in Appendix H). Here, the heterogeneity of Makarov bounds is also smaller than that of potential outcomes, reflecting inductive biases in the real world. We report the out-sample performance of different methods in Table 2 (the Plug-in DKME is omitted due to a too-long runtime). Therein, our AU-CNFs (CRPS) achieve the best performance and, thus, scale well with the dataset size and the dimensionality of covariates. Further, our AU-CNFs (W 2 2 ) improve the performance of CA-CNFs (W 22 ). In general, we observe that the loss based on the CDF distance (i. e., CRPS) has a lower variance and is easier to fit. In Appendix I, we additionally report the results for another popular semi-synthetic benchmark, IHDP100 [46, 114] .Case study. In Appendix J, we provide a real-world case study based on the observational dataset from [7] . Therein, we demonstrate how our AU-learner (AU-CNFs) can be used to estimate the effectiveness of lockdowns during the COVID-19 pandemic. We estimate the probability of benefit from intervention (a special case of Makarov bounds with δ = 0). As expected, we observe a drop in the incidence rate is highly probable after the implementation of a strict lockdown. 
discussion: 
 Low-sample & asymptotic performance. In several experiments, especially in low-sample settings, the CA-learner or even the plug-in approach are performing nearly as well or even sometimes better than the AU-learner. This can be expected, as the best low-sample learner and the asymptotically best learner can, in general, be different [20] , and there is no single "one-fits-all" data-driven solution to choose the former one [19] . This can be explained by too small dataset sizes or the severe overlap violations (as is the case with the IHDP100 dataset; see Appendix I). Yet, only our doubly robust AU-learner offers asymptotic properties in the sense that it is asymptotically closest to the oracle (see Fig. 4 ). We thus argue for a pragmatic choice in practice (i. e., in the absence of ground-truth counterfactuals or additional RCT data) where our AU-learner should be the preferred method for the covariate-conditional Makarov bounds even in low-sample data.Future work. Our work sets a foundation for several extensions to estimate covariate-conditional Makarov bounds. For example, the estimation of the interval probabilities (see Appendix B) of the treatment effect can provide a connection with the existing works on total uncertainty with conformal prediction [1] . Additionally, one might want to study possible extensions of Makarov bounds tailored to high-dimensional outcomes.Limitations & broader impact. Our work is subject to the standard assumptions of the potential outcomes framework. We further make assumptions on the outcome distribution, though these are very mild. Nevertheless, we expect our work to have a positive impact, as it will help to improve the reliability of decision-making in medicine and other safety-critical fields.We are the first to offer a theory of orthogonal learning to quantify the aleatoric uncertainty of the treatment effect at the covariate-conditional level and present flexible neural instantiation.Partial identification and sensitivity models in the potential outcomes framework. Interventional quantities (such as CAPOs and distributions of potential outcomes) and some counterfactual quantities (e. g., CATE) are point identifiable in potential outcomes framework (see Fig. 6 ). However, relaxation of the unconfoundedness assumption renders those quantities partially identifiable or even nonidentifiable identifiable [89] . In this case, additional sensitivity models can be employed. Examples are the marginal sensitivity model [11, 23, 37, 38, 52, 53, 56, 119] and outcome sensitivity model [11, 104] . Other approaches suggest using instrumental variables [4, 45] or noisy proxy variables [42, 120] . Counterfactual quantities, on the other hand, are inherently non-identifiable (e. g., joint distribution of potential outcomes [26] , expected counterfactual outcomes of (un)treated [95] , and distribution of treatment effect [26] ). There is no general approach for deriving sharp bounds for partially identifiable causal quantities, therefore the above-mentioned works are not relevant in our setting. Distributional treatment effects and the distribution of the treatment effect. There are two seemingly similar notions in the existing literature: (a) distributional treatment effects [16, 29, 62, 92, 100] and (b) the distribution of the treatment effect [6, 26, 28, 33, 88] . However, (a) and (b) have two important differences:• (i) Interpretation. Distributional treatment effects represent the differences between different distributional aspects of the potential outcomes (e.g., Wasserstein distances, KL-divergence, or the quantile differences) [62] . Hence, they can answer questions like "How are 10% of the worst-possible outcomes with treatment different from the worst 10% of the outcomes without treatment?". Here, the two groups (treated and untreated) of the worst 10% contain, in general, different individuals. This is problematic in many applications like clinical decision support and drug approval. Here, the aim is not to compare individuals from treated vs. untreated groups (where the groups may differ due to various, unobserved reasons). Instead, the aim is to accurately quantify the treatment response for each individual and allow for quantification of the personalized uncertainty of the treatment effect. The latter is captured in the distribution of the treatment effect, which allows us to answer the question about the CDF/quantiles of the treatment effect. For example, we would aim to answer a question like "What are the worst 10% of values of the treatment effect?". Here, we focus on the treatment effect for every single individual. The latter is more complex because we reason about the difference of two potential outcomes simultaneously. Hence, in natural situations when the potential outcomes are non-deterministic, both (a) the distributional treatment effect and (b) the distribution of the treatment effect will lend to very different interpretations, especially in medical practice. In particular, the distribution of the treatment effect (which we study in our paper) is important in medicine, where it allows quantifying the amount of harm/benefit after the treatment [60] . This may warn doctors about situations where the averaged treatment effects are positive but where the probability of the negative treatment effect is still large. • (ii) Inference. The efficient inference of the distributional treatment effects only requires the estimation of the relevant distributional aspects of the conditional outcome distributions (e. g., quantiles) and the propensity score [62] . However, in our setting of the bounds on the CDF/quantiles of the treatment effect, we also need to perform sup/inf convolutions of the CDF/quantiles of the conditional outcomes distributions. Hence, while the definitions of (a) the distributional treatment effects and (b) the distribution of the treatment effect appear related, their estimation is very different.Alternatives to Makarov bounds. Potential alternatives to Makarov bounds vary depending on the type of aleatoric uncertainty. For example, explicit or implicit sharp bounds were proposed for:• The variance of the treatment effect, Var(Y [1] -Y [0]) [5] . Here, the sharp bounds are explicitly given by Fréchet-Hoeffding bounds [39, 49] and are, in general, different from Makarov bounds. • The interval probabilities, P(δ [33] . The sharp bounds on the interval probabilities are only defined implicitly and are also, in general, different from Makarov bounds (see Appendix B).1 ≤ Y [1] -Y [0] ≤ δ 2 )We are not aware of other bounds for measuring aleatoric uncertainty (e. g., kurtosis, skewness, or entropy). Importantly, the above-mentioned bounds on different measures of uncertainty are orthogonal to our work, and we focus on the bounds on the CDF/quantiles of the treatment effect.Efficient estimators and orthogonal learners. Efficient estimation of causal quantities (target estimands) was studied in the scope of (i) semi-parametric efficient estimation theory and, more general, (ii) orthogonal learning theory. Both theories (i) and (ii) rely on the concept of influence functions (pathwise derivatives) [63, 72, 122] , which allow performing a one-step bias correction for (i) a plug-in estimator of the target estimand or (ii) the population risk containing target estimand, respectively. (i) Semi-parametric efficient estimation theory [10, 73, 74, 107, 122] provides asymptotically efficient estimators for finite-dimensional target estimands (e. g., average treatment effect (ATE) and average potential outcomes (APOs)). On the other hand, (ii) orthogonal learning theory (or debiased ML) [15, 35, 96, 125] was developed for infinitely-dimensional (functional) estimands, like, CATE and CAPO. Orthogonal learning theory (or orthogonal learners) estimates target estimands by minimizing (Neyman-)orthogonal losses that are first-order insensitive to the misspecification of the nuisance functions. Specific examples of orthogonal learners for identifiable quantities include CAPO learners [125] and CATE learners [18, 64, 71, 96, 98] .Estimation of partially identifiable quantities. Orthogonal learners were also proposed for bounds on partially identifiable causal quantities. Examples include different interventional quantities such, e. g., in marginal sensitivity model [99] , in instrumental variables setting [79] , in noisy proxy variables setting [80, 118] . Examples from counterfactual quantities are, e. g., the variance of treatment effect and joint CDF of potential outcomes (Fréchet-Hoeffding bounds) [6] . More generally, estimation and learning for partially identifiable quantities was studied in econometrics as estimation of intersection bounds [14, 112, 113] . In this paper, we construct a novel orthogonal learner targeting at Makarov bounds on the CDF/quantiles of the CDTE.Total uncertainty quantification. Total uncertainty in potential outcomes framework can be generally quantified with two approaches: (i) Bayesian methods (e. g., Gaussian processes) and (ii) conformal prediction framework. (i) Bayesian methods [2, 3, 46] allow to infer posterior predictive distributions or credible intervals and for both potential outcomes and treatment effect, but under additional identifiability assumptions (e. g., an assumption of additive latent outcome noise, which renders treatment effect distribution identifiable). (ii) Conformal prediction aims at providing a valid predictive interval and was applied to quantify total uncertainty of predicting potential outcomes [76, 78, 90, 111] and treatment effect [1, 55, 58, 67, 78] . The latter works either make a similar additive latent outcome noise assumption or 'hide' non-identifiable treatment effect into the predictive interval (however, oracle predictive interval could never be reached in this case). We argue that total uncertainty needs to be split into epistemic and aleatoric explicitly to provide insights on the origin of uncertainty [50] , especially in our setting where aleatoric uncertainty of the treatment effect itself is partially identifiable (see Fig. 2 ).In the following, we provide an additional intuition on how Makarov bounds are inferred from conditional distributions of the potential outcomes, P(Y [a] | x). For example, Makarov bounds for the CDF of the CDTE are a composition of the sup/inf convolutions, applied to the conditional CDFs of the potential outcomes, and the rectifier functions (see Fig. 7 ).Makarov bounds can be inferred (i) analytically or (ii) numerically. (i) Analytic formulas were proposed for very simple distributions (e. g., a normal distribution [26] ). At the same time, the (ii) numerical approach is more flexible. For example, we can discretize the Y-space or [0, 1]-interval and perform maximization/minimization on that grid. Notably, in our experiments, we infer the ground-truth bounds with the approach (i), when the ground-truth conditional potential outcomes distributions are normal; and (ii) otherwise. We also plot the full Makarov bounds with a gray color.The sharpness of Makarov bounds proposed in [26] has one important characteristic. Specifically, although the Makarov bounds on the CDFs, F(δ | X), are CDFs themselves, they are not valid solutions to the partial identification task. This can be easily checked, i. e., the expectation wrt. to the Makarov bounds does not coincide with the point-identifiable CATE, τ(x) = E(∆ | x): - 0 -∞ F(δ | x) dδ+ ∞ 0 1-F(δ | x) dδ < τ (x) < - 0 -∞ F(δ | x) dδ+ ∞ 0 1-F(δ | x) dδ.(15) Therefore, there is an important distinction between so-called pointwise sharpness and a uniform sharpness [33] .In the case of uniform sharpness, a sharp bound coincides with the solution to the partial identification task. This implies that if a bound is uniformly sharp, then the joint bound on the set of quantities, evaluated in two (or more) points of ∆ or [0, 1], is also sharp. Many known bounds (e. g., Fréchet-Hoeffding bounds [39, 49] , and marginal sensitivity model bounds [23, 37, 38, 119] ) are uniformly sharp. Yet, Makarov bounds are only pointwise sharp.Recent works developed uniformly sharp bounds on the CDF of the CDTE [33] . However, their inference requires a special computational routine for every value of δ/α wrt. the CDF/quantiles of the CDTE. Their usage is further complicated by the fact that the CDFs of the uniformly sharp bounds correspond to mixed-type discrete/continuous random variables. We show an example of the uniformly sharp bounds on the CDF for δ = 0 in Fig. 8 . The uniformly sharp bounds may be useful for more complex aleatoric uncertainty quantities (e. g., interval quantities likeP(δ 1 ≤ Y [1] -Y [0] ≤ δ 2 | x))or simultaneous bounds on the variance and the CDF of the CDTE. Nevertheless, in many practical applications, pointwise sharp bounds (Makarov bounds) are enough and we focus on those in our paper. (Makarov) and uniformly sharp bounds on the CDF of the CDTE based on i = 7-th instance of the semi-synthetic IHDP100 dataset [46] . Pointwise (Makarov) bounds are shown in gray . Also, we show uniformly sharp bounds inferred for F(0EQUATIONF(δ | x7) F(δ | x7) τ (x7) = E(∆ | x7) F 0 (δ | x7)| x 7 ) = P(Y [1]-Y [0] ≤ 0 | x 7 ), namely F 0 (δ | x 7 ) (shown in red).We display the lower bound in (a) and the upper bounds in (b). Notably, the expectation wrt. F 0 (δ | x 7 ) coincides with the CATE, τ (x 7 ) = E(∆ | x 7 ).In the following, we provide a connection of the Makarov bounds for the continuous outcome Y (the main target of our paper), with the bounds on the fraction of the negatively affected (FNA) from [60] . The FNA is defined for the binary outcome, Y ∈ {0, 1}, asEQUATIONAs suggested by [60] , sharp bounds on the FNA are given byEQUATIONTo connect the Makarov bounds for continuous outcome and the bounds on the FNA, we formulate the following remark.Remark 1 (Makarov bounds for mixed-type outcome [129] ). Sharp bounds on the CDF of the treatment effect for the mixed-type outcome Y areEQUATIONP 0 = P(Y [0] = • | x) = P(Y = • | x, A = 0) is a probability mass function.Notably, for (absolutely) continuous outcome P 0 = 0, Eq. ( 18) matches Eq. (3) of Sec. 3. Remark 1 immediately yields the bounds on the FNA, proposed by [60] :EQUATIONEQUATIONEQUATIONEQUATIONEQUATIONEQUATION= E min µ 0 (X), µ 0 (X) + 1 -µ 1 (X)EQUATIONIn the following, we develop an improved single-stage learner, namely, an inverse propensity of treatment weighted (IPTW)-learner. First, we revisit the plug-in learner by defining its estimation objective. Then, we introduce the IPTW-learner, which addresses one of the shortcomings of the plug-in learner. At the end, we mention a surprising property of the IPTW-learner, namely, the orthogonality wrt. to target risk aiming at the potential outcome distributions.Empirical risk of the plug-in learner. We assume that the plug-in learner (Sec. 4.1) uses the estimators of the conditional outcome CDFs, F a , that aim at minimizing the following empirical risk:EQUATIONwhere l(•, •) > 0 is a probabilistic loss (e. g., negative log-likelihood, check score, or CRPS with an empirical CDF). Here, the plug-in learner has two possible variants, namely, S-and T-learner, depending on whether the conditional outcome distribution is learned by a single model or two models [20, 69] .IPTW-learner. The IPTW-learner addresses the selection bias of the plug-in learner. For this, it additionally employs an estimated propensity score, π. The estimated propensity score is used to re-weight the original probabilistic loss, l(•, •), in the following way:EQUATIONwhere πa (x) = a π(x) + (1a) (1π(x)). The IPTW-learner up-weights the loss of F 0 and F 1 in the treated and untreated populations, respectively.Orthogonality of the IPTW-learner. Interestingly, the IPTW-learner from Eq. ( 29) can be seen as an orthogonal learner targeting at the following risk:EQUATIONwhere g(•, x) ∈ G is a working model defined as in Sec. 4.2. This target risk aims to find the best approximation of the conditional potential outcome distribution, P(Y [a] | x), with the working model, g ∈ G.The orthogonality of the IPTW-learner wrt. the target risk aiming at the potential outcome distributions was formally proved in [125] . Therein, the authors notice that the target estimand (e. g., the CDF of one of the potential outcomes) coincides with one of the nuisance functions (i. e., F a ). Informally, the orthogonality follows from the fact that the working model, g, simultaneously fits the target estimand and the nuisance function.The orthogonality of the IPTW-learner can also be seen by ( 1) performing a one-step bias-correction of Eq. ( 28) and (2) setting the same estimator for the working model and the nuisance function, i. e., g = F a . For example, if the probabilistic loss is the CRPS with the empirical CDF, l(Y, g(•, X)) = Y (1{Y ≤ y}g(y, X)) 2 dy, the (1) one-step bias-corrected loss is given by [6, 43] LDR,a(g, η = (π, Fa)) =PnEQUATIONwhere πa (x) = a π(x) + (1a) (1π(x)). Then, after step (2), g = F a , we immediately yield the IPTW-learner from Eq. ( 29). Similarly, for the negative log-likelihood loss (NLL) l(Y, g(•, X)) = log g(Y, X), the (1) one-step bias-correction is given by [65, 94] LDR,a(g, η = (π, P)) =PnEQUATIONwhere πa (x) = a π(x) + (1a) (1π(x)), and P is an estimator of the conditional density of the outcome. Again, after (2) setting g = P, it is easy to see that the minimization of one-step bias corrected loss in Eq. ( 32) is equivalent to the minimization of the IPTW-learner's objective in Eq. ( 29) (where P is used in place of F a ). This is possible due to two facts: (1) both entropy terms in Eq. ( 32), -Y log g(y, X) g(y, X) dy, only require the minimization wrt. to the working model g under the logarithm; and (2) these entropies are minimal as the cross-entropy for any distribution is minimal when evaluated with itself.In the following, we provide the main theoretical results of our paper. We use the following additional notation: δ{•} is a Dirac delta function, a ≲ b means there exists C such that a ≤ C • b. Also, in the following theorems, we use red color to show the nuisance functions of P that are influencing the target estimand.We start with deriving the efficient influence functions for the average Makarov bounds and, afterwards, for the target risks. For that, we make two mild assumptions: (1) one the conditional outcome distributions and (2) another on the set of δ where linear rectifiers are differentiable. These assumptions allow us to (1) handle sup-/inf-convolutions as max-/min-convolutions with a finite number of argmax/argmin values and to (2) get a derivative of the linear rectifiers.Assumption 1 (Finite argument sets). We assume that the outcome space Y is compact. Also, we assume that conditional outcome CDFs, F a (y | x), are continuously differentiable and consist of a finite number of strictly concave / convex regions.Assumption 1 implies that sup-/inf-convolutions are achieved by some finite set of values in Y. Furthermore, Assumption 1 is a special case of the margin assumption (Assumption 3.2) from [112] .Yet, we find our version to be more interpretable and many regular distributions satisfy it, e. g., an exponential family, finite mixtures of normal distributions, etc.Assumption 2 (Differentiability of linear rectifiers). Values of δ ∈ ∆ are considered to satisfy differentiability of linear rectifier for somex ∈ X if (F 1 * F 0 ) Y (δ | x) ̸ = 0.It is easy to see that, if Assumption 1 holds, the new Assumption 2 will hold for almost all δ ∈ ∆ and x ∈ X . Notably, no additional requirements are needed for the values of α ∈ [0, 1]. Therefore, the formulation of the following theorem is not too restrictive. (F 1 * F 0 ) Y (• | x) and (F -1 1 * F -1 0 ) [α,1] (• -0 | x), respectively. Then, under the mild assumption of the finite argument sets (Assumption 1), average Makarov bounds are pathwise differentiable for values of δ ∈ ∆ that satisfy the differentiability of linear rectifiers (Assumption 2) and for all values of α ∈ (0, 1). Further, the corresponding efficient influence functions, ϕ, are as follows:EQUATIONEQUATIONEQUATIONEQUATION)where y * is some value from y * Y (δ | X); u * is some value from u * Proof. We start the proof by employing the properties of the efficient influence functions, namely, product rule and chain rule; and some existing efficient influence functions (e. g., for conditional expectation [63] ).The efficient influence function of the lower averaged Makarov bound is as follows:EQUATIONFurthermore, the inner term isEQUATION)Then, the efficient influence function of the sup-convolution can be derived under the finite argument sets assumption and using the envelope theorem (we refer to Appendix C.2.6 and C.2.7 of [83] for further details):EQUATIONwhere y * is some value from y * Y (δ | x). The latter can be expanded by using product and chain rules and the efficient influence function for the conditional expectation. Specifically, we do the following: EQUATIONEQUATIONdy.The terms (2) and (4) yield well-known efficient influence function for CDFs, i. e.,EQUATION)Let us now consider the remaining terms (1) and (3):EQUATIONwhere the last equality holds due to the properties of the argmax, y * ∈ y * Y (δ | x). Namely, under the necessary condition for a local maximum, we haved dy (F 1 (y | x) -F 0 (y -δ | x)) = P(Y = y | x, A = 1) -P(Y = y -δ | x, A = 0) = 0.In the context of the efficient influence functions, this means that the Makarov bounds are first-order insensitive to the misspecification of argmax/argmin. Interestingly, a similar result was demonstrated for the efficient influence functions of the policy values of the optimal policies [82, 85] .Finally, we obtain an efficient influence function for the lower Makarov bound by expanding the part ( * ) of Eq. ( 40):EQUATIONwhere C(δ, Z; η) is defined in Eq. ( 35).The Makarov bounds on the quantiles are derived analogously. However, several last steps are different. The differences start from Eq. ( 44):EQUATIONwhere u * is some value from u * [α,1] (α | x). Now, again, the Makarov bounds for quantiles are first-order insensitive wrt. argmax/argmin values misspecification. Thus, we can employ the efficient influence function for the quantiles [31, 32] :EQUATIONEQUATIONEQUATION)The latter then yields the final part ( * ) of Eq. ( 40) for the efficient influence function:EQUATIONwhere C -1 (α, Z; η) is defined in Eq. ( 36).Corollary 1 (Efficient influence functions of the target risks). Let the Assumption 1 of the Theorem 1 hold. Then, the efficient influence functions, ϕ, of the target risks in Eq. ( 6) and Eq. ( 7) are as follows:ϕ(L CRPS (g); P) = ∆ 2 F(δ | X; η) -g(δ, X) C(δ, Z; η) + F(δ | X; η) -g(δ, X) 2 dδ (62) -E ∆ F(δ | X; η) -g(δ, X) 2 dδ , ϕ(L W 2 2 (g -1 ); P) = 1 0 2 F -1 (α | X; η) -g -1 (α, X) C -1 (α, Z; η) + F -1 (α | X; η) -g -1 (α, X) 2 dα (63) -E 1 0 F -1 (α | X; η) -g -1 (α, X) 2 dα ,where C(δ, Z; η) and C -1 (α, Z; η) are defined in Eq. ( 35) and (36), respectively.Proof. We start by using the properties of the efficient influence function, namely, chain and product rules. Then, the efficient influence function for the CRPS risk of the lower bound is as followsEQUATION)We then note, that the inner term of the ( * ) can be expanded as:EQUATION)Finally, the term ( * ) of Eq. ( 65) equals toEQUATIONEQUATIONThe derivations for the upper bound and for the W 2 2 target risk is fully analogous. Corollary 2 (One-step bias-corrected estimator of target risks). Let the Assumption 1 of the Theorem 1 hold. Then, the γ-scaled one-step bias-corrected estimator of the target risks Eq. (6) and Eq. ( 7) is given by the following:EQUATIONEQUATIONF DR (δ, Z; η, γ) = F PI (δ | X; η) + γ C(δ, Z; η) and F DR -1 (α, Z; η, γ) = F PI -1 (α | X; η) + γ C -1 (α, Z; η),where C(δ, Z; η) and C -1 (α, Z; η) are given by Eq. ( 35) and (36), respectively; and γ ∈ (0, 1] is a scaling hyperparameter.Proof. The γ-scaled one-step bias-corrected estimator of the CRPS target risk proceeds as follows:EQUATIONEQUATIONEQUATIONEQUATIONThe minimization of the latter wrt. g is then equivalent to the minimization ofEQUATIONThe proof for the γ-scaled one-step bias-corrected estimator of the W 2 2 target risk is analogous.Now, we proceed with the second main theoretical result. Here, we use additional notation. Definition 1 (Neyman-orthogonality [35, 96] ). A risk L is called Neyman-orthogonal if its pathwise cross-derivative equals to zero, namely,EQUATIONwhere [35] , g * = arg min g∈G L(g, η), and η is the ground-truth nuisance function.D f F (f )[h] = d dt F (f + th)| t=0 and D k f F (f )[h 1 , . . . , h k ] = ∂ k ∂t1...∂t k F (f + t 1 h 1 + • • • + t k h k )| t1=•••=t k =0 are pathwise derivativesInformally, this definition means that the risk is first-order insensitive wrt. to the misspecification of the nuisance functions. Notably, the pathwise derivative in the direction of the Dirac delta distribution coincides with to the efficient influence function [34, 63] , i. e., D P F (P)[δ{Z -•} -P(Z = •)] = ϕ(F (P); P), where P(Z = •) is the PDF of the P(Z). Theorem 2 (Neyman-orthogonality of AU-learner). Under the assumptions of the Theorem 1, the following holds for AU-learner from Algorithm 1 with the scaling hyperparameter γ = 1:1. Neyman-orthogonality. Population versions of the empirical risks in Eq. (71) and Eq. ( 71) are first-order insensitive wrt. to the misspecification of the nuisance functions, i. e.,EQUATIONEQUATIONwhere g * = arg min g∈G L AU,CRPS (g, η) and g -1 * = arg min g -1 :g∈G L AU,W 2 2 (g -1 , η). 2. Rate double robustness. The bias from the misspecification of the nuisance functions is of second order and vanishes at the same rate as the fastest estimated nuisance functions. Specifically, with the probability of at least 1δ, the following holdsEQUATIONwhere Rate D (G, δ/2; g, η) is the bound for the estimation error of the second stage of learning, namely,EQUATIONFurthermore, if the nuisance functions are estimated sufficiently fast, i. e., ∥η -η∥ 2 = o(n -1/2 ), then the AU-learner achieves the quasi-oracle property (the estimation error of the second stage with the estimated nuisance functions behaves in the same way as if the ground-truth nuisance functions were used).Proof. 1. Neyman-orthogonality. The Neyman-orthogonality follows by the construction of the AU-learner as a one-step bias-corrected estimator. Specifically, it is easy to verify that the pathwise cross-derivative from Eq. ( 78) is equal to zero.Let us consider the CRPS risk. First, we find a pathwise derivative wrt. working model g:EQUATIONEQUATIONEQUATIONThen, we find derivatives wrt. to different nuisance functions (e.g., the propensity score):EQUATIONEQUATIONEQUATIONEQUATIONEQUATIONEQUATIONEQUATIONwhereEQUATIONEQUATIONEQUATIONEQUATIONEQUATIONEQUATIONwhereI(X; η) = 1 (F 1 * F 0 ) Y (δ | X) > 0 or 1 (F 1 * F 0 ) Y (δ | X) < 0 ; ỹ * (t) is some value from ỹ * Y (δ | X); ỹ * Y (• | X) are the argmax/argmin sets of the convolutions (F 1 + t( F1 -F 1 ) * F 0 ) Y (• | X); and ( * ) = 0 follows from the the same considerations as in Eq. (51) . Analogously, the pathwise derivative wrt. F 0 can be shown to be equal to zero. We refer to the appendices of [96] for more details.The Neyman-orthogonality of the W 2 2 population risk can be proved in a similar fashion. First, the pathwise derivative wrt. g -1 has a similar form to Eq. ( 85), namelyEQUATIONFurthermore, similarly to the CRPS target risk, the cross-derivative wrt. to the propensity score isEQUATIONEQUATIONEQUATIONwhere u * is some value from u * [α,1] (α | X). The cross-derivative for the upper bound follows similarly.Finally, to show that a cross-derivative wrt. F -1 1 , we make use of the following property of the derivative of the quantiles (inverse function rule):EQUATIONwhere P( Ỹ = • ; t) is the density function for a distribution with quantiles F -1 (α) + t( F-1 (α) -F -1 (α)). Then, the pathwise derivative isEQUATIONTherefore, the cross-derivative isEQUATIONEQUATION-1EQUATIONEQUATION-EQUATIONwhere ũ * (t) is some value from ũ *[α,1] (α | X); and ũ * [α,1] (α | X) is the argmin set of the convolution (F -1 1 + t( F-1 1 -F -1 1 ) * F -1 0 ) [α,1] (• -0 | X).The cross-derivatives for the upper bound wrt. F -11and for both upper and lower bounds wrt. F -1 0 follow similarly. 2. Rate double robustness. The result is a direct application of Theorem 1 in [35] . First, it is easy to see that the Assumptions 1-4 from [35] hold for the population versions of the empirical risks of our AU-learner (i.e., CRPS and W 22 ).Fig. 9 provides a visual comparison of different pseudo-CDFs for our AU-learner, with and without scaling γ. Therein, we see that the scaling hyperparameter helps to enforce the constraints on the pseudo-CDFs (i. e., [0, 1]-boundedness and monotonicity). Our AU-CNFs allow us to implement the Algorithm 1 of our AU-learner (see Fig. 10 ) by combining several conditional normalizing flows (CNFs) [105, 121] . It consists of a (i) nuisance CNF and (ii) two target CNFs (upper and lower). ( 1) The nuisance CNF aims to fit the nuisance functions, (π, F 0 , F 1 ) or, equivalently, (π, F -1 0 , F -1 1 ). (2) Upper and lower target CNFs constitute the second stage working models, namely, G and G, and minimize the loss of our AU-learner.(1) Nuisance CNF. The nuisance CNF has three components, similarly to [94] . These are two fully-connected subnetworks (FC 1 and FC 2 ) and a CNF, parametrized by θ. The two subnetworks FC 1 and FC 2 form a hypernetwork, which outputs the conditional parameters, θ = θ(X, A). This allows us to flexibly model the conditional outcome distribution.The nuisance CNF has the following joint loss for the nuisance functions: L N = L NLL + αL π . Here, L NLL is a conditional negative log-likelihood loss, L π is a binary cross-entropy, and α > 0 is a hyperparameter. We additionally employed noise regularization to regularize the conditional negative log-likelihood loss [108] .(2) Upper and lower target CNFs. The upper and lower target CNFs use the pseudo-CDFs / pseudo-quantiles, generated by the nuisance CNF, and then implement a second stage loss of our AU-learner. Both target CNFs have the same structure. Specifically, they have a fully-connected subnetwork, FC 3 , and a CNF, parametrized by β. Analogously, FC 3 serves as a hypernetwork so that the parameters can be conditioned on X: β = β(X).To fit the target CNFs, we use a second stage loss of our AU-learner, namely, Eq. ( 13) or Eq. ( 14). For that, we discretize the Y-space or the [0, 1]-interval of u into n d values and infer argmin/argmax values based on those grids. Then, to approximate the integrals, we do the same for the ∆-space and the [0, 1]-interval of α. The later creates a δ/α-grid with n δ /n α points. Those grids are later used for a rectangle quadrature integration. Furthermore, we also regularize the target CNFs by applying the noise regularization [108] .Implementation. We implemented our AU-CNFs using PyTorch and Pyro. For the CNFs of both stages of learning, we used neural spline flows [24] with a standard normal distribution as a base distribution. Neural spline flows build an invertible transformation based on invertible rationalquadratic splines and, thus, allow the direct inference of the (conditional) log-probability, CDF, and quantiles. Neural spline flows are characterized by two main hyperparameters, namely, a number of knots n knots and a span of the transformation interval, [-B, B]. The number of knots, n knots , controls the expressiveness of the flow. The span B defines the support of the transformation. In our experiments, we tune the number of knots n knots and set the span B via a heuristic depending on sample max/min values (as Y is assumed to be compact).Training. To train our AU-CNFs, we make use of Algorithm 1. However, both first and second stage models are fit on the same training data D without cross-fitting as (regularized) CNFs as neural networks belong to the Donsker class of estimators [123] . Training of our AU-CNFs proceeds as follows: (1) we fit the nuisance CNF; (2) we freeze the nuisance CNF and generate pseudo-CDFs/pseudo-quantiles; and (3) we train the upper and lower target CNFs. The hyperparameters are then as follows:1. First stage. We used stochastic gradient descent (SGD) with a minibatch size b N , n e,N = 200 epochs and a learning rate η N . Furthermore, we set the loss coefficient to α = 1. Both the number of hidden units of FC 1 /FC 2 and the size of the output of the FC 1 are set to 10. The nuisance CNF has the number of knots n knots, N and the span B = max i (ỹ i )min i (ỹ i ) + 5, where ỹi are standard normalized outcomes y i . The intensities of the noise regularization for the input and the output are set to σ 2 x and σ 2 y , respectively. 2. Intermediate stage. We set n d = 200 and n δ = n α = 50. Furthermore, we set the scaling hyperparameters γ = 0.25 for the CRPS loss and γ = 0.01 for the W 2 2 loss. We clipped too low propensity scores (lower than 0.05). x . The number of hidden units of FC 3 is also set to 10. The target CNFs have the number of knots twice larger than the nuisance flow, n knots, T = 2 n knots, N , and the span B = max i (a i ỹi )min i (a i ỹi ) + max i ((1a i ) ỹi )min i ((1a i ) ỹi ) + 5, where ỹi are standard normalized outcomes y i . To further stabilize the training of the target CNFs, we employed an exponential moving average (EMA) of the target CNFs parameters [103] with a smoothing hyperparameter λ = 0.995.We demonstrate the detailed training procedure of our AU-CNFs (CRPS) in Algorithm 2 (AU-CNFs (W 22 ) follow analogously). Hyperparameter tuning. We performed extensive hyperparameter tuning only for the nuisance CNF. The following hyperparameters are subjects to tuning: the minibatch size n b,N , the learning rate η N , the number of knots n knots, N , and the intensities of the noise regularization, σ 2x and σ 2 y . Further details of hyperparameter tuning are provided in Appendix G. The hyperparameters of the target CNFs for all the experiments are either kept fixed or are inherited from the nuisance CNF. (i) 1 (X) 6: Noise regularization: ξx ∼ N (0, σ 2 x ), ξy ∼ N (0, σ 2 y ), R, Ỹ ← R + ξx, Y + ξy H.3 IHDP100 datasetThe Infant Health and Development Program (IHDP100) [46, 114] is a standard semi-synthetic benchmark for treatment effect estimation. It contains 100 train/test splits with n train = 672, n test = 75, and d x = 25. Yet, this dataset contains severe overlap violations, which makes the methods using propensity re-weighting unstable [20, 21] .The IHDP100 dataset samples synthetic outcomes from the conditional normal distribution, N (µ i , 1), where µ i are CAPOs provided in the dataset. Therefore, the ground-truth Makarov bounds are given by the half-normal distributions [26] . Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: Our work does not pose such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: All the datasets, used in our work, are cited accordingly.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.In the Neyman-Rubin potential outcomes framework[109], a potential outcome Y [a] refers to the value an outcome variable would take for an individual under a specific treatment or intervention a. Each individual has multiple potential outcomes -one for each possible treatment condition -but only one of these outcomes is observed.Notably, (a) the distributional treatment effects and (b) the distribution of the treatment effect (our setting) are both different interpretationally and inferentially. That is, (a) provide contrasts between distributions of potential outcomes and are point identifiable, and (b) work on the distribution of the difference of potential outcomes and are only partially identifiable. See Appendix A for further details.Code is available at https://github.com/Valentyn1997/AU-CNFs.Epistemic uncertainty[50] relates to the uncertainty of fitting a model on finite data and reduces to zero as data size grows. In contrast, aleatoric uncertainty originates from the inherent randomness of the outcome and is irreducible wrt. data size.The work in[54] assumes the possibility of finding feasible Kantorovich dual functions to a target functional (e. g., CDF of the treatment effect). By doing so, the authors are able to infer valid partial identification bounds on very general functionals; yet, the sharpness can not be practically guaranteed.Also known as an individual treatment effect. The term 'individual' should not be confused with the term 'covariate-conditional', which refers to the causal quantity and not the random variable itself.By postulating a restricted working model class, G, we might compromise on the sharpness and end up having looser bounds. Yet, this looseness bias is only relevant in the infinite data regime. In the finite-sample regime, the feasibility of the low-error estimation is a much more important problem.The issue of pseudo-outcomes violating constraints was also raised wrt. DR-learner for CATE[71,125] when the outcome space is bounded (e. g., Y = [0, 1]) and for the CDFs of potential outcomes[43]. Yet, the optimal solution remains an open research question.The work of[54] focuses on averaged (population-level) Makarov bounds for fixed values of δ/α and it is unclear how to extend it to our setting.We omitted the multi-modal setting, as the mixture distribution does not provide direct ground-truth quantiles and, thus, Makarov bounds on the quantiles.The data is available at https://github.com/nbanho/npi_effectiveness_first_wave/blob/ master/data/data_preprocessed.csv. 
 
 Appendix:  


Decision: The paper introduces a method called AU-learner to estimate the conditional distribution of treatment effects (CDTE) for binary treatments and continuous outcomes. The authors use Makarov bounds for partial identification and conditional normalizing flows for estimation, developing a doubly robust learner that satisfies Neyman-orthogonality. 

Pros:

+ Addresses an important problem in causal inference by quantifying aleatoric uncertainty in treatment effects

+ Develops a theoretically grounded method with useful properties (double robustness, Neyman-orthogonality)

+ Provides a comprehensive set of experiments and supplementary materials

Cons:

+ The presentation is dense with technical details, making it challenging to read and understand

Reviews: review: paper_summary: paper_summary: The authors study the distribution over the individual treatment effect for binary treatments, continuous outcomes, and observed, potentially high-dimensional confounders. The authors build up on prior work on Makarov bounds, to develop a new method to lower/upper-bound the conditional CDF and the quantile function of the treatment effect. They develop a doubly robust learner and perform a theoretical and empirical evaluation. In a series of experiments, They compare several instances of their method with a baseline based on kernel density estimation. While the bound-based methods clearly outperform the baseline, they yield rather similar results.
strength_weakness: strength_weakness: The plugin estimator based on a conditional normalizing flow seems to perform rather well in the experiments. The CA learner w/o bias correction seems to perform as good as the AU learner (slightly worse in Table 2, slightly better in Table 4 in the appendix). The added complexity of the AU learner, in particular, the one-step bias correction which comes with an additional tuning parameter, renders the practical value questionable. It seems that a lot of work went into the AU learner; but it might still be worth focusing the presentation more on the simpler CA learner (e.g., discussing the role of g in the main doc).
questions: questions: I do not fully understand the role of the working model g \in G in the CA learner. The Makarov bound is a CDF and G may include all possible CDFs; so if G is rich enough, the CA learner may return the same CDF estimate as the plugin estimator. If we restrict G, why does this cause tighter bounds?
limitations: limitations: NeurIPS Paper Checklist is provided; no concerns.
paper_summary: paper_summary: The authors introduce AU-learner, a method to estimate the conditional distribution of treatment effects (CDTE) and hence capture the  variability in the treatment effect. They use Makarov bounds for partial identification and use conditional normalizing flows for estimation. Further, they show that AU-learner satisfies Neyman-orthogonality and double robustness.
strength_weakness: strength_weakness: (W1) I would have liked a more comprehensive discussion of how the proposed bounds compare against existing ones in terms of sharpness. For instance, it is not clear to me if the bounds you get are sharper than [58] in the binary outcomes setting. 

(W2) I think it would be useful to compare the proposed methods against the non-doubly robust estimators. In practical settings both outcome and propensity models are likely misspecified, hence the main advantage of doubly robust estimators are the faster theoretical rates. In this setting (with both nuisances misspecified), it is not clear if the doubly robust version of the bounds is always better in practice and some more experiments would be nice to explore this. 

(Minor) In the extended related work on partial identification and sensitivity models (line 658) you could also mention more recent approaches that incorporate RCT data together with proxy and instrumental variables, e.g. the works of [1] and [2].


[58] Nathan Kallus. “What’s the harm? Sharp bounds on the fraction negatively affected by treatment”. In: Advances in Neural Information Processing Systems. 2022.

[1] Falsification of Internal and External Validity in Observational Studies via Conditional Moment Restrictions. Hussain et al. AISTATS 2023.


[2] Hidden yet quantifiable: A lower bound for confounding strength using randomized trials. De Bartolomeis et al. AISTATS 2024.
questions: questions: (Q1) For a binary outcome, could you comment on how your bounds relate to Kallus [58]? Why don't you consider them as a baseline (say for estimating bounds on the quantiles of the treatment effect distribution)?
limitations: limitations: Yes
paper_summary: paper_summary: In this paper, the authors propose a method to quantify the aleatoric uncertainty of the treatment effect.
For this, authors estimated Makarov bounds on the CDF and quantiles of the CDTE, and then showed, how one can build a learner, which has properties of Neyman-orthogonality and double robustness.
The authors proved the abovementioned theoretical properties of the resulting learner and demonstrated the usefulness of the proposed approach in a series of experiments on synthetic and real-world data.
strength_weakness: strength_weakness: I don't have critical concerns about the paper, but I believe the following can improve the paper:
- I find the paper logically well-structured, but at the same time challenging to read, as it is too concentrated with technical details. I would suggest authors reconsider the narrative, concentrating more on the conceptual ideas in the main part, and moving technical details to the Supplementary;
- The field of treatment effect estimation is quite special and narrow (from my point of view) among the machine learning community, and it is worth adding some clarifications to the terms and notation used. For example, what is $Y$? It is just said that it is a continuous outcome, without any intuitions of what it could be.
- I feel that the paper missing the discussion of the alternatives of the proposed approach. For example, why specifically Makarov bounds were chosen, but not other possible alternatives? 
- In Table 2, it is worth explicitly writing that CNF corresponds to the Plug-in learner (as it was called throughout the paper). I would suggest Plug-in-CNF (like written for IPTW-CNF).
- In Line 268 there is a minor typo: Should be Eq.13 and Eq. 14.
questions: questions: - What are the alternatives for Makarov's bounds?
- Are such options like conformal predictions or plain confidence intervals somehow useful to estimate AU in the context of treatment effect? 
- Compared to the vanilla Plug-in estimator (say Plug-in-CNF), what is the additional computational overhead of AU-learner?
- In lines 212-216, it is mentioned that CA-learned still has shortcomings (a) (and a new shortcoming (c)). In practice, how severe is the shortcoming (a) (selection bias) compared to AU-learner? Can you provide a toy example?
limitations: limitations: One of the limitations, not mentioned by authors, is that the approach requires (conditional) Normalizing Flows, and hence might not work well in the high dimensional scenario.
paper_summary: paper_summary: The authors propose a partial identification of quantiles of the individual treatment effect, which are not point-identifiable in general. The authors justifiably argue that characterizing the distribution of individual treatment effects gives a better idea of the aleatoric uncertainty in a causal-inference problem. The proposed method is doubly robust, works with heterogeneous treatment effects, and requires minimal additional assumptions about the data.
strength_weakness: strength_weakness: * The benefit and/or novelty of the CA-learner is questionable (also reflected in the results Table 2) and I wonder if its exposition is taking up valuable space. It is an insightful point on lines 182--185 that learning Makarov bounds could benefit from inductive biases of lower heterogeneity than the conditional CDFs. However, it is unclear if the CA-learner loss really incorporates that inductive bias and if so, how much.

 * The conclusion is a bit grand. Arguably, previous papers like those highlighted in Table 1 have proposed robust methods for quantifying a version of aleatoric uncertainty of causal effects.
questions: questions: My main question has to do with recent related work. In particular, Ji et al. [52] appear to solve a similar problem, and the quick dismissal of that approach in this paper because they "made special optimization assumptions" needs further discussion. It would be helpful to spell out what these assumptions are in a concrete sense. In doing so, the authors could also discuss whether these two approaches have any fundamental commonalities, or if the partial identifications are expected to be materially different. Finally, it would be nice to see [52] appear as a baseline in the empirical evaluations, although I understand this could be difficult to implement, especially since [52] is relatively recent.
limitations: limitations: N/A
paper_summary: paper_summary: The paper provides the orthogonal estimator for the distributional treatment effect (the conditional CDF of $Y[1]-Y[0]$).
strength_weakness: strength_weakness: Despite the strong technical details, the paper is poorly written in overall. 

__1. Weak motivation__

A current shape of Introduction is weakly motivated. Firstly, the way that the paper motivates the problem is misleading:
> Methods for quantifying the aleatoric uncertainty of the treatment effect have gained surprisingly little attention in the causal machine learning community.

This sentence is incorrect, since there are literatures on quantile regression and semiparametric density estimation, as reviewed in Section 2. 

More importantly, the introduction doesn't provide the motivation of the problem against the following question: _why a community need the proposed estimator, given that the quantile estimator can capture the distributional treatment effect_.  


__2. Difficult to understand due to insufficient information__

Another issue that the paper has (especially in Introduction) is its lack of back-ground information that readers may need to comprehend. Specifically, in Introduction, the paper doesn’t provide any definition or clue what CDTE is. I understand that the CDTE is defined in the caption of Figure 1 as $P(Y[1] - Y[0] \leq \delta \mid x)$. However, $Y[a]$ is undefined, and this key quantity should be in the main body of the text. Also, even if Figure 1 aims to provide a whole summary of the paper, authors at Introduction have insufficient knowledge to comprehend it. In other words, Figure 1 is too detailed to be presented in Introduction section. Since Figure 1 can only be understood by those who entirely digested the paper from the beginning to the end, Introduction section is not the right position where the Figure 1 is located. I understand the goal of Figure 1, but it doesn't achieve the goal because of insufficient background information. The same issue happens to Figure 2 and Table 1. For example, in Table 1, technical terms like AIPTW, hold-out residual, optimization assumptions are undefined, so it's hard to appreciate the contribution of this paper. 

__3. Fuzzy description on contribution__ 

First, the terms like "Aleatoric uncertainty" and "distributional treatment effect" are used for denoting the same target quantity. Given that the "distributional treatment effect" is clearly describing the problem, I don't see why the authors want to use "Aleatoric uncertainty" as a title and employ these two words to denote the same estimand. 

Second, the contribution is wrongly described. Consider this sentence: 
> AU-learner solves all of the above-mentioned challenges 1 – 3 .

The AU-learner doesn't address Challenge 1 and 2. Challenge 1 means that the distributional treatment effect is not identifiable. The Makarov's bound, _NOT_ AU-learner, is employed to address Challenge 1. Challenge 2 is actually the same as Challenge 1, since it means that there are no known nuisances-based functional for the distributional treatment effect. Again, the Makarov's bound is used to address Challenge 2. AU-learners are representing the approximated quantity of the target estimand in terms of nuisance functionals. 

Finally, third contribution "flexible deep learning instantiation of our AU-learner" is scarcely described only in Section 5. The description needs to be much improved. 

__4. Little focus on the real contribution__ 

The real contribution of this paper, compared to the existing works in Table 1, is to provide the doubly robust conditional distributional treatment effect for the bounds of the conditional CDF of treatment effects. However, little focus and efforts have been made for this contribution. For example, if developing a doubly robust estimator is a contribution, then corresponding results such as detailed error analysis, a closed form of estimators, a detailed recipe of the proposed estimator for the specified working model, how to minimize the losses in Equations (8,9), a simple example, and assumptions should be described.
questions: questions: 1. Are rate-doubly-robustness and Neyman orthogonality violated when the scaling hyperpameters are not $1$? 

2. Is the Makarov bound sharp? 

3. What are the practical examples of the distributional treatment effect? 

4. In line 178, is this phenomenon officially termed the selection bias?
limitations: limitations: The paper is limited to the setting where the ignitability holds.

Comments: title: Response
comment: [3/3] 

> "However, the categorization seems to come out of nowhere and requires justification."

Thank you for the question. Our justification for the above categorization of causal quantities is based on the following rationale.
1. **The AU of the potential outcomes:** The distribution of the treatment is non-identifiable, and, hence, we wanted to first distinguish our work from causal quantities that are identifiable. The reason is that the latter can be addressed by point identification, while the former (our problem) must be addressed by partial identification or making stronger assumptions. 
2. **The distributional treatment effects:** Here, we want to distinguish our work from causal quantities that are contrasts between AUs of both the potential outcomes. Thereby, we aim to spell out clearly that the distributional treatment effects and our distribution of the treatment effect are both very different interpretationally and inferentially. 
3. **The distribution of the treatment effect**: Here, we aim to survey works related to our setting, namely, the distribution of the treatment effect. 

Importantly, the above categorization is neither universal nor final but we used it as an informal guidance to structure the related work in our paper. If there is a better way to categorize the causal quantities, we would be happy to incorporate it into our paper. 

**Action**: We will spell out the rationale for the categorization in our Related Work section more clearly. 


Again, we are sorry for misinterpreting your initial question and for the confusion this has caused. We hope that we have addressed all of your concerns and assure you that they can be easily fixed in our revised manuscript. Should you have any further questions, please let us know so -- we would do our best to answer them promptly. Thanks again for reviewing our submission.
[2/3]


### Regarding the categorization of causal quantities


Below, we respond to your original question about the differences in causal quantities and what appears to be your main concern. We are convinced that the problem is easy to fix for the final version of the manuscript.

> “However, as stated in the authors' response, this categorization is not mentioned in the paper, and no clues are provided regarding this categorization.”


We apologize that we mentioned the categorization into the three streams of literature only as plain text, while, after reading your comment, we realized that we should have made it more explicit (e.g., by adding a formal categorization via a table). In general, three different streams are relevant to our work:  
- the AU of potential outcomes (lines 95-97), 
- the distributional treatment effects (lines 98-99), and 
- the AU of the treatment effect (lines 101-122). 

In the submitted version of the paper, we only made a clear cut between the identifiable causal quantities (1.)+(2.) vs. non-identifiable (3.), which was the main distinction we aimed to communicate due to reasons of space. Below, we appreciate the opportunity to explain the rationale behind our categorization. 

**Action:** We will revise our paper and make the above categorization in our related work section more explicit.  

> "Furthermore, this categorization seems counterintuitive -- does it really make sense to state that "distributional treatment effects" and "distribution of treatment effects" are very different?"

Thank you for asking this important question. There are indeed two major differences between these streams:
1. **Interpretation**. _Distributional treatment effects_ represent the differences between different distributional aspects of the potential outcomes [3]. Hence, they can answer questions like “How are 10% of the worst-possible outcomes _with treatment_ different from the worst 10% of the outcomes _without treatment_?”. Here, the two groups (treated and untreated) of the worst 10% contain, in general, **different individuals**. This is problematic in many applications like clinical decision support and drug approval. Here, the aim is not to compare individuals from treated vs. untreated groups (where the groups may differ due to various, unobserved reasons). Instead, the aim is to accurately quantify the treatment response for each individual and allow for quantification of the personalized uncertainty of the treatment effect.    The latter is captured in the distribution of the treatment effect, which allows us to answer the question about the CDF/quantiles of the treatment effect. For example, we would aim to answer a question like “What are the worst 10% of values of the treatment effect?”. Here, we focus on the treatment effect **for every single individual**. The latter is more complex because we reason about the difference of two potential outcomes simultaneously. Hence, in natural situations when the potential outcomes are non-deterministic, both (a) the distributional treatment effect and (b) the distribution of the treatment effect will lend to _very different_ interpretations, especially in medical practice. In particular, the distribution of the treatment effect (which we study in our paper) is important in medicine, where it allows quantifying the amount of harm/benefit after the treatment [4]. This may warn doctors about situations where the averaged treatment effects are positive but where the probability of the negative treatment effect is still large. 
2. **Inference**. The efficient inference of the distributional treatment effects only requires the estimation of the relevant distributional aspects of the conditional outcomes distributions (e.g., quantiles) and the propensity score [1]. However, in our setting of the bounds on the  CDF/quantiles of the treatment effect, we also need to perform sup/inf convolution of the CDF/quantiles of the conditional outcomes distributions. Hence, while the definitions of (a) the distributional treatment effects and (b) the distribution of the treatment effect appear related, their estimation is very different.  
As you can see above, the distributional treatment effects and the distribution of the treatment effect are related to different questions in practice and help in different situations.

**References**:
- [3] Kallus, Nathan, and Miruna Oprescu. "Robust and agnostic learning of conditional distributional treatment effects." International Conference on Artificial Intelligence and Statistics. PMLR, 2023.
- [4] Nathan Kallus. “What’s the harm? Sharp bounds on the fraction negatively affected by treatment”. In: Advances in Neural Information Processing Systems. 2022.
[1/3]

Thank you for the quick response! We sincerely appreciate the time and effort you have taken to provide us with valuable feedback. We are very sorry for the ambiguities in the terminology across existing research research streams and for our original manuscript not providing enough context to resolve them. We further apologize if our previous comments were unclear or could be interpreted as implying a lack of understanding – this was not what we meant, and **we apologize if our response may have come across in the wrong way**. Rather, we feel grateful for having received such thorough and knowledgeable reviews that both commended the quality of our paper and also made further suggestions to improve our work for the camera-ready version.  



### Regarding the unanswered Q3 and W1

> “Using your framework, my questions and concerns (Q3 and W1) were about what practically interesting scenarios "Aleatoric uncertainty of the treatment effect" can capture that other research streams cannot.” 

> “why a community need the proposed estimator, given that the quantile estimator can capture the distributional treatment effect.”

**We apologize for misunderstanding your questions**. In our paper, we work with the aleatoric uncertainty of the treatment effect in the form of the CDF/quantiles of the treatment effect at the covariate-conditional level. More specifically, we aim at inferring the probability that the treatment effect is less than or equal to a certain value ($\delta$), conditional on the covariates. For $\delta=0$, the latter becomes the covariate-conditional probability of the treatment harm (benefit), i.e., $\mathbb{P}(Y[1] \le Y[0] \mid x )$. 

_Why is the distribution of the treatment effect relevant?_ Here are three examples of how the covariate-conditional probability of the treatment harm (benefit) is useful in practice and how other research streams are unable to provide this information:  

1. **Medicine**. In cancer care, for instance, the conditional average treatment effect may suggest whether a treatment is beneficial _on average_ while it can not offer insights into how probable negative outcomes are. For example, consider a patient with a tumor and a drug for which the conditional average treatment effect is larger than zero, suggesting that the treatment has a benefit _on average_ (=  the tumor size reduces on average).  However, the average treatment effect does not tell us how likely such a reduction is. Given the randomness of the potential outcomes, there could be a chance that the tumor size will increase after the treatment. Hence, medical practitioners are often interested in understanding the probability of treatment benefit or harm [1,2], which is captured by the distribution of the treatment effect. This allows us to answer questions such as: _what is the probability that the treatment effect is larger than zero_? In the above example, this means how probable is a reduction of the tumor after treatment? Crucially, such questions cannot be answered by distributional treatment effects and require knowledge of _the distribution of the treatment effect_, which is the focus of our paper.  
   
2. **Public health**. As another practical example, we refer to our case study analyzing the effectiveness of lockdowns during the COVID-19 pandemic. Here, policy-makers are interested in knowing the probability that the incidence after a strict lockdown will be lower than or equal to the incidence without it (=probability of treatment benefit) (see Appendix I).  

3. **Post-approval monitoring of drugs.** Understanding the aleatoric uncertainty of the treatment effect is also relevant when monitoring the efficacy of drugs post-approval. Here, substantial increases in the aleatoric uncertainty serve as an early warning mechanism for when treatments are not working well for certain subgroups of patients or when the pharmacodynamics are not fully understood for all parts of the patient population.  

In sum, there are many examples – especially in medicine – where the distribution of the treatment effect is relevant for practice. Importantly, the distribution of the treatment effect is necessary in order to understand the probability of treatment benefit (or of treatment harm). Below, we also discuss why the distributional treatment effects are very different from the distribution of the treatment effect, and why only the latter can answer the above questions. 

**References**:
- [1] Bordley, Robert F. "The Hippocratic Oath, effect size, and utility theory." Medical Decision Making 29.3 (2009): 377-379.
- [2] Nicholson, Kate M., and Deborah Hellman. "Opioid prescribing and the ethical duty to do no harm." American journal of law & medicine 46.2-3 (2020): 297-310.
Thanks for clarifying my question!

I read the other reviews and the authors' response. I agree with some feedback asking to further clarifying the setting and how it compares/differs to related settings. I have no concerns about the validity and novelty of the method, and stick to my rating.
I thank the authors for their helpful response and maintain my current score.
Overall, I am embarrassed by this response. The response is written as though the authors' categorization—(1) Aleatoric uncertainty of potential outcomes, (2) Distributional treatment effects, and (3) Aleatoric uncertainty of the treatment effect—has been clearly established in the paper, and my questions/concerns arise from a lack of understanding of this framework (e.g., " the distributional treatment effects are not the target of our paper,", where the term "distribution treatment effect" is just defined by this response, not the paper). 

However, as stated in the authors' response, this categorization is not mentioned in the paper, and no clues are provided regarding this categorization.

Furthermore, this categorization seems counterintuitive -- does it really make sense to state that "distributional treatment effects" and "distribution of treatment effects" are _very different_? I do understand that they are different, and that the distribution of the treatment effect, P(Y[1] - Y[0]), is not pointwise identifiable due to the fundamental problem of causal inference. However, the categorization seems to come out of nowhere and requires justification.

I do understand which research streams this work is located in. Using your framework, my questions and concerns (Q3 and W1) were about what practically interesting scenarios "Aleatoric uncertainty of the treatment effect" can capture that other research streams cannot. I don't believe your response has fully addressed this question/concern yet.

I adjusted my score based on this response.
I would like to thank the authors for their detailed response to my review.

My concerns and questions were well addressed. I believe that the revised version of the paper, with the incorporated changes, will be easier for readers to follow.

As a result, I am raising the score from 6 to 7.
I thank the authors for their response, the comparison to previous bounds cleared my doubts. I maintain my original score of accept.