Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
69,455 characters · 11 sections · 77 citation commands
Partial Identification of Causal Effects for Endogenous Continuous Treatments
Causal inference in observational studies commonly proceeds under the assumption of no unmeasured confounding (NUC), also known as no endogeneity: treatment assignment is effectively randomized within strata defined by a sufficiently rich set of observed covariates. Although this assumption enables identification of causal parameters and effective inference through propensity score matching rosenbaum2010design, inverse-probability weighting horvitz1952generalization and doubly robust and double-machine learning estimators robins1994estimation, van2006targeted, chernozhukov2018double — it is fundamentally untestable from the observed data without an alternative assumption and may be violated in practice.
Lack of sufficient control for confounding can lead to different conclusions from observational and experimental studies, and hence can threaten the validity of conclusions obtained from the former. For example, observational studies had raised concerns about the possible harm associated with the use of calcium channel blockers like nifedipine (psaty1995risk, pahor2000health) with respect to myocardial infarction (commonly known as heart attack), only to be settled definitively about a decade later via randomized controlled trials that nifedipine was safe (psaty2004contemplating). In-depth analysis revealed that calcium channel blockers were used to treat hypertension, which in itself was a risk factor associated with heart attacks (see rutter2007identifying for more details).
To assess the robustness of causal conclusions about violations of NUC, researchers increasingly rely on sensitivity analysis to quantify the impact of potential unmeasured confounding. In response to concerns about unmeasured confounding dating back to fisher1958cigarettes about the impact of smoking on lung cancer, cornfield1959smoking conducted the first known formal sensitivity analysis to establish that an unmeasured confounder would need to be nine times more prevalent among smokers than non-smokers to explain the observed association, a magnitude deemed implausible. Since then, there have been immense interest in developing methods for sensitivity analysis, the works of robins1999association, rosenbaum2010design, vanderweele2011bias, imbens2015causal, chernozhukov2021omitted, kallus2019interval, zhao2019sensitivity, bekerman2024planning, nabi2024semiparametric are some notable works in this regard.
However, most existing works have primarily focused on binary exposures. Continuous exposures provide unique challenges to causal inference problems, since (i) the probability of observing any $\{T=t\}$ is zero, and (ii) there would be infinitely many potential outcome means $\{\mathbb{E}[Y(t)]\}_{t\in\mathcal{T}}$ to identify. Existing literature on causal inference with continuous exposure have been considered in lu2001matching, fogarty2021biased, zhang2023statistical etc. for matching methods, and kennedy2017non, colangelo19052025, schindl2024incremental etc. for double robust estimation, but these works have not included sensitivity analysis to NUC. CATE estimation with continuous covariates also deploys similar techniques as continuous exposures, as illustrated in kennedy2023towards and semenova2021debiased.
Sensitivity analysis with continuous exposures in matched studies has recently been considered in zhang2024sensitivity and frazier2024bias, and interestingly, the first work establishes hardness results for sensitivity analysis in matched studies with continuous exposure and outcome. A more general class of $f$-sensitivity models have also been considered in jin2022sensitivity and frauen2024neural, where they bound a more general convex-link $f$ of the odd's ratio, instead of the identity link, to accommodate the unbounded support of unmeasured confounders. bonvini2022sensitivity also studied sensitivity analysis in marginal structural models with continuous exposures, introducing approaches based on propensity scores (see also jesson2022scalable, marmarelis2023partial), outcome regression, and subset confounding that differ from standard formulations; notably, they bound density ratios in the former model rather than odds ratios (see Section (ref) and Appendix (ref) for further details), and leverage low-dimensional structure to achieve $\sqrt{n}$-rate inference. In this work, we deviate from their strategy and instead consider double-robust partial identification of the dose-response curve, and we directly estimate corresponding bounds under less restrictive non-parametric regression conditions.
Although the need for sensitivity analysis is widely accepted, there remains a broad range of sensitivity models to choose from. We focus on continuous exposure extensions of two well-established canonical binary exposure models: Rosenbaum's model (rosenbaum2002observational), and the marginal model of Tan (tan2006distributional). Both models, originally devised for binary exposures, bound the odds ratio of receiving exposure for different levels of an unmeasured confounder $U$. Moreover, both models include the NUC data-distribution within permissible distributions for any non-trivial value of their sensitivity parameter. We extend these models to the continuous exposure case, and derive novel identification bounds for the dose-response curve under both models. Our work builds upon the existing extensions to multi-level exposures of rosenbaum1989sensitivity, however, we develop our framework for partial inference to respect key inequality constraints anchoring our sensitivity analyses. It may be worthwhile to mention that both aforementioned models belong to a class of $L^\infty$ models, which bounds the maximum deviation of the data-distribution from NUC. Recent works have also explored $L^2$ class of models, which arguably being less interpretable than $L^\infty$ models (zhang2024linftyl2sensitivityanalysiscausal, huang2025variance), produce narrower bounds for the potential outcome mean, and could hence be appealing to researchers. Nevertheless, in this work, we focus on the specified $L^\infty$ models in the case for continuous exposures and defer the development of analogous bounds under $L^2$ or more general models to future work.
We delineate our contributions in this paper as the following:
The remainder of the paper is organized as follows. Section (ref) we introduce the sensitivity analysis setup, and formally introduce the notion of sensitivity functions and the Rosenbaum and Marginal sensitivity models. Section (ref) produces bounds on the dose-response curve using a common regression formulation. Section (ref) derives the double robust mapping using semi-parametric techniques, while Section (ref) establishes the $L^2$-consistency and asymptotic normality in the second-step of counterfactual regression. In Section (ref), we dive into a comparative discussion between the two models and identify contexts we believe are appropriate for each. Finally in Section (ref), we evaluate the performance of our methods through simulations, and apply our method to measuring the effect of second hand lead exposure on blood lead levels in children.
Glossary. We denote the CDF associated with any random variable $A$ as $F_A$, a conditional CDF of $A|B=b$ by $F_{A|B}(\cdot|b)$, and the conditional quantile function $Q_{A|B=b}(\tau) = \inf\{x\in \mathbb{R}:F_{A|B}(x|b)\ge \tau\}$ The density of a random variable $A$ is denoted by $f_A$, and the analogous notation holds for joint and conditional densities. We also define $a_+ =\max\{a,0\}$, $a_-=-\min\{a,0\}$, and $a_\pm^2 = (a_\pm)^2$, for any real number $a$. For any vector $a\in \mathbb{R}^d$, let $\|A\| = (\sum_{i=1}^d a_{i}^2)^\frac 12$ denote its vector norm. Also, for any random variable $A$ (either scalar or a vector in $\mathbb{R}^d$) let $\|A\|_{p} = (\mathbb{E} \|A\|^p)^{\frac 1p}$, for $p = \mathbb{R}^{\ge 0}\cup\{\infty\}$. We shall generally represent nuisance functions by $h$ and their estimates by $\hat h$, and for any $p$ as above, define $\|A\|_{p|\hat h} = (\mathbb{E}[\|A\|^p|\hat h])^{\frac 1p}$. Also, for any matrix $B\in \mathbb{R}^d\times \mathbb{R}^d$, let $\|B\|$ denote the operator norm, ie, the largest eigen value of $B$. Finally, we write two sequences $\{a_k\}_1^\infty$ and $\{b_k\}_1^\infty$ to be $a\lesssim b$ if $\exists c>0$ and $K>0$ such that $a_k\le Cb_k$ for $k\ge K$, and $a\simeq b$ if $a\lesssim b$ and $b\lesssim a$; and for random variables $A$ and $B$, we denote $A\lesssim_P B$ if $A = O_p(B)$, and $A\simeq B$ if $A \lesssim_P B$ and $B\lesssim_P A$.
We assume access to $n$ independent and identically distributed samples of $(W_1,\cdots, W_n)$, where $W = (X, T, Y)$, with $X\in \mathcal{X}$ denoting a rich set of measured confounders, $T\in \mathcal{T}$ is a continuous exposure (or treatment, used interchangeably), and $Y\in \mathcal{Y}$ is the outcome of interest. We adopt the potential outcome framework to delineate our causal quantities, denoting $Y(t)$ as the potential outcome that would have been observed had $T$ been set to $t$ (neyman1923application,rubin1974estimating) Throughout this paper, we are interested in the counterfactual quantity $\{\mathbb{E}[Y(t)]\}_{t\in\mathcal{T}}$, referred to as the average dose-response curve.
Suppose $T|X$ admits a density with respect to some base measure $\lambda_b$. We make two standard causal assumptions of consistency and positivity.
The consistency assumption bridges the observed outcomes $W$ to the counterfactual quantity of interest $Y(\cdot)$, and implicitly assumes away complex situations like interference and spillover effects across units, and multiple versions of the same treatment. While the no-unmeasured confounding (NUC) assumption, alluded to in Section (ref) and formalized as $\{Y(t)\}_{t\in\mathcal{T}}\perp \!\!\! \perp T|X$, in conjunction with Assumption (ref) allows us to identify the potential outcome mean, NUC is generally not plausible in observational studies. Instead, one may imagine an unmeasured confounder $U\in\mathcal{U}$ such that, one may generalize the confounding condition, as in Assumption (ref).
This is a generalization of the standard no-unmeasured confounding assumption (also known as treatment exchangeability or ignorability in the causal inference literature), whereby, while conditioning on $X$ alone does not ensure exchangeability, further conditioning on $U$ does. Assumption (ref) is quite mild, since $U$ is not measured; however, the counterfactual dose response curve is no longer identifiable under the condition, without a different condition.
To address unmeasured confounding by $U$, one strategy is to impose additional structure on their distribution, possibly involving auxiliary variables have been measured, such as instrumental variables or proxies, which may, under stringent conditions, lead to point identification of the counterfactual dose response curve (heckman2010comparing, newey2003instrumental, miao2018identifying). In contrast, we adopt a different approach: we aim to partially identify the potential outcome mean under interpretable assumptions about the relationship between the unmeasured confounders and the exposure. A common strategy is to specify a sensitivity model, following the conceptual framework introduced by Rosenbaum (rosenbaum2002observational). In this work, we consider a generalized version of such a model that accommodates continuous exposures.
Assumption (ref) bounds the ratio of conditional densities of $T$ given $X$ and two distinct values of $U$ by a prespecified function of $t$. Note that taking $\Gamma_{t,t'} = 1$ is equivalent to no unmeasured confounding, since the conditional density of $T|U=u,X$ would no longer depend on $u$. In case of binary $T\in\{0,1\}$, replacing the conditional density of $T$ by its natural conditional probability analogue yields the classical Rosenbaum sensitivity model on the odds ratio scale, with $\Gamma_{1,0} = \Gamma \ge 1$. In fact, the model imposed by Assumption (ref) can be considered a generalization of the sensitivity model with many-ordered treatments introduced in rosenbaum1989sensitivity. It is worth emphasizing that, although we follow the binary-treatment literature in omitting $X$ from the sensitivity function $\Gamma_{t,t'}$, our theory and methods extend directly to sensitivity functions that also depend on $X$, without requiring any additional technical assumptions.
Another relatively more recent sensitivity model is the Marginal sensitivity model of tan2006distributional, also see zhao2019sensitivity, dorn2023sharp, dorn2024doubly, introduced in the context of binary exposures. A natural generalization for continuous exposures entails:
Similar to Assumption (ref), the Marginal Sensitivity Model in Assumption (ref) also bounds the ratio of density ratios, but while the Rosenbaum model compares the density ratios between different strata of the unmeasured confounder $U$, the Marginal model on the other hands compares the full-data conditional density ratio to its marginal observed data counterpart . We provide a detailed discussion of the merits and possible use-cases of each framework in Section (ref).
It is important to note that, unlike usual sensitivity models, both Assumption (ref) and Assumption (ref) are not parametrized by a scalar sensitivity parameter, but rather a general sensitivity function. Replacing a scalar sensitivity parameter with a sensitivity function as ours allows several standard parametric and semi-parametric models to be considered -- for instance, if $U$ has a compact support and $(A,Y)|X,U$ has a multivariate normal distribution, then no scalar parameter exists that satisfies the bounds of either Equation (ref) or (ref). One could in principle potentially produce scalar bounds by truncating the exposure level, but this approach may not be appropriate for unbounded exposure, and even if appropriate, the approach could yield unnecessarily wide confidence intervals, by replacing a naturally varying sensitivity function with a uniform worst-case bound. However to make our models useful, we need to obtain a rich class of sensitivity functions well-suited to our goal. In Proposition (ref) we achieve this goal, and we describe a rich class of admissible sensitivity functions which are provably compatible with Assumptions (ref) and (ref), from which an analyst can draw good candidates for conducting a given sensitivity analysis.
The proof of Proposition (ref), as well as various potential options of functions $\Upsilon_t$ to generate $\Gamma_{t,t'}$ is discussed in Appendix (ref) and Table (ref).
Although Assumption (ref) effectively restricts the degree of unmeasured confounding bias, it does not in of itself provide point-identification of the potential outcome mean. The key identification challenge stems from the fact that: $$\mathbb{E}[Y(t)] = \mathbb{E}_X\left[\int\mathbb{E}[Y(t)|T=t',X] \,dF_{T'|X}(t'|X)\right],$$ but since $Y(t)$ is not observed unless $t'=t$, the above expression cannot be identified under Assumptions (ref) and (ref). yadlowsky2022bounds and dorn2023sharp considered the Rosenbaum model and the Marginal model respectively for binary treatments, and derived partial identification strategies for the potential outcome mean. In this section, we extend their results to the case of continuous exposures. To simplify the exposition, our theoretical developments primarily focus on identifying lower bounds for $\mathbb{E}[Y(t)]$ under the Rosenbaum and Marginal models, respectively; corresponding upper bounds can be be obtained by analogy, upon redefining $Z(t) := -Y(t)$ and obtaining the negative value of the infimum of $\mathbb{E}[Z(t)]$ over all counterfactual distributions under the respective sensitivity models.
Let $\theta_{t'}(x,t) := \inf\{\mathbb{E}_Q[Y(t)|T=t',X=x]: Q\in \mathcal{R}_x\}$, where $\mathcal{R}_x$ is the set of all laws of $(\{Y(t)\}_{t\in\mathcal{T}}, T, U)$ conditional on $X=x$ satisfying Assumptions (ref), (ref) and the Rosenbaum model (ref). Theorem (ref) provides a formal characterization of $\theta_{t'}(x,t)$ in terms of the observed data on which our empirical approach shall be based.
Theorem (ref) produces partial identification of $\inf_{Q\in\mathcal{R}_x}\mathbb{E}[Y(t)]$ by marginalizing over $\theta_{T}(X,t)$, which is identified by the conditional expectile of the distribution $Y|T=t,X$ as a function of $\Gamma_{t,t'}$. Interestingly, expectiles have recently emerged in the literature as an alternative to quantiles (philipps2022interpreting), and the result in Theorem (ref), analogous to the case of binary treatments in yadlowsky2022bounds, underscores the importance of expectiles as a critical causal quantity under sensitivity models. The theorem, proved in Appendix (ref), is the consequence of a key restriction on the likelihood ratio of $L_{t,t'}(y,x) = \frac{d\mathbb{P}(Y(t)|T=t',X)}{d\mathbb{P}(Y(t)|T=t,X)}$, given by $L_{t,t'}(y,x)\le \Gamma_{t,t'}L_{t,t'}(\tilde y,x)$ for all $y,\tilde y$, which is imposed by the Rosenbaum Model (ref).
Analogously, the Marginal model also imposes a restriction on the likelihood ratio, given by $L_{t,t'}(y,x) \in [\Lambda_{t,t'}^{-1},\Lambda_{t,t'}]$. Let us denote $\zeta_{t'}(x,t) := \inf\{\mathbb{E}_Q[Y(t)|T=t',X=x]: Q\in \mathcal{M}_x\}$, where $\mathcal{M}_x$ is the set of all full data laws conditional on $X=x$ that satisfy Assumptions (ref), (ref) and the Marginal model (ref). Then the aforesaid restriction imposed by the marginal model provides an identification of $\zeta_{t'}(x,t)$, as illustrate in Theorem (ref) and proved in Appendix (ref).
Thus if the above condition holds, Equation (ref) could be used to estimate $\zeta$ as a linear combination of the conditional mean, and the truncated tail mean function $\mathbb{E}[Y|Y>q_{t'}(X,T),T,X]$ function, sometimes referred to as the Conditional Value-at-Risk (CVar) function (rockafellar2000optimization). It is interesting to note that Theorem (ref) recovers the adversarial regression formulation of dorn2024doubly in the case of continuous exposure, despite our very distinct proof strategies. Note that both Theorem (ref) and (ref) use outcome regression formulations to identify the estimand of interest, which is a convenient strategy for continuous exposures. An adversarial propensity score construction strategy, like that in dorn2023sharp, is difficult to apply in the continuous exposure setting, due to the appearance of Dirac-delta at $t$ terms when attempting to partially identify $\mathbb{E}[Y(t)]$. Some recent work has considered kernel based localization to circumvent this issue (eg: colangelo19052025) in the context of NUC; but we defer exploring such an approach in our setting to future work. Although we emphasize that while the identification strategies in Theorems (ref) and (ref) bypass the need to incorporate propensity scores (ie, conditional density in this case), they still play a key-role in efficiently estimating the estimands defining our bounds, as highlighted in the following sections.
Let $r(t) = \mathbb{E}[\theta_T(X,t)]$ and $m(t) = \mathbb{E}[\zeta_T(X,t)]$ be the lower bounds on $\mathbb{E}[Y(t)]$ under the Rosenbaum and Marginal models respectively. If $r(t)$ and $m(t)$ were parametrized by finite-dimensional structures, standard tools from semi-parametric efficiency theory (van2000asymptotic, tsiatis2006semiparametric, hines2022demystifying) could be used to derive their efficient influence functions (EIFs). However, under milder smoothness conditions, these functionals are no longer pathwise-differentiable if $T$ is absolutely continuous with respect to the Lebesgue measure, and no $\sqrt{n}$-consistent estimators can exist (bickel1993efficient, diaz2013targeted).
To address this challenge, one can adopt a pseudo-outcome approach analogous to kennedy2017non,yang2023forster, chernozhukov2024conditional. We define $\psi_R = \mathbb{E}[r(T)]$ and $\psi_M = \mathbb{E}[m(T)]$, where the outer expectations are taken with respect to the marginal distribution $F_T$ of the observed exposure. Unlike $r(t)$ and $m(t)$, the population level functionals $\psi_R$ and $\psi_M$ are now both pathwise-differentiable under standard regularity conditions and therefore admit influence functions under an appropriately defined semiparametric model, which can in turn be used to define so-called pseudo-outcomes as corresponding un-centered influence functions. The key insight of pseudo-outcomes is that under conditions of Theorem 2 of yang2023forster, the influence function-based estimators of $\psi_R$ and $\psi_M$ inherit a conditional second-order bias property even for their conditional counterparts, such as $r(t)$ and $m(t)$, thereby providing improved estimators for such nonregular, i.e. non $\sqrt{n}$-estimable, functionals.
With the EIFs for $\psi_R$ and $\psi_M$ at hand, we construct pseudo-outcomes based on the recipe from yang2023forster and dalal2024anytime, as the estimated uncentered EIF, which we fit to produce counterfactual regression estimates of $\hat r(t)$ and $\hat m(t)$ by deploying the following two-step procedure:
Step 1: Divide the dataset into three non-overlapping parts $I_1, I_2$ and $I_3$. From $I_3$, we estimate the relevant nuisance parameters $\hat\theta$, $\hat\nu$, $\hat f(t|x)$ for $\hat r(t)$, and $\hat q$, $\hat \alpha$ and $\hat f(t|x)$ for $\hat m(t)$.
Step 2: For $i\in I_1$, construct the following pseudo-outcomes
where the integrals are computed numerically.
Step 3: For $(W,T)\in I_1$, regress $\hat Y_r(W) \sim T$ and $\hat Y_m(W)\sim T$ using non-parametric ordinary least squares (OLS) on a user-specified basis system, to obtain $\hat r(t)$ and $\hat m(t)$.
We elaborate further on the non-parametric OLS in Section (ref). It is worth mentioning that cross-fitting rather than sample-splitting could be used to improve efficiency. Specifically, one can divide the dataset into $K\ge 3$ non-overlapping parts $J_1,\cdots, J_K$, and setting $I_2 = J_{i_2}$, $I_3 = J_{i_3}$ and $I_1 = \{1,\cdots, n\} - I_2 - I_3$, where $\{i_1,i_2,\cdots, i_K\}$ is a permutation of $\{1,\cdots, K\}$. The ensuing $\binom K2$ estimates could then be averaged to obtain analogous theoretical guarantees on the averaged estimators with full-sample efficiency.
In this section, we elaborate on the non-parametric OLS regression of Step 2 in Section (ref), and provide theoretical guarantees for the convergence of $\hat r(t)$ and $\hat m(t)$ to their respective population analogues.
Let $\lambda$ denote the Lebesgue measure, and let $\Psi = \{\phi_1(\cdot) \equiv 1, \phi_2(\cdot),\phi_3(\cdot),\cdots \}$ be a sequence of functions such that linear combination of these functions are dense in $L_2(\lambda)$. Examples of such sequences include polynomial series, Fourier series, regression splines (huang2003asymptotics), local polynomial partition series (cattaneo2013optimal), wavelet series, etc. The idea is to use these basis functions to approximate any $L^2$ function from the sample at hand. For any $J\ge 1$, let $\bar \phi_J = (\phi_1,\cdots, \phi_J)$ be a collection of the first $J$ basis functions. Also, for any $f\in L^2(\lambda)$, define $E_J^\Psi(f):= \inf_{b\in \mathbb{R}^J}\|f - b^T \bar \phi_J\|_{L^2(\lambda)}$ as the error of approximating $f$ by the first $J$ functions of $\Psi$. The choice of a dense basis ensures that $\mathbb{E}_J^\Psi(f)\to 0$ as $J\to\infty$ for any $f\in L^2(\lambda)$.
Let $Q = \mathbb{E}[\bar\phi_J(T)\bar \phi_J(T)^T]$, and let $\hat Q_{|I_1|} = \frac{1}{|I_1|}\sum_{i\in I_1} \bar \phi_J(T_i)\bar\phi_J(T_i)^T$ be its empirical counterpart in $I_1$. The non-parametric OLS estimator for $a(t)$, $a\in \{r,m\}$ is given by $\hat a_J(t) := \bar \phi_J(t)^T\hat\beta_a,$ where $$\hat\beta_a :=\arg\min_{b\in \mathbb{R}^J} \dfrac{1}{|I_1|}\sum_{i\in I_1} (\hat Y_a(W_i) - b^T\bar \phi_J(T_i))^2 = \hat Q_{|I_1|}\cdot\dfrac{1}{|I_1|}\sum_{i\in I_1} \bar\phi_J(T_i)\hat Y_a(W_i), \ a\in \{r,m\}.$$
To consider the theoretical properties of this estimator, let us define $\xi_J :=\sup_{t\in\mathcal{T}} \| \bar\phi_J(t) \|$, which plays a key role in the choice of $J$. When $\mathcal{T}$ is compact, for polynomial bases, we have $\xi_J \lesssim J$, whereas for bases like the Fourier series, splines, wavelets, local polynomial partitions etc., $\xi_J \lesssim \sqrt{J}$ (belloni2015some). We establish our results under the following technical conditions.
Assumption (ref) (i) ensures that the components of $\bar \phi_J$ are not too co-linear, while maintaining that any linear combination of $\bar \phi_J$ maintains a bounded variance. Such a stability condition can be simply ensured if the components of $\Psi$ are orthonormal on $(\mathcal{T},\lambda)$, and $dF_T/d\lambda$ is bounded above and away from zero (belloni2015some). (ii) is a mild regularity condition required to quantify the $L^2$-error of estimation. (iii) quantifies the rate of growth of $J$ in comparison to the sample size. Condition (i) and (iii) allows $\hat Q_{|I_1|}$ to concentrate around $Q$, and hence obtain necessary rates of consistency. We note that the requirement that $\xi_J^2\log J/|I_1|$ could be relaxed to a lesser degree using potentially different learners than the OLS, for instance the Forster-Warmuth learner considered in yang2023forster. Here, we primarily focus on the OLS regression due to its familiarity and well-established theoretical properties. Condition (iv) allows us to express the $L^2$ rate of convergences in terms of the approximation bias $E_J^\Psi$ with respect to the Lebesgue measure.
Condition (v) guarantees sufficient smoothness of the mapping $\hat\theta_{t'}(X,T)\mapsto \mathbb{E}[\psi^{T,t'}_{\theta_{t'}(X,T)}(Y)|X,T]$ or that of $q_{t'}(X,T)\mapsto \mathbb{E}[\rho^{T,t'}_{q_{t'}(X,T)}(Y)|X,T]$. On closer inspection, the proof of forthcoming Theorems (ref) and (ref) show that Condition (ref) (v) could be relaxed: if $\theta_{t'}(X,T)$ and $q_{t'}(X,T)$, alongwith their estimates, have ranges $\mathcal{A}_r$ and $\mathcal{A}_m$ respectively, then one could instead have ${\text{ess}\sup}_X \sup_{y\in \mathcal{A}_a(X)} f_{Y|T,X}(y|t,x) <\infty,$ which is satisfied when $Y$ is binary and $0<\mathbb{P}(Y = 1|T,X)<1$ for $y\in\{0,1\}$, as the estimates $\hat\theta$ and $\hat q$ will eventually be inside $(0,1)$ and $Y$ would have no mass on any Borel set not containing 0 or 1.
The proof of Theorem (ref) is given in Appendix (ref). Interestingly, one can consider Theorem (ref), in conjunction with Remark (ref) in the case of smooth function classes like $s$-H\"older or $s$-Sobolev spaces (refer to Appendix (ref) for definitions). Suppose $\Psi$ is such that $\xi_J \lesssim \sqrt{J}$ (Fourier, spline, wavelet, local polynomial series, etc.). For $a\in \{r,m\}$, if $a$ belongs to such an $s$-smooth function class, and $\eta_J = O(J^{-2s/d})$, where $d$ is the dimension of $T$, then choosing $J_{|I_1|} \simeq |I_1|^{\frac{d}{2s+d}}$, and assuming the bias due to the estimation of nuisance functions is $\lesssim_P |I_1|^{-\frac{s}{2s+d}}$ yields, $\|\hat a_J-a\|_{2|\hat h} \lesssim_P |I_1|^{-\frac{s}{2s+1}}$, which matches the oracle minimax rate for this problem. The decay condition on $\eta_J$ holds true for H\"older smooth function classes $\Sigma_d^H(s,L)$ and Sobolev spaces of order $s$ for the bases mentioned above. On the other hand, if one uses bases for which $\xi_J$ grows faster than $\sqrt{J}$, for instance, classical polynomial basis where $\xi_J \lesssim J$, then to allow $J_{|I_1|}$ to simultaneously be of the order of $|I_1|^{d/(2s+d)}$ to achieve the minimax rate, while maintaining $\xi_J^2\log\xi_J/|I_1|\to 0$ requires $d/(2s+d)<1/2$ or $s>d/2$. Thus a stricter tradeoff is required when using polynomial bases with non-parametric OLS.
Theorem (ref) also highlights the utility of using pseudo-outcomes constructed in Section (ref). For any basis dense in $L^2(\lambda)$, taking $J\to\infty$ and $J/|I_1|\to 0$ allows achieving consistency, as long as the nuisance estimation errors vanish suitably. In particular, for $\hat r_J$, it suffices that either $\|\hat \theta_T(X,T') - \theta_T(X,T')\|_{4|\hat h}\to 0$, or $\|\hat f(T|X) - f(T|X)\|_{4|\hat h}\to 0$ and $\|\hat \nu_T(X,T') - \nu_T(X,T')\|_{4|\hat h}$. Similarly, for $\hat m_J$, it is sufficient for either $\|\hat q_T(X,T')-q_T(X,T')\|_{4|\hat h}\to 0$, or $\|\hat\zeta_T(X,T') - \zeta_T(X,T')\|\to 0$ and $\|\hat f(T|X) - f(T|X)\|_{4|\hat h} \to 0$. These conditions reflect the lower order bias (more precisely second order bias) built into the pseudo-outcome. In particular, the second order terms appearing in Equations (ref) and (ref) help facilitate the attainment of minimax rates: they allow the estimation errors due to nuisance functions to converge at the rate $ \|I_1\|^{-\frac{s}{2s+d}}$, even when the nuisance functions themselves are less smooth than the target curves $r(t)$ and $m(t)$.
Next, we focus on pointwise limit theory for $\hat r(t)$ and $\hat m(t)$, and establish pointwise asymptotic normality results. For that, define $Y_a(W)$ as the oracle analogue of $\hat Y_a(W)$ in Step 2 of Section (ref); in which, we replace nuisance function estimators $\hat h$ by their true counterparts, and let $\varepsilon_a = Y_a(W) - a(T)$. Then, we formulate Assumption (ref) required to achieve pointwise asymptotic normality.
Assumption (ref)(i) and (ii) are regularity conditions to control the behavior of the target functional $a(t)$ as well as the moments of the nuisance functions. It imposes restrictions on the tails of the regression errors, so that a Lyapunov-type condition can be established to produce normality, even when $J$ varies with sample size. (iii) ensures that the target functionals are non-trivial. (iv) ensures that the pointwise linearization of $\hat a_J(t) - a(t)$ holds under non-parametric OLS, with unknown design matrix $Q$. These conditions lead us to the asymptotic normality result we establish in Theorem (ref).
The proof of Theorem (ref) can be obtained in Appendix (ref). Note that the rate conditions imposed upon the nuisance function estimates are of second order, as was obtained in Theorem (ref), again a consequence of the influence function-based pseudo-outcome construction. The condition $a(t) - \bar\phi_J(t)\beta_J^a = o(\|s(t)\|)$ can be perceived of as an undersmoothing condition. Typically, one would expect $s(t) \simeq \sqrt J$, and in that case the undersmoothing condition can be replaced by $l_J = o(\sqrt{J/|I_1|}).$ If one considers $a$ belonging to a H\"older or Sobolev classes of smoothness $s$, then $l_J \lesssim_P J^{-s/d}\log J$, and thus we require $\sqrt{I_1}J^{-\frac{s}{d} -\frac 12}\log J\to 0$.
Theorem (ref), proved in Appendix (ref), allows for consistent estimation of the asymptotic variance, and thus $\Omega$ in Theorem (ref) can be consistently replaced by $\hat\Omega$ to produce pointwise confidence intervals for $r(t)$ and $m(t)$. Uniform confidence bands may also be obtained using a similar procedure under slightly stronger conditions, analogous to Section 4.3 of belloni2015some.
A comment is warranted about the estimation of nuisance functions. Despite the dependence on both $t$ and $t'$, the estimation of the expectile $\theta_{t'}(x,t)$ and $q_{t'}(x,t)$ can be considered an $M$-estimation problem, as highlighted in Equation (ref) for the former, and an analogous pinball loss-formulation for the latter conditional quantile. Conditions and rates for convergence for such general $M$-estimation problems using sieve estimators have been extensively considered in chen1998sieve, chen1999improved, chen2007large (Chapter 76). Learning expectile regression as in $\theta_{t'}(x,t)$ have also been considered using kernel methods (farooq2019learning) and support vector machines (farooq2017svm). Having obtained $\theta_{t'}(x,t)$ and $q_{t'}(x,t)$, the additional nuisance functions $\nu_{t'}(x,t)$ and CVar could again be formulated as $M$-estimation problems, and sieve-based estimators could in principle likewise be used. Estimation for the CVar function has also previously received some attention, for instance see cai2001weighted. Estimation of conditional density functions has also been well studied in the literature, for instance one could deploy Nadaraya-Watson estimators (see wand1994kernel). In practice, one could also use machine learning to estimate nuisance functions for additional flexibility and potentially improved empirical performance To implement such methods, we propose dividing $I_3$ into a partition $\{I_3^1, I_3^2\},$ where we learn the first stage nuisance functions $\theta_{t'}(x,t)$ and $q_{t'}(x,t)$ using ML techniques, and learn $\nu_{t'}(x,t)$ and $\zeta_{t'}(x,t)$ from $I_3^2$ and estimates $\hat\theta_{t'}(x,t)$ and $\hat\zeta_{t'}(x,t)$. The integrals required to compute the pseudo-outcomes in $\hat Y_r$ and $\hat Y_m$ may be evaluated numerically using Gaussian quadratures (stoer1980introduction).
Both the Rosenbaum and marginal sensitivity models impose restrictions on the ratio of density ratios, which -- unlike the approach in bonvini2022sensitivity -- can be interpreted analogously to odds ratios in the binary treatment setting (chen2007semiparametric, tchetgen2010doubly). In contrast, models that directly bound the likelihood or density ratios by a scalar—such as those in bonvini2022sensitivity and jesson2022scalable—impose overly strong restrictions on the conditional distribution of the exposure, often making the assumption difficult to satisfy in practice. Even when the scalar is replaced by a sensitivity function, as in Assumptions (ref) and (ref), the interpretation remains opaque. For example, if the exposure space $\mathcal{T}$ is truncated— as is common for technical reasons in nonparametric regression—the meaning of the sensitivity parameter becomes entangled with the choice of truncation threshold, complicating interpretation. We elaborate on this issue further in Appendix (ref).
In comparing the Rosenbaum and Marginal models, the key distinction is that the former restricts the density ratios for pairs of $u$ and $u'$, whereas the latter is instead constraining the conditional and unconditional distribution of the exposure $T$ within a strata of $X$. There is an interesting geometric interpretation to these models, as delineated in Proposition (ref).
Proposition (ref), proved in Appendix (ref), highlights the key distinction between the Rosenbaum and Marginal sensitivity models (analogous distinctions still hold in the case of binary exposures). The Rosenbaum model restricts the diameter of the set characterizing the relationship between the exposure and the unmeasured confounder $U$, whereas the Marginal model instead restricts its `radius' from a reference point $u^*$ (depending on the marginal $T|X$). This also allows us to rediscover familiar relationship between the Rosenbaum and marginal models, as pointed out in zhao2019sensitivity (reformulated to accommodate continuous exposures in (ref)), using simple geometric arguments (Appendix (ref)).
The Rosenbaum sensitivity model has another important feature worth emphasizing: it is invariant to the marginal distribution of $U$, since it only specifies a relationship in the treatment mechanism $T|U,X$, and is always compatible with the observed data distribution $(Y,T,X)$, as highlighted in Theorem (ref).
Theorem (ref) is related to a result in osius2009asymptotic, though we provide a fully self-contained proof in Appendix (ref). The theorem underscores a key distinction of the two models - not only is the Rosenbaum sensitivity model data-compatible, a property it shares with its Marginal counterpart, but it also leaves the marginal distribution of $U$ unrestricted, making it particularly suited for addressing certain sensitivity questions.
To elaborate upon that point, one can think of a sensitivity analysis to NUC in the following distinct ways - (a) as a global statement\footnote{leamer1985sensitivity uses the term “global sensitivity analysis" in a related but distinct concept; our interpretation in inspired by, but not identical to, his.}: `How robust are the study's conclusions to any potential unmeasured confounder?', or (b) as a targeted statement: `How would the results of the study vary if a specific, plausible unmeasured confounder would have been included?'. Because the Rosenbaum model is agnostic to the marginal law of $U$, it is well-suited for global statements of type (a), by remaining agnostic to the nature of the confounder and seek a general sense of robustness to causal conclusions. On the other hand, the marginal sensitivity model, as elucidated in Proposition (ref), measures the distances from a particular confounder level $u^*$, which depends on the marginal distribution of $U$ (through that of $T$). Thus it may be more informative when one has external information about the unmeasured confounders, either through knowledge of say the prevalence of a binary $U$ in the population, negative control proxies, instrumental variables, etc. In such cases, the marginal model may yield tighter bounds than the Rosenbaum model (since if one believes the Marginal model with $\sqrt{\Lambda}$, it still needs the Rosenbaum model with parameter $\Lambda$), subject to the validity of the implied marginal distribution of $U$.
The distribution of $U$ could also be informed by other important aspects of a study, like a negative control outcome or an instrumental variable, and hence can thus be incorporated into a marginal model to produce tighter bounds on the estimands of interest. While the Rosenbaum model can still be used in such analyses, the Marginal model with function $\Lambda$ being a subset of the Rosenbaum model is expected to yield tighter bounds, subject to the credibility of the marginal distribution of $U$.
We next demonstrate the performance of Marginal and Rosenbaum sensitivity models in simulated instances and a case study.
We evaluate the finite-sample performance of both sensitivity models using the following data-generating process:
We are therefore evaluating the methods under NUC conditions. The distribution for $T$ is similar to that in kennedy2017non, and the conditional variance of $Y$ varies with both $X_4$ and $T$ incorporate heteroskedasticity. We used sensitivity functions $\Lambda_{t,t'} = \exp(\log(5)|t-t'|)$ and $\Gamma_{t,t'} = \exp(\log(25)|t-t'|)$, in line with the relationship between the two models from Corollary (ref).
Conditional density is estimated via binning with neural networks, while outcome regression is estimated via xgboost. Quantiles in the Marginal model are fit using natural splines in $X$ and polynomial sieve in $T$. CVar function has been estimated by xgboost, and $\hat\zeta_{t'}(x,t)$ is constructed via Equation (ref). For the Rosenbaum model, expectiles are learned across a grid of $\tau$ values, and the closest $\tau$ to $(1+\Gamma_{t,t'})^{-1}$ has been used to approximate $\theta_{t'}(x,t)$. Both expectiles and conditional tail CDF are estimated by $\texttt{xgboost}$, and $\nu_{t'}(x,t)$ is estimated via plugins of the same.
We use a sample size of $n= 5000$, split into three parts: roughly 2500 for nuisance estimation, and the rest split in half for the aggregation ($I_2$) and evaluation ($I_1$) sets. The final regression step uses a degree-5 polynomial regression, and the estimates have been cross-fitted to gain efficiency.
The experiment is repeated 500 times, and Figure (ref) shows boxplots of the resulting dose-response and sensitivity bound estimates across a grid of $T$. The NUC estimates (from semenova2021debiased), as well as the true-curve are overlayed for reference. As expected, the estimated curves generally agree with the truth. However, despite the moderately large sample size, we observe finite-sample smoothing bias in the estimates, which is more pronounced in the bounding curves than in the curve based on NUC. These biases arise from nuisance estimation as well as from the undersmoothing inherent in approximating the curves with polynomial bases. Overall, the simulations highlight the tradeoff between robustness and flexibility in nonparametric sensitivity analysis. While correct parametric or semiparametric specifications could potentially yield faster convergence rates and thus reduce finite sample bias, they could be inconsistent if incorrect. On the other hand, nonparametric estimation of high-dimensional nuisance functions using machine learning methods generally requires larger sample sizes for reliable convergence -- especially when estimating tail behavior.
In this subsection we apply our proposed methods to estimate the effect of second-hand smoke exposure on blood-lead levels in children. While the link between active smoking and elevated blood lead levels is well-established, we focus here on the impact of passive exposure—a critical public health question given the heightened vulnerability of children’s developing nervous systems. Elevated blood lead in children has been linked to reduced intelligence, stunted growth, and anaemia (national1993measuring).
mannino2003second found an association between second-hand smoke exposure and higher blood lead levels using data from National Health and Nutrition Examination Survey (NHANES). However, the observational nature of the data raises potential concerns about unmeasured confounding. zhang2020calibrated addressed this by conducting a sensitivity analysis, finding strong evidence for a causal effect even after accounting for moderate unmeasured confounding. Their analysis used cotinine— a metabolite of nicotine— as a biomarker, dichotomizing its levels to define exposure.
We extend this line of work using more recent NHANES data from 2003–2004 and 2011–2016. Following mannino2003second, we restrict to children aged 4-16, and having cotinine levels at most 15.0 ng/mL to exclude active smokers, resulting in 6616 children to be included in our dataset. Measured confounders include poverty-income ratio, age, sex and number of rooms at home. Our method estimates the effect of relative changes in cotinine on blood lead levels. Plausible unmeasured confounders in our analysis includes school type (more active smokers or older lead-based paint in public school buildings american2005lead), or source of drinking water (gibson2020children).
Figure (ref) shows a steady rise in log blood lead levels with rising cotinine, consistent with prior findings. For lower cotinine levels, the NUC curve is tightly enveloped by the upper and lower bounds from both the Marginal and Rosenbaum sensitivity models, suggesting robustness to even strong unmeasured confounding. At higher cotinine levels, the bounds begin to widen, possibly due to increased heterogeneity in blood lead levels or smaller sample sizes (reflected in wider standard error estimates) in that exposure range.
This work extends the scope of sensitivity analysis in causal inference by introducing generalizations of the Rosenbaum and Marginal sensitivity models to the setting of continuous exposures. Our framework replaces classical scalar sensitivity parameters with sensitivity functions that vary with the exposure level, enabling richer modeling of exposure-confounder relationships and sharper identification bounds for the dose–response curve. In addition, we offer a flexible collection of sensitivity functions—each compatible with the observed data—allowing researchers to tailor their choice based on the contextual demands of the application.
We obtained novel partial identification results for the dose-response curve $\{\mathbb{E}[Y(t)]\}_{t\in\mathcal{T}}$. Our results indicate that both the sensitivity bounds depend on tail behavior of the outcome distribution-- expectiles for the Rosenbaum model and quantiles for the Marginal models. We also establish influence function-based estimators for the sensitivity bounds using a debiased pseudo-outcome formulation, and the resulting estimators are shown to be consistent, achieve the minimax rate, and asymptotically normal under weaker second-order convergence rates of the nuisance functions.
A common theme we observe across both sensitivity models is that the resulting bounds depend on the tail behavior of the conditional outcome distribution—through expectiles, tail CDFs, quantiles, or CVaR. This matches our intuition --- unmeasured confounding can cause the likelihoods of $Y(t)|T,X$ and $Y|T,X$ to differ markedly, with the worst-case departure quantified by the sensitivity model. The corresponding worst-case likelihood ratio therefore relates to the behavior of the outcome at the extremes, making the functionals tail-dependent. This highlights a key challenge in estimating sensitivity bounds --- as the gap between the bounding sensitivity functions widens, the estimand increasingly depends on tail properties of the outcome distribution, making nonparametric estimation more difficult and data-intensive.
We also provide a novel geometric connection between the Rosenbaum and Marginal models, relating to the diameter and radius of a set of density ratios. Moreover, we explore the practical implications of using one model over the other --- the Rosenbaum model is more suited to global sensitivity statements, while the Marginal model provides a targeted approach in settings where one has auxillary information or substantive domain knowledge about the nature of the unmeasured confounder.
Overall, we hope that this work contributes toward a more flexible and interpretable toolkit for sensitivity analysis in the continuous treatment setting, and provides a foundation for future methodological and applied developments in this domain.