{'fraction_non_alphanumeric': 0.05689522400229514, 'fraction_numerical': 0.015250577559001615, 'mean_word_length': 4.21152030217186, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sensitivity Analysis under the f -Sensitivity Models: Definition, Estimation and Inference

The digitization of the economy has witnessed an explosive growth of available observational data across a variety of domains, which has provided an exciting opportunity for causal inference that aims to discern the effectiveness of different treatment options. However, unlike in experimental data – where the designer has the full control to randomize treatment choices based on observable covariates – causal inference using observational data faces the key challenge of confoundedness: there might exist unobservable covariates that affect both the treatment choice and the outcome. In this paper, we aim to add to the emerging line of tools for sensitivity analysis that assesses the robustness of causal conclusions against potential unmeasured confounding. We propose a new sensitivity model where, in contrast to uniform bounds on the selection bias in the literature, we assume the selection bias is bounded “on average”. It generalizes the widely used Γ-marginal selection model of Tan (2006) and is particularly suitable for situations where the selection bias grows unbounded locally but, due to the impact of such region being small overall, is controlled at the population level. We study the partial identification of treatment effects under the new sensitivity model. From a distributional-shift perspective, we represent the bounds on counterfactual means via certain distributionally robust optimization programs. We then design procedures to estimate these bounds, and show that our procedure is doubly robust to the estimation error of nuisance components and remains valid even when the optimization step is off. In addition, we establish the Wald-type inference guarantee of our procedure that is again robust to the optimization step. Finally, we demonstrate our method and verify its validity with numerical experiments.


Introduction
Causal inference using observational data is a necessity in many fields for several reasons. First, experimental data is often limited and costly to obtain. Second, modern advances in digitization technology have brought forth a growing availability of observational data across a variety of domains. Third, even in settings where experimentation is an economic possibility, observational data can help support initial inquiries and inform effective designs of experiments, by, for instance, identifying the unknown/least confident regimes that require further experimental probing. As such, developing efficient -both statistically and computationally -causal inference schemes using observational data has immense practical value.
In contrast to experimental data where the experimenter has the luxury to randomize treatment options, the key difficulty with observational data is the (untestable) presence of unobserved confounding covariates (also known as confounders). A confounder is a variable U that affects both the treatment T and the outcome Y , and whose presence could invalidate the causal relationship drawn from the inference. For instance, one might observe that patients that are sent to ICUs (intensive care units) have a high death rate. However, one would be wrong to conclude that admission to ICUs (a treatment option) causes high likelihood of death, because the patients who are sent to ICUs are those who are already very sick. As such, in this example, the patient condition is a confounder, 1 which impacts both the treatment selected (ICU/non-ICU) and the medical outcome (death/recovery): ignoring this confounder will lead to faulty conclusions.
To address this challenge, we start by describing the problem of causal inference with observational data and stating the assumptions that will be used throughout the paper. We follow the potential outcome model (Rosenbaum and Rubin, 1983b;Rubin, 1990) and posit an underlying data-generating distribution P on (X, U, T, Y (1), Y (0)), where X ∈ X is the observed covariate vector and X is compact, U ∈ U is the unobserved confounding vector, T ∈ {0, 1} is the treatment option, 2 and Y (1) ∈ R and Y (0) ∈ R are the two potential outcomes. We assume access to a dataset {(X i , T i , Y i )} n i=1 of n triples, where X i is unit i's covariate vector, T i is the treatment applied to unit i and Y i = Y i (T i ) is the observed outcome for unit i under treatment T i . Note that although there are two potential outcomes Y i (0) and Y i (1), only Y i (T i ) -the potential outcome corresponding to the actual applied treatment T i -is observed. 3 Note that without loss of generality (since U is arbitrary), P satisfies We are interested in estimating the following quantities: where in all quantities the expectation is taken with respect to the underlying joint distribution. To make progress, we assume that there is sufficient exploration in the dataset, known as the overlap assumption in the literature. We define the observed propensity score e(x) P(T = 1 | X = x).
Assumption 1 (Overlap). 0 < e(x) < 1 for P-almost all x ∈ X . 4 Under the overlap assumption, classical causal inference (see Imbens (2004); Imbens and Rubin (2015) for two excellent surveys) has extensively treated the topic of estimating ATC and ATT (along with other treatment effect quantities) under the unconfoundedness assumption (a.k.a. strong ignorability (Rosenbaum and Rubin, 1983b)), where it is assumed that Y (1), Y (0) |= T | X (i.e. no confounders U exist). As previously mentioned, this assumption is not testable and hard to justify in observational data, where the statistician often has little knowledge and/or control over how the dataset was collected. The key technical challenge arising therefrom is that the counterfactual distribution P X,Y (1) | T =0 (and similarly P X,Y (0) | T =1 ) is not identifiable using observed data. To address this problem, a prominent approach, dating back to Cornfield et al. (1959); Bross (1966), is sensitivity analysis, which studies how robust a certain causal conclusion drawn from observational data is against unmeasured confounding. In particular, a sensitivity model often makes assumptions on how much the treatment selection probability when conditioned on the confounders deviates from the treatment selection probability when unconditioned on the confounders. From such an assumption on the deviation, the statistician then proceeds to estimate the range -rather than a single value -of the underlying treatment effect quantity, where the width of the range depends on the deviation bounds. Such an approach is valuable in that it provides the statistician with a quantitative understanding how much the treatment effect is under varying severities of confoundedness.
Within sensitivity analysis, Rosenbaum's Γ-selection sensitivity model (Rosenbaum, 1987) provides a pioneering model on the selection bias that has instantly become a classic. Subsequently, it was generalized by the marginal Γ-selection sensitivity model proposed in Tan (2006), which imposes weaker assumptions and is widely adopted. Based on Tan (2006), a series of work have developed various treatment effects estimation and inference schemes (Zhao et al., 2017;Kallus et al., 2019;Lee et al., 2020;Nie et al., 2021;. More specifically, Tan (2006) assumes uniform upper and lower bounds on selection bias, in terms of the odds ratio of receiving treatment conditional on both unmeasured confounders and observed covariates versus conditional only on observed covariates, specified as 1 Γ ≤ P(T = 1 | X = x, U = u)/P(T = 0 | X = x, U = u) for P-almost all x ∈ X and u ∈ U for some Γ ≥ 1. Note that when Γ = 1, this assumption recovers the unconfoundedness assumption, and the larger the Γ, the more confoundness the model tolerates.

Motivation of our approach
Despite being widely used, the marginal Γ-selection model (2) can be limited. We illustrate this point by a simple and natural parametric example, which further motivates our new sensitivity model.
Example 1. Let us consider a simple example where there are no covariates, and the statistician would like to infer the counterfactual mean of Y (1) in the control group (i.e. E[Y (1) | T = 0]). We assume the observed probability of treatment is P(T = 1) = 1/2. Instead of believing in the strong ignorability condition (that is, all units receive treatments with the same probability in our example), the statistician would like to estimate the range of E[Y (1) | T = 0] if the observational data is confounded to some extent. She might think of a logistic model T | U ∼ Bern exp(δU ) 1+exp(δU ) for some δ ∈ (0, 1), where the confounder U ∼ N (0, 1); hence, where the second-to-last inequality follows by symmetry of U , and agrees with the observed probabilities. For simplicity, we also assume that Y (1) = g(U, δ) for some measurable function g, so that Y (1) | T = 1 agrees with the observations regardless of the choice of δ (and hence also Y (1) |= T | U , satisfying (1)). By construction, the odds ratio characterizing the selection bias caused by U is OR(U ) := P(T = 0 | U ) P(T = 1 | U ) · P(T = 1) P(T = 0) = e −δU .
With an unbounded U , the above odds ratio can be unbounded and hence the marginal Γ-selection assumption (2) does not hold. Consequently, estimating the range of counterfactual means under a hypothesized uniform bound on OR(U ) is problematic even in this simple stylized example. More generally, this example highlights a critically inadequate aspect in the marginal Γ-selection model: it may well be that some region of the sample space could be drastically confounded; considering such a possibility, the statistician might be reluctant to imagine a uniform bound on the selection bias. However, one could still imagine that even though the selection bias can have a huge scale in some region, the impact of such region might still be small if it only takes a tiny proportion of the whole population. As such, the locally-severe confounding might not have an overall big impact at the population level. A desirable sensitivity model should still allow us to conduct estimation and inference in such scenarios. We do mention that the earlier work Imbens (2003) postulates a parametric logsitic model for treatment that includes the confounders, which contains the simple example discussed here as a special case and hence allows for unbounded local confounded effects. However, its key limitation is that the proposed confounding model is highly specialized to the logistic form, whereas the marginal Γ-selection criterion provides a nonparametric model that is quite general, which is an important merit. As such, motivated by achieving the best of both worlds, we are naturally led to the following research questions: Can we formulate a general, non-parametric sensitivity model that can deal with unbounded local confounding effects with an overall bounded impact at the population level? If so, can we efficiently-both statistically and computationally-conduct estimation and inference on average treatment effects?

Our proposed f -sensitivity model
We answer the above questions in the affirmative in this paper. In this subsection, our focus is on providing a quick answer to the first question, where we briefly describe our proposed f -sensitivity model, characterized by the following (f, ρ)-selection condition.
for P-almost all x.
Remark 1. If P satisfies the marginal Γ-selection condition (2), then it automatically satisfies the (f, ρ)selection condition with any qualified f and ρ = max{f (1/Γ), f (Γ)}. This can easily be checked by noting that (2) implies by convexity of f , thereby leading to the bound on the second integral in d f (P): The first integral in d f (P) can be bounded by exactly the same way by the symmetry of (2).
Remark 2. In the definition of the (f, ρ)-selection condition, we take the maximum of two integrals, each from one direction. This is mainly to keep the condition symmetric with regards to the choice of the treated or control groups, and align with the convention in the sensitivity models in the literature (Tan, 2006). That said, as we would see shortly in Section 2, it might be more natural to only work with one of them (i.e. assume one of them be bounded by ρ) when our primary interest is to estimate and perform inference on one of the counterfactuals.
We now show that the proposed (f, ρ)-selection model addresses the unbounded confoundness in Example 1: even though the selection bias (the odds ratio) is not uniformly bounded, it is controlled overall; in this case" our proposed model can be a better characterization of the selection bias.
Example 1 (Continued). We take f (t) = t log t, a convex function with f (1) = 0. Continuing the computation in Example 1, the first term in d f (P) of Definition 1 can be computed as The right-handed side is approximately 0.2 if we take δ = 1 and 0.6 if we take δ = 2. Note that f OR(U ) dP U | T =1 can be interpreted as the overall deviation of OR(U ) from 1 in the treated (observed) group, which is bounded, even though OR(U ) → ∞ when U → −∞.
While the uniform bound on selection bias is infeasible in this case, the procedures proposed in this paper still yield meaningful bounds on the counterfactual mean E[Y (1) | T = 0]. We apply our procedure to Example 1 for a sequence of δ and the corresponding upper bounds in (3). The results are summarized in Figure 1, where we assume Y (1) | T = 1 ∼ N (0, 1) for simplicity.

Our contributions
Our contributions are threefold. First, we propose the f -sensitivity model, a general non-parametric model that allows for unbounded local confounding with a limited overall impact at the population level. In Section 2, we show a distributional shift interpretation of our f -sensitivity model that is perhaps more illuminating: under the (f, ρ)-selection condition, the (unobservable) counterfactual distribution P Y (1) | X=x,T =0 is close to the observable distribution P Y (1) | X=x,T =1 under the f -divergence measure, whose precise definition is provided shortly. This allows us to represent the partial identification lower (resp. upper) bound on the unidentifiable E[Y (1) | T = 0] based on the observable (identifiable) distribution P Y (1) | X=x,T =1 with a distributionally robust optimization problem we provide in Section 3.1, where we infer the smallest (or largest) objective within the f -divergence ball of radious ρ centered at P Y (1) | X=x,T =1 .
Building on this distributional shift perspective, we provide in Section 3 an estimation procedure for the bounds on counterfactual means and treatment effects. We first transform the distributionally-robust optimization problem into an equivalent dual problem, through which the bounds are in fact connected to the optimal objective of a risk minimization problem. We then utilize an empirical risk minimization (ERM) procedure to solve for the risk minimizers. Finally, we propose a cross-fitted estimator for the bounds, where we adjust for the estimation error of the plugged-in empirical risk minimizer and other nuisance components.
Third, in Section 4, we establish the estimation and inference guarantees of our procedures for counterfactual means, demonstrating certain robustness properties. We show that our estimators are doubly robust to the nuisance estimation when the ERM step is consistent, and remain valid even when the ERM step is misspecified. Furthermore, we show that our cross-fitted estimator yields Wald-type inference guarantees even when all the nuisance components can only be estimated at a slower-than-parametric rate; such inferential guarantee is again robust to the misspecification of the ERM step. Based on the the bounds on counterfactual means, we provide confidence intervals for the treatment effects.

Related work
In addition to the sensitivity analysis literature mentioned before, several other works exist that take the approach of modelling the impact of unmeasured confounding is through bounds on selection bias. Remarkably, Rosenbaum and Rubin (1983a) studies the impact of selection bias among matched pairs, which is further extended by a series of works (Rosenbaum, 1987;Gastwirth et al., 1998;Rosenbaum, 2002b,a) to sensitivity models that uniformly bounds the selection bias among samples with matching covariates.
Also related to sensitivity analysis under uniformly bounded selection bias is Yadlowsky et al. (2018), which works under an extention of Rosenbaum's sensitivity model that is similar to Tan (2006). Besides modelling the selection bias among treatment, observed covariates and unmeasured confounders, we mention in passing that Ding and VanderWeele (2016) considers the sensitivity analysis of other metrics than treatment effects.
Our novel sensitivity models provide a perspective of studying the distributional shift from observations to the counterfactuals, which echos the ideas of several previous works Yadlowsky et al., 2018;. In particular, we formulate the estimand via an optimization problem under constraints on distributional shifts, similar to Yadlowsky et al. (2018); . However, as they work under different sensitivity models, the techniques employed in our works are in sharp distinction from theirs.
The f -divergence bounds derived from our sensitivity models are related to a series of works that characterize discrepancy between distributions using f -divergence (Rényi, 1961;Morimoto, 1963;Csiszár, 1964;Liese and Vajda, 2006). For example, Christensen and Connault (2019) places bounds on the marginal distributions of hidden variables in structural equation models; Andrews et al. (2020) studies parameter estimation when the joint distributions of variables are within an f -divergence ball; Rahman (2016) estimates the f -divergence between conditional and unconditional distributions to measure variable importance; we differ from them both in terms of the target (bounding selection bias), the estimand (estimation of tretment effects) and the resulting formulation (f -divergence between conditional distributions, while the shifts between X-distribution are identifiable). In particular, these existing f -divergence bounds are not sufficient to capture the impact of unmeasured confounding in our problem setting.
Our work also finds connection to a line of work estimation, inference and learning under various types of distributional shifts, e.g., f -divergence, outside the task of inferring causal effects, including Duchi and Namkoong (2021) which studies the empirical risk minimization problem when the joint distribution of (X, Y ) shifts within an f -divergence ball, Si et al. (2020) which studies the policy learning for contextual bandits under unknown marginal distribution shifts, Gupta and Rothenhäusler (2021) that studies the estimation and inference of statistical parameters under distributional shifts in certain directions, etc. As we work under distinct type of (conditional) distributional shifts, we need a different set of tools than theirs.
Finally, the post-hoc adjustment method we propose to deal with nuisance components is also generally connected to a vast body of missing data literature with unknown missing mechanisms, although in different contexts and with different details. In particular, we use the cross-fitting (Schick, 1986;Zheng and Laan, 2011;Chernozhukov et al., 2018) technique to mitigate the error in nuisance component estimation; the technique of bias correction using a different dataset under covariate shift is also related to and inspired by Jin and Rothenhäusler (2021).

Distributional Shifts under the f -Sensitivity Model
The f -sensitivity model leads to contraints on the distributional shifts from the observations to the counterfactuals; this identifies a new class of robust inference problems that are new to the literature. We can thus view sensitivity analysis through the lens of distributioanlly-robust inference. Take the counterfactual inference on Y (1) for T = 0 as an example (the case of Y (0) is completely symmetric): since we are able to observe Y (1) in the treated group, the distribution of (X, Y (1)) | T = 1 is identifiable, and decomposes as The counterfactual distribution of (X, Y (1)) in the control group is where part (a) is identifiable but part (b) is not without the strong ignorability condition. Put another way, the distributional shift from P X,Y (1) | T =1 to P X,Y (1) | T =0 has an identifiable component as a function of X, and an unidentifiable one concerning the conditional distribution of Y | X; the latter, however, is bounded under the (f, ρ)-selection condition in terms of f -divergence defined as follows.
Definition 2 (f -divergence). Let P and Q be two probability distributions over a space Ω such that P is absolutely continuous with respect to Q. For a convex function f such that f (1) = 0, the f -divergence of P from Q is defined as Popular examples for f -divergence in the literature include the Kullback-Leibler (KL) divergence with f (t) = t log t, the total variation (TV) distance with f (t) = |t − 1|/2, Pearson χ 2 -divergence with f (t) = (t − 1) 2 , and the Cressie-Read family of f -divergences (Cressie and Read, 1984) parametrized by k, where f k (t) := t k −kt+k−1 k(k−1) . Throughout this paper, we work with generic forms of f -divergence, and provide discussions on concrete examples where proper conditions are satisfied for our methodology.
for P X | T =1 -almost all x; that is, the f -divergence between the conditional distributions in the two groups are bounded by ρ for almost all X in group T = 1.
Proof of Lemma 1. Suppose a distribution P over the (X, U, T, Y (0), Y (1)) satisfies the (f, ρ)-selection condition. By the data-processing inequality, Above, the likelihood ratio can be decomposed as where step (a) is due to condition (1). Combining the above, we have where the last inequality is due to the (f, ρ)-selection condition.
In Proposition 1, we characterize the counterfactual distributions induced by all super-populations that agrees with the observables and satisfies our sensitivity models.
Proposition 1. Let P sup be the true unknown super-population over (X, U, T, Y (0), Y (1)) and let P be the set of all distributions over (X, U, T, Y (0), Y (1)). Let P obs X,Y,T be the joint distribution of all observable random variables (X, Y, T ). Define Q 1,0 to be the set of all counterfactual distributions that agrees with the observables and satisfies the (f, ρ) selection condition, i.e., Q 1,0 = P X,Y (1) | T =0 : P ∈ P, P X,Y,T = P obs X,Y,T , P satisfies Definition 1 .
Then P sup X,Y (1) | T =0 ∈ Q 1,0 , and where r 1,0 (x) = (1−e(x))p1 e(x)(1−p1) , and e(x) = P obs (T = 1 | X = x), p 1 = P obs (T = 1). We defer the proof of Proposition 1 to Appendix A, where we prove a stronger version of Proposition 1 that gives a tight characterization of Q 1,0 . Symmetrically, we can also define Q 0,1 as the identification set of P X,Y (0) | T =1 ; the tight characterization of Q 0,1 is also given in Appendix A. Throughout this paper, we mainly work under the ambiguity sets in Proposition 1 to emphasize the more general distributionally-robust estimation aspect of this problem.
Proposition 1 identifies a new class of robust inference problems, where the target distribution has an identifiable X-shift and the unidentifiable conditional distribution is restricted in f -divergence ball; it is similar to the ambiguity set studied in . This model is closely related to, but quite distinct from other ambiguity sets considered in the literature that involves f -divergence (Duchi and Namkoong, 2021;Si et al., 2020;Andrews et al., 2020) that only concern marginal distributions.
Remark 3 (Relation to other f -divergence bounds). Previous works in the literature (see Section 1.4 for a summary) often work under the f -divergence ball around the marginal distribution P X,Y , characterized by While in our formulation, the ambiguity sets take the form where r(x) is a known or identifiable function. Such distinction is similar to Remark 3 of . Instead of bounding the overall shift as in (4), the constraint in (5) actually allows freedom in the shift of X: the sets in (5) can be small as long as ρ is small. For counterfactual inference under the strong ignorability condition, the set (5) is a singleton even if P X and P X are drastically different, while (4) would require a large ρ to hold. More generally, when there is a known (or estimable) large shift in P X but a relatively small shift in P Y | X , (5) provides a tighter range of the target distributions, and the methods we develop in this paper can be directly applied.

Bounding treatment effects
We first focus on the bounds on counterfactual means, which serve as building blocks for bounding various treatment effects. For example, is the likelihood ratio that lies within the ambiguity set in Proposition 1, and P X,Y (1) | T =1 is a distribution that one has access to. Proposition 2 explicitly characterizes the bounds on under the (f, ρ)-selection condition; they are provided as solutions to convex optimization problems that only involve identifiable quantities.
Proposition 2. Let µ − 1,0 (resp. µ + 1,0 ) be the optimal objective function of the convex optimization problem min(resp. max) where all the expectations are induced by the observed distribution.
Proposition 2 immediately implies bounds on the ATC: denote the observable group-wise means as µ obs Switching the role of 1 and 0 in Proposition 2, one can obtain bounds on E[Y (0) | T = 1]; let µ + 0,1 and µ − 0,1 denote the upper and lower bound on E[Y (0) | T = 1], respectively, we then get bounds on the ATT: By the decomposition of ATE (average treatment effects), we also have the representation of lower and upper bounds for Under the (f, ρ)-selection condition, we have These bounds thus boil down to the estimation of µ ± t,1−t under our sensitivity models, the optimal objective value of the convex optimization problems in Proposition 2.
Remark 4. As we mentioned before, optimal objectives of the problems in Proposition 2 are not necessarily tight bounds for counterfactual means. To align with the literature and keep a relatively clean formulation of dual problems, we only account for the direction For completeness, we include the discussion of the tight bounds on counterfactual means, hence ATT and ATC, in Section 6. We also note that combining sharp bounds on ATT and ATC does not necessarily lead to sharp bounds on ATE, as they might be attained by different super-populations. We leave the investigation of sharp bounds on ATEs for future pursuit.

From the primal to the dual
The infinite-dimensional optimization problem (6) is hard to solve directly. We address this issue by transiting to its dual form that is easier to tackle with. In the following, we primarily focus on µ − 1,0 , the lower bound on E[Y (1) | T = 0], and the same idea carries over to other quantities. Proposition 3 represent µ − 1,0 with a dual formulation, whose proof is in Appendix C.1.
Proposition 3. The optimal objective of (6) is given by The formulation (9) can be viewed as a risk minimization problem, where the minimizer (α * (x), η * (x)) is also the minimizer of per-x conditional risk; as a result, we might employ tools in empirical risk minimization to estimate α * (·), η * (·), and then estimate µ − 1,0 by plugging in the empirical risk minimizer. However, when it involves functions α(·), η(·), the estimation might not converge at a parametric rate; furthermore, another particular challenge is the nuisance component r 1,0 (x) that relies on e(x) = P(T = 1 | X = x), which typically needs to be estimated in observational studies. Similar to the situations for treatment effect estimation without confounding, naive plug-in estimator might suffer from slow rate of convergence. We are to show that, with an adjustment term to the naive plug-in estimator that utilizes covariates in the control group, we can still obtain Wald-type inference when the nuisance components are estimated at a lower-than-parametric rate.
Before proceeding to our procedures, we show a result on the behavior of the optimizer α * (x): it is positive as long as P Y (1) | X,T =1 does not have a large point mass at its essential infimum and the function f satisfy some regularity conditions in the limit. The proof of Proposition 4 is deferred to Appendix B.1.
In particular, the conditions on f hold for a large variety of functions; concrete examples include KL divergence, where f (x) = x log x and f * (x) = e x−1 , as well as In the following, we assume throughout that the conditions of Proposition 4 hold, hence there exists some > 0 such that α * (x) > for P X | T =1 -almost all x by the compactness of X . An important implication is that, as α * (x) lies in the interior of [0, ∞), the gradient of the risk function is typically mean-zero at (α * (x), η * (x)); this would play an important role in our inferential guarantees.

The estimation procedure
We start with splitting the treated and control groups into three folds, denoted as I respectively. For each j = 1, 2, 3, we use the fold I to obtain an estimator r (j) for r 1,0 , and run an empirical risk minimization problem to obtain α (j) and η (j) for (α * , η * ) on the same fold (the estimation details are to be discussed right after the introduction of the procedure). In particular, the optimization of α (j) and η (j) does not require the knowledge of r 1,0 here. We then define the function Then on fold I (j+2) 1 , we run a regression algorithm to obtain an estimator h The above procedure is repeated for each j = 1, 2, 3, and we average the three estimators to obtain The estimation procedure is summarized in Algorithm 1.
Algorithm 1 The estimation procedure for µ − 1,0 Input: The set of treated samples I 1 ; the set of control samples I 0 ; the algorithm R for estimating r; the optimization algorithm O solving for (α, η); the regression algorithm Reg for obtainingh.

Solving for α(·) and η(·)
We take a moment to elaborate on the estimation of θ (j) := ( α (j) , η (j) ). From Proposition 3, θ * := (α * , η * ) is also (i.e., removing r 1,0 (x)) the population risk minimizer of The empirical risk is correspondingly (recall that we run with fold I We can thus consider a function class Θ, and solve for the empirical risk minimization (ERM) problem. While the downstream inference is different, this approach is similar to Yadlowsky et al. (2018), where they express the bounds of conditional expectations of counterfactuals themselves as empirical risk minimizers.
To be specific, we use the method of sieves (Geman and Hwang, 1982); we consider an increasing sequence Θ 1 ⊂ Θ 2 ⊂ · · · of spaces of smooth functions, and let We consider the following two examples inspired by Yadlowsky et al. (2018), where we also truncate the functions away from zero for α(x): we note that, if α * (x) is always positive (implied by the minimality of the risk function) and continuous (satisfied if P Y (1) | X=x,T =1 is smooth in x) and X is a compact set, then there exists a positive > 0 such that inf x∈X α * (x) ≥ . In practice, we can set to be small enough, or let = n decays slowly to zero; this does not hurt the capability of function class or the convergence.
Example 2 (Polynomials). Let Pol(J) be the space of J-th order polynomials on [0, 1]: and let Pol(J, ) be the space of J-th order polynomials on [0, 1] truncated at > 0: Then we define the sieve Example 3 (Splines). Let 0 = t 0 < · · · < t J+1 = 1 be knots that satisfy max 0≤j≤J (tj+1−tj ) min 0≤j≤J(t j+1 −t j ) ≤ c for some c > 0. We define the space for r-th order splines with J knots as and the truncated space for r-th order splines with J knots as We consider the classes of sufficiently smooth functions; for p 1 = p − 1 and p 2 = p − p 1 , we define To ensure non-negativeness, we also define the truncated function class Λ p c (X , ) := x → max{f (x), } : f ∈ Λ p c (X ) , obtained by thresholding Λ p c (X ) away from zero. To obtain convergence of the estimators, we impose the following assumptions on the true optimizer and regularity conditions of the loss function.
We include a detailed discussion of Assumption 1 in Appendix B.2, where we provide concrete examples and justifications. In Assumption 1, we assume the true optimizer is sufficiently smooth, so that function approximator can learn it well. It can be satisfied if the conditional distribution P Y (1) | X,T is sufficiently "smooth" in x. We require the strong convexity of the conditional risk function at its minimizer θ * (x); it is typically the case if Y (1) is not deterministic given X. The stability condition at θ * (x) can be satisfied if Y is not heavy-tailed. We also assume that the population risk is stable in terms of Under the above regularity conditions, we obtain convergence rates of the empirical risk minimizers ( α (j) , η (j) ). The proof of the following theorem is in Appendix C.2.
Theorem 1. Suppose Assumption 1 holds. We set J n = ( log n n ) 1/(2p+d) for the sieve estimators in Examples 2 and 3, and suppose θ . Then employing the function classes given in Examples 2 or 3, we have θ The above theorem shows that under reasonable smoothness of the optimizer and regularity conditions on the loss function, the empirical risk minimizer α (j) , η (j) converges to the truth. Besides the examples and guarantees we provide, similar results might be obtained for other function classes like wavelets (Daubechies, 1992), and the conditions in Assumption 1 might be weakened or modified to account for more generality. Such extension is beyond the scope of this work.

Theoretical Guarantees
In this section, we provide the theroetical guarantees of our procedure in Section 3.3. We first show the double robustness and one-side validity regarding consistency. We then present inferential guarantees: we achive Wald-type inference for µ − 1,0 under slower-than-parametric convergence rates of the nuisance component estimation; moreover, even when the empirical risk minimization is not consistent to the optimum, our inference procedure can still be valid. Finally, we show how to leverage our procedure to construct bounds for treatment effects.

Double consistency and one-side validity
Our first result is regarding the consistency of our estimator: we show that µ − 1,0 is doubly robust to nuisance estimation. Even more interestingly, our estimator is robust to the ERM step: given that either r (j) or h (j) is consistent, it converges to the true bound if the ERM step is consistent; otherwise, our estimator converges to a conservative but still valid lower bound of µ − 1,0 . We call this "one-side validity". We start with a mild assumption on the convergence of ERM step; note that we do not assume the convergence to the true minimizer (α * , η * ). For notational convenience, we denote the risk function (θ, x, y) = af * y+b −a + b + aρ for θ = (a, b) as in Proposition 3. When there is no confusion, we equivalently use (θ, x, y) = ((α(x), η(x)), x, y) when θ = (α, η) is a function.
In Assumption 2, we additionally assume a mild regularity condition on the first-order expansion at the limit; it is satisfied if the loss function is differentiable or locally Lipschitz. The second moment condition is also mild and standard. The following theorem shows the double robustness as well as one-side validity of our estimator, whose proof is in Appendix C.3.
Theorem 2. Suppose Assumption 2 holds for some fixed θ = (α , η ). Assume either ( . Then the following holds: if θ = θ * , i.e., the ERM step is consistent, then µ − 1,0 = µ − 1,0 + o P (1); otherwise, µ − 1,0 = µ 1,0 + o P (1) for some constant µ 1,0 ≤ µ − 1,0 . A similar property of one-side validity has been documented by a recent work of , where they work under the marginal sensitivity model of Tan (2006) and develop a similar property based on an exact characterization of the worst-case scenario. However, in our setting, the one-side validity is a relatively straightforward consequence of the interplay between primal and dual problems. It would be interesting to find connections between our results; for example, whether their result can also be implied by the duality.
4.2 Wald-type inference for µ − 1,0 We now turn to inferential guarantees. We show that our procedure yields valid Wald-type inference under slow convergence rates of nuisance estimations. We begin with some regularity conditions on the risk function.
In Assumption 3, we require the risk function to be differentiable and admits a Taylor expansion near some optimizer, as well as a regularity condition on the exchangeability of differentiation and conditional expectation. These are mild conditions that are commonly adopted in the literature ( Van der Vaart, 2000). The risk function is assumed to be stable, so that plugging in estimators of α * , η * won't cause large errors, which is also a mild condition that can be satisfied under a first-order Taylor expansion condition.
We assume the following convergence rates, where we assume the ERM step is o P (n −1/4 ) consistent, and the nuisance estimation error of r (j) and h (j) has a product of order o P (n −1/2 ).
In Assumption 4, the rate of r (j) depends on the estimation of e(x) = P(T = 1 | X = x), a standard regression problem. The estimation of h (j) is also a regression problem viewing H (j) as fixed. Convergence rate guarantees for such conditional mean estimation problems are well-established in the literature (Stone, 1982;Mallat, 1999;Pagan et al., 1999;Shen and Wong, 1994;Wasserman, 2006;Simonoff, 2012). The estimation of ( α (j) , η (j) ) has been discussed in Section 3.3.
Under the above two assumptions, we show that our estimator is asymptotically normal, where the estimation error of nuisance component is negligible. The proof of Theorem 3 is deferred to Appendix C.4.

Theorem 3. Suppose Assumptions 3 and 4 hold. Then
Here p 1 = P(T = 1) = 1 − p 0 , and we define H(x, y) = (θ * , x, y), h(x) = E H(X, Y (1)) X = x, T = 1 . All the expectations (variances) are induced by the observed distribution. Furthermore, define Similar results can be obtained for µ + 1,0 , if we simply flip the sign of Y (1) and flip back after running the same procedure. The above procedure can also be generalized to the inference of µ ± 0,1 ; the simplest way might be just switching the two groups. We summarize these results in Appendix B.3 for completeness.

Robustness to misspecification of ERM
Our inferential guarantee in Theorem 3 relies on consistency of both the nonparametric regression and the ERM steps. While these conditions are relatively mild, in this part, we take a step further and note that our estimator is in particular robust to the ERM step.
In theorem 4, we only require the convergence of ( α (j) , η (j) ) in L 2 (P · | T =1 )-norm any pair of fixed functions. This might happen, for example, if the function class we employ does not approximate (α * , η * ) very well, but our estimators still converge to a fixed in-class risk minimizer. In this case, our estimator converges to a conservative lower bound of the counterfactual mean and still yields valid inference.
In parallel to the mean-zero gradient property of θ * , we assume a local first-order condition for θ restricted to Θ n , which is crucial for the double robustness to the estimation error. This condition is satisfied as long as θ is the population risk minimizer (with weight r(X)) among Θ n . To obtain an estimator that converges to θ , we might slightly change the procedure: fit r (j) (x) on one fold and and run the ERM with the fitted r (j) (x) on a new fold. The convergence of the empirical risk minimizer can be satisfied if Θ n is not too complex and r (j) is consistent with a slow rate.
We also note that, as implied by Theorem 4, plugging in any fixed function into our procedure without ERM (or equivalently, setting Θ n = {θ} for some fixed θ that satisfy the regularity conditions) also yields a valid lower bound. However, this is uninteresting as it may be way too conservative.

Inference for treatment effects
With the above estimator for counterfactual means in place, we briefly discuss the construction of confidence intervals for treatment effects. Let us first start with ATT/ATC. Following the preceding example of µ − 1,0 , in view of (8), we can construct an estimator for the lower bound of ATC, defined as where µ − 1,0 is constructed as in Section 3.3. Theorem 3 directly implies the following result of double robustness and asymptotic normality for τ − ATC , and the proof is omitted for brevity. Corollary 1. Under the same conditions of Theorem 3, is a lower bound for ATC under the (f, ρ)-selection condition, and Similar to Theorem 3, a consistent estimator σ 2 ATC,− can also be constructed for Var(φ − ATC (X i , Y i , T i )), enabling Wald-type inference. Based on the results of Section 4, we can similarly construct other bounds for ATT and ATC and combine them to obtain bounds on ATE. For example, let µ + 0,1 estimate an upper bound on E[Y (0) | T = 1] with influence function φ + 0,1 (see Appendix B.3 for details). We may construct , where τ − ATE is a lower bound for ATE under the (f, ρ)selection condition, and the influence functions are φ

Numerical Experiments
We illustrate the performance of our procedure on simulated datasets. We focus on the estimation of the counterfactual mean E[Y (1) | T = 0] given confounded observational data and take f (t) = t log t.
In this setting, the confounder is entirely driven by U := Y (1). The odds ratio is and we obtain an upper bound for the f -divergence as ρ = δ 2 /2. The same bound can be obtained for the other odds ratio of the control group. The observed dataset is thus Intuitively, δ drives the direction and magnitude of confounding: when δ > 0, larger values of Y (1) has larger probability of getting treated even conditional on X; as a result, the observed Y (1) in the treated group is actually shifted to larger values, leading to overestimate of treatment effects if confounding is not accounted for. The larger δ is, the more severe the impact of confounding is. On the other hand, when δ < 0, inference under the strong ignorability assumption tends to underestimate the treatment efects. In this setting, although we anticipate OR(X, U ) to be controlled overall, it does not admit a uniform upper bound; we plot several quantiles of OR(X, U ) in the treated group in Figure 2. We apply Algorithm 1 to obtain bounds and confidence intervals of E[Y (1) | T = 1]. The detailed implementation is as follows: the regression algorithm we use for both r and h is Random Forest Regressor from scikit-learn Python library (Pedregosa et al., 2011); we use cubic spline in Example 3 to approximate α * , η * , where we threshold at = 0.001 to guarantee positiveness of α * (yet our estimates turn out to be strictly larger than this threshold). We employ the Nelder-Mead optimizer implemented in SciPy Python library (Virtanen et al., 2020) to optimize the coefficients in the spline approximation.

Sensitivity analysis with one dataset
We first illustrate the estimators and confidence intervals we obtain under a fixed confounded data generating process. To be specific, we fix δ = 0.5 (hence ρ = 0.125) to generate the data, and apply our procedure to the fixed dataset for a series of ρ ∈ {0.05, 0.1, . . . , 0.95, 1.0}. We obtain 0.975-lower confidence bound (LCB) for the lower bound of ATC and 0.975-upper confidence bounds (UCB) for the upper bound of ATC (i.e., the bounds in (7)), which together form a 0.95-CI for ATC under a hypothesized confounding level ρ. The results are plotted in Figure 3.
Without accounting for confounding, reweighting on the covariates tend to overestimate the ATC (indicated by the estimators for small ρ). The LCB crosses the ground truth at ρ = 0.1; this can be viewed as a lower confidence bound for the true confounding level ρ = 0.125 (we elaborate on this in the discussion when the ground truth is zero). Finally, the LCB hits zero at ρ 0 = 0.65; we can thus conclude with 0.95-confidence that ATC is non-negative as long as the true confounding level does not exceed ρ 0 .
The estimators for bounds of counterfactuals over N = 500 runs for each ρ are plotted in Figure 4 (they are evenly spaced on the x-axis). The simulation results show the sharpness and acuracy of our estimators: they are quite close to the ground truth, especially for small values of ρ; they get a bit conservative and have a larger variance when ρ is as large as 1. Interestingly, there are also a few outliers when ρ is very small, and the estimators seem to be the most stable for an medium scale of ρ (around 0.18 to 0.5). The actual value of E[Y (1) | T = 0] in our design, represented by the red triangles, are very close to the lower solid line, the ground truth of µ − 1,0 ; this means our simulation design is close to the worst case. To further validate our inference procedure, we compute the empirical coverage of CI lower and CI upper for µ ± 1,0 over N = 500 runs. We also compute the ground truth of E[Y (1) | T = 0] under our design as a baseline, and compute the empirical coverage of CI mean . They are plotted in Figure 5. Our empirical coverage is close to the nominal level 0.95 in almost all settings, showing the validity of our inference procedure.   Our theory shows that even though the ERM is off, these C.I.s still have valid asymptotic coverage; such robustness is also supported by empirical evidence.

Discussion
In this work, we propose a new sensitivity model based on the f -divergence that characterizes the average effect of confounders on selection bias. Under the f -sensitivity model, we offer a scheme for the estimation and inference on the counterfactual and the ATE. We close the paper by a discussion on possible extensions.
Tightness. As mentioned before, the optimal value of (6) is not necessarily the tightest lower bound for E[Y (1) | T = 0]: the sharp one under (f, ρ)-selection condition is given by where Q 1,0 is the identification set of all distributions that agree with the observed distribution and satisfy the (f, ρ)-selection condition. The constraints in (6) define a superset of Q 1,0 , potentially leading to conservativeness.
Using the exact characterization of Q 1,0 provided in Proposition 5, we can represent the sharp lower (resp. upper) bound of E[Y (1) | T = 0] under the (f, ρ)-selection condition as the optimal value of min(resp. max) With the same argument, we can also develop the optimization problems for sharp bounds on the ATT and the ATC. Compared to the dual problems in Proposition 3, the additional constraints in the last line above leads to a dual form that is not as clean. Developing an efficient algorithm that solves this tight bound remains an interesting avenue for future research.
Sensitivity analysis. In this paper, we have focused on conducting inference on the counterfactuals and treatment effects under the (f, ρ)-selection condition, with a prescribed confounding parameter ρ. Based on this, we can make robust causal conclusions and conduct sensitivity analysis by inverting the confidence intervals as follows. Suppose the goal is to detect if there is a nonzero ATE; we can consider a increasing sequence of ρ, and construct a level 1 − α confidence interval C(ρ) for the ATE using the method introduced in this paper at each value of ρ; finally let ρ be the smallest ρ such that C(ρ) contains zero. We can interprete the results as either there is a nonzero ATE, or there is a confounder as large as ρ to explain away the observed treatment effects.
More rigorously, let ρ * denote the true confounding level and suppose the constructed confidence intervals C(ρ) are nested in ρ: for any ρ 1 ≤ ρ 2 , C(ρ 1 ) ⊂ C(ρ 2 ). We then have if C(ρ * ) is an asymptotically valid confidence interval for the ATE. In words, when the ATE is indeed zero, ρ is an asymptotic level-(1−α) confidence lower bound for ρ * . Similar to the case of , here only point-wise validity is necessary, i.e., we only need our CIs to be asymptotically valid for each fixed ground truth of ρ. Finally, we note that the monotonicity of the confidence intervals is satisfied with a reasonable estimation procedure; one can also force the confidence intervals to be monotone by enlarging some of them to conform to those for smaller values of ρ, without hurting the asymptotic validity.
Implications for the conditional average treatment effect (CATE). The mothodology proposed in this paper also provides bounds on CATE under the (f, ρ)-selection condition. For example, the proof of Propositions 2 and 3 implies that a lower bound for E[Y (1) | X = x, T = 1] is given by the optimal value of The dual form of the above optimization problem is Note that the optimizer (α * (x), η * (x)) defined in (10) is exactly the optimizer of (11). In fact, µ(x) := r (j) (x) h (j) (x) where r (j) , h (j) are defined in Algorithm 1 is an estimator for the optimal objective in (11).
These quantities are repeatedly estimated on distinct folds of data as intermediate steps of our procedure. While such sample splitting does not compromise the efficiency of inference due to the final averaging step, how to efficiently estimate these CATE bound functions with statistical guarantees might call for distinct considerations from ours. We leave this for future investigation.
Marginal (f, ρ) selection condition. We might even relax the per-x uniform bound on the f -divergence in Definition 1 to a marginal fashion, so that the selection bias is controlled averaged over both U and X. More formally, we might consider the constraint that In this setting, the odds ratio can be very large for a small proportion of X | T = 1, but still controlled in the average sense. This type of marginal (f, ρ)-selection model leads to a larger class of distributional shifts than the (f, ρ)-selection condition here, and a different optimization problem for bounds on counterfactual means. Following similar arguments here, we see that the dual formulation, parallel to Proposition 3, can still be viewed as a risk minimization problem; however, the risk function would involve the unknown X-shift r 1,0 , which might make the estimation and inference more complicated. The estimation and inference under this marginal f -sensitivity model is an ongoing work.

A.1 Proof of Proposition 1
Here we prove a stronger result that directly implies Proposition 1. The following proposition is a tight characterization of the identification set induced by the f -sensitivity models.
Proposition 5. Let (X, Y (0), Y (1), U, T ) ∼ P sup be the true unknown super-population over all random variables of interest. Let P be the set of all distributions over (X, Y (0), Y (1), U, T ). Let P obs X,Y,T be the joint distribution of all observable random variables (X, Y, T ). Let t ∈ {0, 1}. Define Q t,1−t as the set of all counterfactual distributions that agrees with the observables and satisfies the (f, ρ) selection condition, i.e., Q t,1−t = P X,Y (t) | T =1−t : P ∈ P, P X,Y,T = P obs X,Y,T , P satisfies Definition 1 .
Proof of Proposition 5. Fix t = 1. For any P X,Y (1) | T =0 ∈ Q 1,0 , since P X,Y,T = P obs X,Y,T , On the other hand, where the last inequality is due to the (f, ρ)-selection condition. Combining the above, we establish the "⊂" direction. It remains to prove the reverse. We show the proof for the case of t = 1 here, and the t = 0 case follows from similar arguments. Given any Q ∈ Q 1,0 , we aim to find a distribution P sup over (X, • P sup is compatible with P obs X,T,Y ; • P sup satisfies the (f, ρ)-selection condition; To construct P sup , we first set P sup X,T = P obs X,T . Then we specify the distribution of Y (1) | X, T via So far the joint distribution of (X, T, Y (1)) has been determined. We let U = Y (1) be the unobserved confounder. Finally, the distribution of Y (0) | X, T, Y (1), U is specified via Having constructed P sup , we proceed the check that it satisfies the conditions. Conditional on X and U , Y (1) becomes deterministic and the distribution of Y (0) only depends on X. Hence (Y (1), Y (0)) |= T | X, U . By construction, it is straightforward to see that P sup is compatible with P obs X,T,Y . For any x, again by the construction of P sup , where the last inequality is due to the definition of Q. On the other hand, where the last equality is because U = Y (1) under P sup . Combing the above, we have Similarly, Therefore, the super population P sup satisfies the (f, ρ)-selection condition. By construction, P sup Y (1) | X,T =0 = Q Y | X . It remains to show that P sup X | T =0 = Q X . For any measurable set A, Since the above holds for any measuable set A, P X | T =0 = Q X .
Finally, switching the role of 1 and 0 completes the proof.

B.1 Proof of Proposition 4
Given X = x, suppose instead α * (x) = 0. We consider the following two cases: where step (a) uses the fact that liminf n→∞ a n + b n ≥ liminf n→∞ a n + liminf n→∞ b n and step (b) follows from Fatou's lemma and the condition that Above, step (a) is due to the fact that f * (x) is bounded when x ≤ 0 and the dominated convergence theorem; step (b) is because f * (x)/x → 0 as x → −∞.
Combining the two cases above, we conclude that η * (x) = −y(x) and the optimal value of the dual problem is −y(x). By the strong duality, the optimal value of the primal objective function is E r 1,0 (X)y(X) | T = 1 . As an implication, there exists a feasible L(x, y) such that E[Y (1)L(X, Y (1)) | T = 1] = E r 1,0 (X)y(X) | T = 1 . Let Q Y | X=x denote the measure induced by L(x, y): This is a valid transformation of measure because L(x, y) is feasible. Then Y (1) = y(X) a.s. under Q Y | X . Consequently, Again since L is feasible, This is a contradiction to the condition. Hence α * (x) > 0.

B.2 Discussions on Assumption 1 on sieve estimation
We provide additional discussion on Assumption 1 for sieve estimators in the context of (X, Y (1)) | T = 1. In particular, we first justify the smoothness of the optimizers when the conditional distributions are sufficiently smooth. We then verify the technical conditions for two choices of f -divergences: KL-divergence and χ 2divergence. Then we discuss some considerations of relaxing the conditions with implementations in practice.
Smoothness of the optimizers. We first provide some justifications for assuming the optimizers are continuously differentiable. By the strong convexity of f , its conjugate f * is continuous, hence without loss of generality we always assume the differentiation and expectation are exchangeable. We also assume the conjugate f * is sufficiently smooth, which is the case for many popular choices of f -divergence. As we discussed in Proposition 4, under mild conditions, the optimizers (α * (x), η * (x)) lies in the interior of R + × R.
The optimizers are thus the solutions to where the right-hand side takes the form for some differentiable or smooth function g decided by f * and its derivative (f * ) . Thus F (a, b, x) is smooth in (a, b) when f * is sufficiently smooth. Now let us assume the conditional distribution P Y (1) | X=x,T =1 is smooth; for example, for some h ∈ X , P Y (1) | X=x+th,T =1 = P Y (1) | X=x,T =1 + t · P h for some measure P h on Y; and similar for higher-order expansions. This is a reasonable assumption if we are willing to assume that the conditional distributions of Y (1) are close for similar covariates. Concretely, such condition holds when Y (1) | X = x, T = 1 is a normal distribution with homoskedastic noise and a smooth mean function, or heteroskedastic noise with a smooth mean function and smooth standard deviation function, etc. When the conditional distributions are smooth in x, the function F (a, b, x) is also smooth in x by the linearity of conditional expectation. Finally, if the derivatives with respect to a, b is always invertible (which is the case under mild conditions for the examples we discuss shortly) and smooth, invoking the Implicit Function Theorem (Rudin et al., 1976), the minimizer can be smooth in x.
KL-divergence. A popular choice for the function f is f (x) = x log x, which leads to the KL-divergence (Kullback and Leibler, 1951). The dual function in this case is f * (y) = e y−1 , and the loss function becomes The conditional expectation is We first look at the strong convexity assumption. The conditional expectation is twice differentiable, with Therefore, a simple calculation shows that as long as Y (1) is not deterministic at (α * (x), β * (x)), the Hessian matrix is non-singular. Also, if the underlying distribution P Y (1) | X=x,T =1 is continuous in x, the above derivatives, hence the eigenvalues of the Hessian matrix is continuous; since X is compact, there exists a positive uniform lower bound for the smallest eigenvalue of the Hessian matrix, leading to strong convexity.
We then consider the continuity condition | (θ, x, y) − (θ * , x, y)| ≤¯ (x, y) θ(x) − θ * (x) 2 for θ(x) − θ * (x) 2 < for some sufficiently small > 0, where θ(x)−θ * (x) 2 is the Euclidean norm, and sup x∈X E[¯ (x, Y (1)) 2 | X = x, T = 1] < M for some constant M > 0. By Taylor expansion, we have where θ(x) lies between θ(x) and θ * (x). We note that ∇ θ is also a smooth function of θ, and the gradient is uniform bounded for θ(x) within a neighborhood of θ * (x) in terms of Euclidean L 2 -norm. In particular, For any (a, b) − θ * (x) 2 ≤ for sufficiently small , we can take¯ (x, y) as the uniform upper bound of the Euclidean norm of the gradient, which has finite second moment if Y (1) is not too heavy-tailed. Finally, the last condition is that there exists a constant C 1 such that is sufficiently small. Similar to arguments in the proof of Theorem 1, sufficiently small θ − θ * L2(P · | T =1 ) implies sufficiently small θ − θ * ∞ for this function class. Therefore, we can consider θ ∈ Λ p c (X ) × Λ p c (X ) such that θ − θ * ∞ is sufficiently small. With a taylor expansion of the conditional expectation of the risk at (α * (x), η * (x)), we have since the gradient is zero, where θ(x) lies between θ(x) and θ * (x). Previous derivations have shown that the Hessian is continuous; also, by the compactness of X and continuity of θ * (x), there is a uniform lower bound > 0 for α * (x). Thus, when θ − θ * ∞ is sufficiently small, the Hessian is also bounded. Again by the compactness of X , this bound can be taken to be uniform for x ∈ X , which leads to the desired condition.
Another popular choice is f (x) = (x − 1) 2 , so that f * (y) = 1 4 ((y + 2) 2 + − 1). The conjugate function is a quadratic function on [−2, ∞) and zero on (−∞, −2], with continuous gradient (f * ) (y) = ( s 2 + 1) + , and second-order derivative (f * ) (y) = 1 2 1{y > −2}; the latter is almost-everywhere (under Lebesgue measure) except y = 2. We now proceed to verify the conditions. The loss function is Assuming Y (1) does not have point measure, the differentiation and expectation are exchangeable, and Also, the gradient is given by which are both zero at (a, b) = (α * (x), η * (x)). The form of the loss function implies that α * (x) > 0 for almost all x; hence there is a uniform lower bound > 0 by the compactness of X . By Cauchy-Schwarz inequality, the Hessian at (a, b) = (α * (x), η * (x)) is positive P( Y (1)+η * (x) −a > −2} = 0 almost surely for some c(x) ∈ R. By the optimality condition, the former is impossible, and the latter is also impossible if Y (1) is not deterministic conditional on X = x. Thus, as long as Y (1) | X = x is not deterministic for almost all x, the Hessian is positive definite for all x ∈ X . By compactness of X and the continuity, we know that the minimial eigenvalue of the Hessian is uniformly lower bounded away from zero, hence the strong convexity follows. The other two conditions are easy to verify in this case: the conjugate function f * is a truncation of a quadratic function. Since truncation is a contraction map, these results hold easily by the uniform boundedness of second-order derivatives. We've thus verified the conditions in Assumption 1 for χ 2 -divergence.

C.1 Proof of Proposition 3
Proof of Proposition 3. We first claim that solving (6) amounts to solving the following problem for each x: To be specific, denoting the optimal objective of (12) as µ(x) and that of (6) as µ − 1,0 , we are to show that To see why it is the case, suppose L * is the optimizer of (6), then it is measurable with respect to X and Y (1) and satisfies the constraints of (6). Then L(x)(·) := L * (x,) is measurable with respect to Y (1), and satisfy the constraints of (12). As a result, On the other hand, suppose L * (x)(·) is measurable with respect to Y (1) and is the minimizer for (12) for P X | T =1 -almost all x. We let L(x, y) = L * (x)(y), so that it is measurable with respect to (X, Y (1)) and satisfy the constraints of (6).
Combining the two directions leads to the equivalence.
In the following, we solve (12) and write E x in place of E[· | X = x, T = 1] for simplicity. Invoking Luenberger (1997, Theorem 8.6.1) to this convex problem, we have where the Slater's condition is satisfied and strong duality holds, and ϕ(α, η, x) = inf L≥0 measurable L(α, η, L, x), The minimum of L(α, η, L, x) is thus given by Now we write α(x) and η(x) to emphasize its dependency on x. Therefore, by the equivalence discussed in the beginning, we have With a change-of-variable from α(x) to α(x)r 1,0 (x), we have The minimum of (6) can thus be written as where for P X | T =1 -almost all x, it holds that Therefore, we complete the proof of Proposition 3.

C.2 Proof of convergence of sieve estimator
Proof of Theorem 1. We analyze the behavior of θ (j) for each fold j. As |I (j) 1 | n, we take the generic notation of θ and sample size n, so that For some fixed b > 0, we denote the sequence δ n := inf δ ∈ (0, 1) : where N ( , Θ n , · L2(P · | T =1 ) ) is the -covering number of Θ n in the L 2 -norm under P · | T =1 . We employ the established convergence results for sieve estimators adapted from Chen (2007, Theorem 3.2) and Yadlowsky et al. (2018, Lemma B.3), stated in Lemma 2.
Similar to the results in Yadlowsky et al. (2018), we have δ n J d n log n n .
We finally bound the approximation error using Θ n . Note that we take Θ n to be truncated at . However, since the population minimizer θ * is uniformly bounded above , since truncation is a contraction map, where the last inequality is a wellestablished result, see, e.g., Timan (2014). We now set J n = (n/ log n) 1/(2p+d) , so that θ − θ * L2(P) = O P ((log n/n) p/(2p+d) ). This completes our proof.
It thus remains to show (13). We prove the results when either r or h is consistent.
Consistent r. We first show the case where r is consistent for r 1,0 , but not necessarily the regression function h. Recall that H(x, y) = ( θ(x), x, y). Note that The first summation can be controlled as (where the expectation is implicitly conditional on other folds except I (j) Invoking Lemma 3, we can drop the conditioning and the first summation is o P (1). On the other hand, since the covariate shift between P X | T =1 and P X | T =0 is exactly r 1,0 , we know that E[ h(X)r(X) | T = 1] = E[ h(X) | T = 0], where we still implicitly condition on other folds. As a result, where both terms are unbiased. Thus by Cauchy-Schwarz inequality, invoking the assumption that h as finite second moment. Drop the conditioning by Lemma 3, we know that this summation is also o P (1), hence Finally, once θ − θ ∞,P X | T =1 = o P (1), by the local expansion around θ (x), we have hence (implicitly conditioning on other folds) we have since E[M (X, Y (1)) 2 | T = 1] ≤ M for some constant M > 0. We've thus completed the proof of (13).
Consistent h. We then show the results when h is consistent, but not necessarily r. In this case, h − h L2(P X | T =1 ) = o P (1), whereh = E[ H(X, Y (1)) | X = x, T = 1] viewing H as fixed. Note that The first summation is unbiased conditional on other folds, hence Var r(X){ H(X, Y (1)) −h(X)} T = 1 = o P (1) due to the finite second moments. By Cauchy-Schwarz inequality, the second summation satisfies due to the boundedness of r. Similarly, we know that 1 |I0| i∈I0h (X i ) − h(X i ) = o P (1). Consequently, invoking Lemma 3 we drop the implicit conditioning and arrive at µ = 1 |I 0 | i∈I0h (X i ) + o P (1) = E h (X) | T = 0] + o P (1) further using the finite second moment ofh (or H), where we implicitly condition on other folds and view h as fixed. Finally, denoting h (x) = E[ (θ (x), x, Y (1)) | X = x, T = 1], we note that by Jensen's inequality, By the same argument as the previous case and due to the uniform boundedness of the covariate shift r 1,0 (·), the above term is o P (1). Therefore, by the change-of-measure with r 1,0 , we have µ = E h (X) | T = 0] + o P (1) = E r(X)h (X) | T = 1] + o P (1) = E r(X) (θ (X), X, Y (1)) T = 1 + o P (1) by the tower property of conditional expectations. We thus complete the proof of two cases and conclude the proof of Theorem 2.

C.5 Proof of Theorem 4
Proof of Theorem 4. The proof follows exactly the same arguments as the proof of Theorem 3 with θ in place of θ * , where all the errors are controlled in the same way; the only difference is to show that E r(X)∇ θ θ (X), X, Y (1) [ θ(X) − θ (X)] T = 1 in parallel with (15) in the proof of Theorem 3. We note that this is directly implied by our local condition. Therefore, under the conditions of Theorem 4, The terms in the summation has expectation µ 1,0 := −E r 1,0 (X i )H (X i , Y i (1)) T = 1 .

D Technical lemmas
Lemma 3. Let F n be a sequence of σ-algebra, and let A n ≥ 0 be a sequence of nonnegative random variables. If E[A n | F n ] = o P (1), then A n = o P (1).
Proof of Lemma 3. By Markov's inequality, for any > 0, we have → 0 as i → ∞. By the dominated convergence theorem, we have E[B n k i ] → 0, or equivalently, P(A n k i > ) → 0. Therefore, for any subsequence {n k } k≥1 of N, there exists a subsequence {n ki } i≥1 ⊂ {n k } k≥1 such that A n k i P →0 as i → ∞. By the arbitrariness of {n k } k≥1 , we know A n P →0 as n → ∞, which completes the proof.