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.
63,801 characters · 23 sections · 100 citation commands
A sensitivity analysis for the average derivative effect
Drawing causal inferences from observational studies is challenging for a variety of reasons. Two prominent challenges are a) the treatment (or exposure) of interest is continuous rather than binary or discrete, and b) the treatment/exposure assignment is not random, and so treated and control groups may substantially differ on the basis of unmeasured confounders. We are motivated by the intersection of these two challenges.
In observational studies, continuous exposures are widespread. In epidemiology, researchers are often interested in quantifying the effect of lead exposure or air pollution on health outcomes Dominici2006, Gibson2022. In the social sciences, researchers may be interested in quantifying the effect of household's income on a child's later life educational or economic outcomes Lundberg2023. The vast majority of causal inference methods focus on treatments that are binary or take on a finite number of values. A popular approach when the exposure is continuous is to simply dichotomize the exposure at some threshold. However, this is undesirable because dichotomizing discards potentially useful statistical information from the exposure level and can invalidate the downstream inference VanderWeele2013, or at least muddies its interpretation Lee2024. When the exposure of interest is continuous, some popular causal estimands include the dose-response curve (also known as the exposure response function) kennedy_dose_response, a (projection) parameter of the dose-response curve Bonvini2022SensitivityModels, average derivative effects Newey1993EfficiencyModels, Klyne2023AverageLearning, and stochastic intervention effects Schindl2024, among others.
In this paper, we focus on the average derivative effect (ADE), which has been extensively studied in economics Hardle1989, Newey1993EfficiencyModels, with recent work in statistics clarifying its causal interpretation Rothenhausler2019, Hines2021ParameterisingEffects. It has many other names, and has been referred to as an incremental effect Rothenhausler2019, average partial effect Klyne2023AverageLearning, or an average causal derivative Chernozhukov2022LongLearning. Roughly speaking, the ADE measures the average effect of infinitesimally increasing the exposure level for all units in the population. The ADE corresponds to well-known quantities in economic applications. For example, when one is interested in the effect of the price of a good on demand, the ADE corresponds to the average price elasticity of demand. In cases where one is interested in the effect of increased disposable income on consumption, the ADE corresponds to the average marginal propensity to consume bruns2025two. As elucidated in Hines2021ParameterisingEffects, the ADE estimand has some attractive properties, relative to the widely studied dose-response curve. For example, the ADE is a single number summary of the causal effect; it is difficult to summarize the dose-response curve by a single number. Relatedly, estimation of the ADE, a scalar, is a simpler statistical task than estimating the dose-response curve, which is infinite-dimensional. The ADE also relies on a weaker version of the overlap/positivity assumption than that required for identification of the dose-response curve. Finally, as Hines2021ParameterisingEffects explain, it may be difficult to envision an intervention that sets exposure to the same level for everyone, which is what the dose-response curve measures. Of course, there are settings where the ADE is less appropriate than the dose-response curve. Notably, if the causal effect is not monotonic, an average of positive and negative derivatives could cancel and result in a near zero ADE, rendering the ADE an inadequate summary of the causal effect. Data where the treatment is continuous almost exclusively comes from observational studies, as continuous treatments are rare in randomized experiments. As a result, the ADE is typically only a relevant estimand in observational studies, where unmeasured confounding can never be ruled out. It is therefore of interest to develop methods for ADE estimation that take into account potential unmeasured confounding.
A popular way to alleviate concerns about unmeasured confounding in an observational study is to perform a sensitivity analysis. Dating back to Cornfield2009, a sensitivity analysis acknowledges the existence of unmeasured confounding, but asks how strong it must be to overturn a qualitative causal conclusion. There are now a wide array of sensitivity analysis methods developed, especially for binary exposures. Some sensitivity models constrain the effect of the unmeasured confounders affect on the treatment assignment rosenbaum_obs, Tan2006AScores. Others constrain how far the potential outcome distribution can be from the observed outcome distribution Robins2000, Diaz2013, Nabi2024. Some focus on worst-case departures from the no unmeasured confounding assumption Yadlowsky2022,Dorn2023, and others on average-case departures Huang2024, Zhang2022a. Some are parametric in nature Imbens2003, Cinelli2020, Zhang2022 while others are nonparametric. Researchers have also proposed methods that constrain the proportion of unmeasured confounding Bonvini2022.
For the ADE, we introduce a new sensitivity analysis model that bounds worst-case departures of the treatment assignment from no unmeasured confounding. The model can be thought of as a generalization of the marginal sensitivity model for binary treatments due to Tan2006AScores, and is closely related to the model of rosenbaum1989sensitivity, which itself is a generalization of Rosenbaum's sensitivity model for binary treatments rosenbaum_obs. Under the sensitivity model, for a fixed value of the sensitivity parameter, we derive closed-form upper and lower bounds for the ADE. We consider both binary and continuous outcomes, which lead to different bounds. We introduce nonparametric, robust estimators for the bounds, as well as corresponding confidence intervals. These estimators can leverage data-adaptive nuisance estimators that may converge slower than the parametric rate (though not too slowly). One attractive property of our approach is that the closed-form bounds lend themselves to conducting simultaneous inference over an interval of sensitivity parameters in a particularly simple way.
There is a growing body of literature on sensitivity analyses for non-binary treatments. For the multi-valued treatment case, Basit2023SensitivityTreatments generalize the model and approach of Zhao2019SensitivityBootstrap to derive bounds on linear combinations of potential outcome means. For the continuous treatment case, most sensitivity analysis methods have been aimed at bounds on (aspects of) the dose-response curve/exposure response function. For example, Bonvini2022SensitivityModels propose sensitivity analyses for parameters of marginal structural models for any type of treatment and an array of estimands, including the dose-response curve. Jesson2022, Marmarelis2023 Frauen2023, and Baitairian2024 all propose distinct sensitivity analyses/models for the dose-response curve. One the models considered in dalal2025partial is closely related to ours, but the estimand of interest is the dose-response curve. The recent work of Levis2024 proposes a suite of sensitivity analysis methods for stochastic intervention estimands. Meanwhile, Zhang2024a and Zhang2024b propose sensitivity analysis methods for matched studies with continuous treatments under a related sensitivity model.
Chernozhukov2022LongLearning study omitted variable bias (analogous to sensitivity analysis) within a general class of estimands that are continuous, linear functionals of a conditional expectation. Their results apply to a general class of estimands, with the ADE being one example. They assume bounds on the $L_2$ distance between the “short” and “long” Riesz representers and conditional expectation functions. We take a different approach by appealing to a sensitivity model that imposes a bound on the $L_\infty$ distance between the odds ratio of the observed and unobserved generalized propensity score at any pair of points in the support of the exposure. This gives rise to different interpretations of the sensitivity parameters and also different bounds on the ADE. We give a more detailed comparison to previous sensitivity models in Remarks (ref) and (ref) in Section (ref).
The paper is organized as follows. We formally introduce notation, assumptions, and the causal estimands in Section (ref). In Section (ref), we introduce the sensitivity model and discuss its relation to previous models. In Section (ref), we compute worst-case bounds on the ADE under the sensitivity model. We then propose efficient and robust estimators for the bounds in Section (ref). We then apply the methods in a simulation study (Section (ref)) and two real data applications (Section (ref)).
The underlying data are assumed to be a vector of independent and identically distributed samples from some unknown distribution $(A,X,U,Y(a)_{a \in \mathcal{A}}) \sim F$, where $A$ is a continuous exposure (also referred to as a treatment or dose), $X$ are observed pre-exposure confounders, $U$ some unobserved pre-exposure confounders, $Y(a)$ are potential outcomes, and $\mathcal{A}$ is the support of the exposure. The potential outcomes represent the outcome quantity that would have been observed if the exposure had been externally set to $a$. In the observed data, we only see one of the potential outcomes. Specifically, the observed outcome satisfies $Y = Y(a)$ when $A = a$ (Assumption (ref)). Thus, the observed data consists of $n$ i.i.d. samples of the data vector $(X, A, Y)$. Throughout the paper, we will consider the two scenarios where $Y$ is continuous and $Y$ is binary separately. We will use $'$ and $\partial_a$ to denote taking a partial derivative with respect to exposure $a$, which will be clear from the context. $\norm{\cdot}$ represents the $L_2$ norm, $\perp \!\!\! \perp$ denotes statistical independence, and $\mathbbm{1}$ denotes the indicator function.
We now formally introduce the causal estimand, which is well-defined for both binary and continuous (potential) outcomes. The ADE is defined as
It is useful to unpack the meaning of $\theta$. For some small $\delta > 0$, note that $E[Y(A+\delta)-Y(A)]$ captures an average difference in outcomes in two worlds. In the first world, every unit's observed exposure is increased by $\delta$ from its natural value haneuse2013modified. In the second, there is no intervention on the exposure, i.e. exposure takes its natural value and $E[Y(A)] = E[Y]$. $\theta$ then captures the limiting difference between average outcomes under the small $\delta$ shift intervention and no intervention, scaled by the shift hines2025learning. When $Y$ is continuous, $\theta$ also matches the $E[Y'(A)]$ estimand introduced in Rothenhausler2019 under additional smoothness conditions that ensure $E[Y'(A)]$ is well-defined. The required regularity conditions are made explicit in Appendix (ref). We now introduce causal assumptions that allow the ADE to be written as a statistical functional of the full data, some components of which are not observed.
Consistency requires that the potential outcome corresponding to an exposure $a$ matches the observed outcome when the observed exposure matches $a$. It also prohibits interference, where the exposure of one individual can affect the outcome of another.
Latent ignorability essentially requires that the $(X, U)$ vector contains all common causes of the exposure $A$ and outcome $Y$. This assumption is notably weaker than the typical ignorability or no unmeasured confounding assumption, which requires $\{Y(a)\}_{a \in \mathcal{A}} \perp \!\!\! \perp A \mid X$. Moreover, it is a relatively mild assumption since $U$ is unobserved, and there are no restrictions placed on $U$.
In contrast to a typical overlap assumption that would require $f(a \mid x, u) > 0$ for all $a \in \mathcal{A}$ (this would be needed for the dose-response curve), the local overlap assumption only requires that at any level of confounders $(u, x)$ for which $f(a \mid x, u) > 0$, the conditional density must be positive in a small neighborhood around $a$ as well. We only require this weaker notion of overlap as we focus on the ADE. We will also refer to the conditional density $f(a \mid x, u)$ as the generalized propensity score throughout the remainder of the paper Imbens2000. The next assumption collects additional regularity conditions on the conditional densities and potential outcomes.
These regularity assumptions ensure that certain statistical objects are well-defined. In essence, they place smoothness restrictions on the potential outcomes and conditional densities.
Here, $-s(A\mid X, U)$ is the Riesz representer for the statistical functional $E[\partial_a E[Y \mid A , X , U ]]$ Powell1989. Two statistical examples where the interpretation of $\theta$ is fairly simple are the partially linear model and the single index model.
For the remainder of the paper, we focus on the statistical functional $\theta$ in the nonparametric model. In Appendix (ref), we briefly outline how our results could be extended to weighted average derivative effects, i.e. $\theta_w \equiv E[w(A, X)\partial_a E[Y(A) \mid X, U]]$, for weights $w(A,X)$ that are nonnegative and such that $E[w(A, X)] = 1$. Hines2021ParameterisingEffects and Hines2023OptimallyEffects discuss the interpretation of such weighted average derivative effects.
Based on the previous section, we can write down the causal estimand in terms of the statistical functional $\theta$. However, $\theta$ is a function of $E[Y \mid A, X, U]$, where $U$ is unobserved. Thus, $\theta$ cannot be estimated from observed data without further assumptions. In this section, we introduce a sensitivity analysis model that, through a parameter $\gamma$, limits the strength of association between the unmeasured confounder $U$ and the exposure $A$ conditional on $X$. The sensitivity model then facilitates deriving upper and lower bounds on the estimand of interest based on the allowable amount of unmeasured confounding $\gamma$. Recall $\theta$ can be expressed by $\theta = E[\partial_a E[Y \mid A , X , U ]]$ or as a function of its Riesz representer, $\theta = E[-s(A \mid X, U)Y]$. The sensitivity model we consider is a generalization of the model of Tan2006AScores to the continuous exposure case and is related to Rosenbaum's semiparametric model for continuous doses rosenbaum1989sensitivity. The model restricts the odds ratio of the generalized propensity score (including $u$) vs. the generalized propensity score (marginalizing over $u$) at any two dose levels:
This model generalizes Tan's marginal sensitivity model for binary treatments, as it exactly reduces to the model of Tan2006AScores when $a$ and $a'$ are replaced by 0 and 1, and the $\Gamma$ from Tan2006AScores is set to $\exp(\gamma)$. The Tan model has been studied extensively Zhao2019SensitivityBootstrap, Dorn2023, and has recently been generalized to longitudinal settings bruns2023robust, Tan2025. dalal2025partial consider a more general formulation to bound dose-response curves. The connection between the Tan and Rosenbaum models has been previously discussed when $A$ is binary and continuous Zhao2019SensitivityBootstrap, dalal2025partial. We examine the connection between model (ref) and Rosenbaum's model for continuous doses specific to our context in Appendix (ref). In Appendix (ref), we discuss implications of instead imposing $g(-\gamma(|a-a'|)) \leq \frac{f(a' \mid x, u)f(a\mid x)}{f(a\mid x, u)f(a' \mid x)} \leq g(\gamma(|a-a'|)) \ \forall a, x, u$ for a smooth, nonnegative, strictly increasing function $g$ such that $g(0) = 1$. We now comment on the relationship between model (ref) in relation to other sensitivity models for continuous exposures.
We next establish an implication of the sensitivity model that will facilitate reducing the problem to an optimization problem with a mathematically tractable form.
This result demonstrates that the sensitivity models constrain the Riesz representer $s(a \mid x, u)$ to be within $\gamma$ of $s(a \mid x)$, on the additive scale and in a symmetric fashion. Next, we introduce a statistical restriction (entirely separate from the sensitivity model) on $s(a\mid x,u)$. The restriction is a well-known property of score functions.
Thus, combining the previous two lemmas with the fact that $\theta = E[-s(A \mid X, U)Y]$, we can formulate the sensitivity analysis as an optimization problem with (ref) and (ref) as constraints. As a result, a valid sensitivity analysis under model ((ref)) for the ADE would solve the following optimization problem:
It is straightforward to see that the optimization can be conducted in each stratum $(A = a, X = x)$ separately. Thus, we focus on solving the following formulation:
These optimization problems formulation follow a Lagrangian formulation, and they resemble other optimization problems in the causal inference literature, for example those in Jin2022SensitivityPerspective, Zhang2022a, Dorn2023, among others. Equipped with this formulation of the optimization problem, we will aim to obtain closed-form solutions for the cases where $Y$ is continuous or binary.
We first consider the case where the outcome $Y$ is continuously distributed, i.e. for all $a, x$, the distribution $Y \mid A, X$ has no point masses.
One might observe that the solution has a Neyman-Pearson flavor. This flavor of solution in sensitivity analysis has been observed before, for example for the average treatment effect in the binary treatment case Dorn2023, Zhang2022a. In addition, the optimal values of the optimization problem equal $E[-s(A \mid X)Y]$ (what one would estimate if the unmeasured confounding due to $U$ is ignored), plus or minus $\gamma$ times a nonnegative correction term. Thus, the bounds are symmetric around $E[-s(A \mid X)Y]$. It is clear that the correction term is nonnegative, since it is exactly the average difference between outcomes above and below the conditional median.
The above formulation and closed form solution required the outcome $Y$ to be continuous. This can be seen from the fact that the optimal choices for $s$ depend on $Y$ being above or below some median cutoff point. In the binary case, for strata of $(A,X)$ where $0 < P(Y = 1 \mid A, X) < 1$, one cannot simply take $s^*(A \mid X, U)$ to match the form in Proposition (ref), replacing $M( A, X)$ with 1/2, as this will lead to a violation of the constraint on the score in Equation (ref) unless $P(Y = 1 \mid A, X) = 1/2$ exactly. Of course, the bound obtained by replacing $M(A,X)$ in the solution from Proposition (ref) with 1/2 would still be valid, but potentially conservative. Instead, we can show the following result for the binary outcome case:
Again, the solution has a Neyman-Pearson flavor and the optimal value takes the form of $E[-s(A \mid X)Y]$ (what one would estimate if they ignored $U$), plus or minus $\gamma$ times a nonnegative correction term, which increases as $P(Y = 1 \mid A, X)$ approaches 1/2. In the binary case, however, the correction term is not as simple to estimate. We observe that the optimal value for the binary outcome involves an absolute value (or maximum), because of the term $\gamma E[1/2 - |P(Y = 1 \mid A,X) - 1/2|] = \gamma E[1/2 - \max\{P(Y = 1 \mid A,X) - 1/2, 1/2 - P(Y = 1 \mid A,X) \}] = \gamma E[\min\{1 - P(Y = 1 \mid A,X) , P(Y = 1 \mid A,X) \}]$. Such a term is not smooth when the probability that $P(Y = 1 \mid A,X) = 1/2$ is not zero. An approach that is popular when trying to estimate such non-smooth quantities is to instead target a smooth approximation that bounds the true quantity of interest, which we describe in detail in Section (ref).
For estimation and inference for the closed-form bounds, we appeal to semiparametric efficiency theory Tsiatis2006. As alluded to previously, the bounds for a binary outcome can be non-smooth, so we instead target a smooth approximation. The central object in semiparametric efficiency theory is the efficient influence function, whose variance equals the semiparametric efficiency bound, and is unique in a completely nonparametric model. The rest of this section is devoted to characterizing the efficient influence functions of the (smoothed) bounds for continuous and binary outcomes in the nonparametric model, which will motivate construction of estimators. For convenience, we will at times refer to the efficient influence function even when the precise terminology would be the uncentered efficient influence function. For a recent review of semiparametric theory, we refer the reader to Kennedy2022.
In this subsection, we will derive the efficient influence function for the bounds on the ADE under the sensitivity model with continuous outcomes. Recall that these were $\psi_{\max} = E[-s(A \mid X)Y] + \gamma E[Y (\mathbbm{1}_{\{Y > M( A, X)\}} - \mathbbm{1}_{\{Y < M( A, X)\}})]$ and $\psi_{\min} = E[-s(A \mid X)Y] - \gamma E[Y (\mathbbm{1}_{\{Y > M( A, X)\}} - \mathbbm{1}_{\{Y < M( A, X)\}})]$. The efficient influence function for the functional $E[-s(A \mid X)Y]$ was derived in Newey1993EfficiencyModels. Thus, by a linearity property of efficient influence functions Kennedy2022, it only remains to find the efficient influence function of $\gamma E[Y (\mathbbm{1}_{\{Y > M( A, X)\}} - \mathbbm{1}_{\{Y < M( A, X)\}})]$.
The efficient influence function of $E[-s(A \mid X)Y]$ derived in Newey1993EfficiencyModels takes the following form:
where $\mu(A, X) \equiv E[Y \mid A, X]$, and $\mu'(A,X) \equiv \partial_a \mu(A,X)$. We then get the immediate corollary:
Equipped with the efficient influence functions, it is straightforward to propose estimators with desirable properties. The estimator will require estimating the unknown nuisance functions $\mu, \mu', s, M$. To ease the notational burden for the theoretical analysis, we simply analyze a sample split estimator where the nuisance functions are estimated on one split of the data, and on the second split, those estimates are plugged in to the efficient influence function at each data point, i.e.
To make use of all of the data, one can employ the now commonly utilized cross-fitting technique Chernozhukov2018b, where the data is randomly split into $K$ roughly equally sized folds $D_1, \ldots, D_K$. For each $k = 1,\ldots, K$, we compute nuisance estimates $\widehat{\eta} = (\widehat{\mu}, \widehat{\mu'}, \widehat{s}, \widehat{M})$ for $\eta = (\mu, \mu', s, M)$ on all folds except $D_k$, and plug in these estimates on fold $D_k$ (as in (ref)). The result of Theorem (ref), which outlines conditions under which the sample split estimator achieves asymptotic normality, will also apply to an analogous cross-fitted estimator.
By virtue of using an estimator based on the efficient influence function, the bias of the estimator only involves second-order nuisance estimation errors. Thus, $\sqrt{n}$ consistency is possible even if the nuisance functions can be estimated at the (slower than parametric) $n^{-1/4}$ rate. This makes it possible to conduct valid inference even when using nonparametric or data adaptive estimates of the nuisance functions, provided they are not converging too slowly to the truth. Based on the asymptotic normality of the estimators, it is straightforward to construct Wald-style confidence intervals. Since we are bounding upper and lower bounds, it is reasonable to construct one-sided confidence intervals. Explicitly,
are asymptotically valid $1 - \alpha$ confidence upper and lower bounds for $\psi_{\text{max}}$ and $\psi_{\text{min}}$, provided the variance estimates $\widehat{\sigma}_{\text{max}}$ and $\widehat{\sigma}_{\text{min}}$ converges to the true variances. Here, plug-in variance estimates can be used:
These plug-in estimates are consistent under the assumptions of Theorem (ref), and so the confidence bounds from Equation (ref) will be asymptotically valid.
As derived in the previous section, the bounds in the binary outcome case involve the term $\gamma E[\min\{1 - P(Y = 1 \mid A,X) , P(Y = 1 \mid A,X) \}]$. The presence of the minimum makes this quantity potentially non-smooth, so we instead rely on a smooth approximation. Specifically, when minima or maxima are involved, the LogSumExp (LSE) function is a popular choice (see Levis2023 for a recent example in causal inference). For a minimum of $k$ quantities, and any fixed $t > 0$,
Specialized to our setting, where we take a minimum of $p$ and $1-p$, we define
Thus, we will instead estimate (for a fixed $t$)
From Equation (ref), it is immediate that $\psi_{\max}^B \leq \psi_{\max,h_t}^{B}$ and $\psi_{\min}^B \geq \psi_{\min,h_t}^{B}$ for any $t > 0$, and the inequality gap shrinks as $t$ increases. At the same time, $h_t$ becomes less smooth as $t$ increases, and consequently, $E[h_t(P(Y = 1 \mid A,X))]$ becomes harder to estimate. Thus, in choosing $t$, there is a trade-off between approximation error and statistical estimation. A rigorously justified “optimal” choice for $t$ is outside the scope of this paper, but we refer the reader to Levis2023 for some additional discussion.
In the remainder of this section, we will derive the efficient influence function for the smoothed lower and upper bounds for the ADE under the sensitivity model with binary outcomes. Recall that these were $\psi_{\text{max},h_t}^{B}$ and $\psi_{\text{min},h_t}^{B}$. As in the continuous outcome case, the efficient influence function for the functional $E[-s(A \mid X)Y]$ was derived in Newey1993EfficiencyModels. Thus, it only remains to find the efficient influence function of $E[h_t\left(P(Y = 1 \mid A,X)\right)]$.
As before, we get the immediate corollary:
Again, we propose estimators based on the efficient influence function, and present sample-split versions that estimate the nuisances $\mu, \mu', s$. They are as follows: (since $Y$ is binary, $\widehat{\mu}(A_i,X_i)$ and $\widehat{P}(Y = 1 \mid A_i,X_i)$ are equivalent):
Theorem (ref) establishes asymptotic normality of the estimators.
Similar to the continuous outcome case, we can construct asymptotically valid Wald-style confidence intervals for $\psi_{\max}^B$ and $\psi_{\min}^B$ (rather than the smoothed $\psi_{\text{max}, h_t}^B$ and $\psi_{\text{max}, h_t}^B$), if we account for the approximation error from the LSE function. The respective upper and lower bounds for the $1 - \alpha$ confidence intervals are
Again, plug-in variance estimates can be used:
These estimates will be consistent under the assumptions of Theorem (ref), and so the confidence bounds from Equation (ref) will be asymptotically valid.
The previous subsections introduced estimators and pointwise confidence intervals for a fixed value of $\gamma$. In this subsection, we briefly outline how to conduct simultaneous inference when we wish to conduct the sensitivity analysis over a bounded interval of values, i.e. $\gamma \in [\gamma_l, \gamma_u]$. Conveniently, the nature of the resulting estimands take the form $a \pm \gamma b$, and we can easily construct Wald confidence intervals for $a$ and $b$ under the same assumptions as in the previous subsections. Thus, by the union bound, a uniform $(1-\alpha)$% confidence band can be straightforwardly constructed for $\gamma \in [\gamma_l, \gamma_u]$ by constructing $(1-\alpha/2)$% confidence intervals for $a$ and $b$ and concatenating accordingly. In our setting, $a$ corresponds to $E[-s(A \mid X)Y]$ and $b$ corresponds to either $E[Y (\mathbbm{1}_{\{Y > M( A, X)\}} - \mathbbm{1}_{\{Y < M( A, X)\}})]$ (continuous $Y$) or $E[h_t(P(Y = 1 \mid A,X))]$ (binary $Y$). Wald-style $(1-\alpha/2)$% confidence intervals for $a$ and $b$ can be constructed using the respective efficient influence functions. Alternatively, a multiplier bootstrap approach could be implemented (see Kennedy2019a and Zhang2022a for recent applications in causal inference). However, we do not pursue that direction as the proposed approach is much simpler.
We now evaluate the finite sample performance of the proposed methods through a simulation study. In the simulation, we assess the coverage of confidence intervals of the true ADE when there is unmeasured confounding based on sensitivity analyses at different choices of $\gamma$.
We conduct simulations corresponding to two different dose distributions, two outcome types (binary and continuous), and three different strengths of unmeasured confounders, yielding $ 2 \times 2 \times 3 = 12$ different settings. For all 12 settings, we draw the confounders $X \sim \text{Unif}[0,1]^d$, $d = 5$, and draw $U \mid X \sim \text{Bern}(\Phi(\sin(X_1+X_2)))$. We draw the dose from a conditional density that is Gaussian or Gamma. For the Gaussian case, $A \mid X, U \sim N(\theta^TX+\zeta U, 1)$. For the Gamma case (shape and rate parametrization), $A \mid X, U \sim \text{Gamma}(13, 8 + \theta^T X - \zeta U)$. $\zeta$ is set to $\log(2)$ for all simulations, and here $\theta$ are randomly drawn coefficients from $N(0,1)$, redrawn at each iteration. For the outcome model, we consider different settings for a binary outcome and for a continuous outcome. For the continuous outcome, we draw $Y \mid A, X, U \sim N(\eta A + \beta^T X + \delta U + \eta_{AX}AX, 1)$. For the binary outcome, we use a probit model and draw $Y \mid A, X, U \sim \text{Bern}(\Phi(\eta A + \beta^T X + \delta U + \eta_{AX}AX))$. $\delta$ is varied in $\{2, 3, 4\}$. For both outcome models, the $\beta$ coefficients are randomly drawn at each iteration from $N(-1,1)$. The interaction coefficients $\beta_{AX}$ are randomly drawn from $N(0,1/4)$, redrawn at each iteration. In the simulation, a higher $U$ leads to a higher chance of both a higher dose and outcome. The derivatives of the conditional expectations $E[Y \mid A, X, U]$ are available in closed form, and so the “ground truth” average derivative effects are approximated by drawing $10^7$ Monte-Carlo samples from the joint distribution of $(A, X, U)$ and computing the sample average of the derivative of $E[Y \mid A, X, U]$.
The nuisance estimates required to compute the estimator include the score function $s(a \mid x)$, the conditional mean $\mu(a,x)$ and its derivative $\mu'(a,x)$, and the conditional median $M(a,x)$ for the continuous outcome case. For the binary outcome case, we set $t = 50$ for computing the LSE function approximation. We use the R packages drape and xgboost for nuisance estimation. Specifically, we estimate scores $s(a \mid x)$ and conditional mean derivatives $\mu'(a,x)$ using adaptations of the methods from the drape package Klyne2023AverageLearning. The methods proposed by Klyne2023AverageLearning can re-smooth any first-stage regression $\widehat{\mu}(a,x)$ estimator to produce a differentiable version to obtain an estimate $\widehat{\mu}'(a,x)$, and model the conditional distribution $f(a \mid x)$ through a location-scale model to estimate $s(a \mid x)$. Hyperparameters for these methods were chosen in the same manner as the simulations in Klyne2023AverageLearning. To fit conditional means $\mu(a,x)$ and medians $M(a,x)$, we use gradient boosted trees as implemented in the xgboost package with default hyperparameters and the appropriate loss function -- squared error for estimating the conditional mean of a continuous variable, logistic loss for estimating the conditional mean of a binary variable, and absolute error for estimating the conditional median of a continuous variable. To make use of the full data sample, we implement 5-fold cross-fitting.
Table (ref) collects coverage results for pointwise 95% confidence intervals of the sensitivity analysis procedures at varying levels of $\gamma$. One can verify that the $\gamma$ at which the sensitivity analysis model (ref) holds (and thus the procedure will be valid) is between $0.5\log(2)$ and $\log(2)$. Therefore, it is not surprising to see in Table (ref) that the 95% sensitivity analysis confidence intervals can severely undercover when $\gamma$ is taken to be $0$ (no unmeasured confounding) or $0.25 \log(2)$, as both of these are less than the lower bound $0.5\log(2)$. This also gives some reassurance that although we have not established sharpness of the analytic bounds, the bounds can still be informative. In addition, one may notice that as the strength of the unmeasured confounders impact on the outcome, measured through $\delta$, increases, the sensitivity analysis intervals cover less. This is expected, as the sensitivity model we consider only restricts $U$'s impact on the treatment. Thus, the sensitivity analysis must protect against arbitrary dependence between $U$ and the (potential) outcomes, i.e. arbitrarily large values of $\delta$. Consequently, it is reasonable to expect that if $\delta$ were to be increased further, the coverage rate of the sensitivity analysis bounds for $\gamma \geq 0.5 \log(2)$ would move closer towards but not necessarily reach the nominal level.
We illustrate the methodology for binary outcomes using an empirical example studying the extent to which household income affects a child's educational attainment Lundberg2023. The data we use comes from the National Longitudinal Survey of Youth 1997 cohort (NLSY97). The NLSY97 is a dataset consisting of a probability sample of U.S. youths ages 12–17, starting in 1997, who were followed up through 2019. We largely follow Lundberg2023 in pre-processing the data. The treatment variable of interest is reported total gross household income in 1996, when the respondents were age 12-17. Lundberg2023 logged and adjusted these measures to 2022 dollars using the Consumer Price Index. Those without income measurements are dropped, as are households coded as the maximum and minimum income values, as these represent upper and lower cutoffs, not actual incomes. The outcome of interest is a report of enrollment in any college up to age 21, which is binary. Those that did not complete a survey at ages 19–21 are omitted. Following Lundberg2023, four measured confounding variables are adjusted for: race, gender, parents’ education, and wealth. The racial categories from 1997 were Hispanic, Non-Hispanic Black, and Non-Hispanic white or other. Parents’ education is categorical, with the 3 categories no parent completed college, one parent completed college, or two parents completed college. Wealth is the log of household net worth reported by the parent in 1997, also adjusted to 2022 dollars. There are 5219 individuals in the final, processed dataset. Unfortunately, there may be confounders beyond race, gender, parents' education, and wealth that are not measured that affect both household income and propensity to pursue higher education. These might include things like innate ability or geographic location, both of which could be strongly related to household income and propensity to attend college. Thus, we implement our sensitivity analysis to assess the impact of hypothetical unmeasured confounders on the statistical conclusions. We use the same nuisance estimators as in the simulations for a binary outcome. The results of the sensitivity analysis are reported in Figure (ref). Each black dot represents a point estimate of an upper or lower bound at some value of $\gamma$. The shaded regions depict 95% confidence intervals. Assuming no unmeasured confounding, the point estimate for the ADE is $0.108$, and the 95% confidence interval $[0.075, 0.141]$. Recall that the unit of the outcome is a percentage, and the treatment is the log of household income. In words, this means that on average, an increase of income of $\delta$ on the log scale might be expected to increase the propensity of attending any college by $\delta \times 10$ percentage points, assuming no unmeasured confounding. As $\gamma$ increases, lower point estimates and confidence bounds decrease. The point estimate ultimately crosses 0 at $\gamma = 0.323$, the 95% pointwise confidence interval crosses 0 at $\gamma = 0.222$, and the 95% uniform confidence interval crosses 0 at $\gamma = 0.197$.
We now illustrate the methodology for continuous outcomes using an empirical example studying the extent to which petrol prices affect the demand for petrol, which was also studied in Chernozhukov2022a. The data come from the Canadian National Private Vehicle Use Survey. We preprocessed the data in an identical fashion to Chernozhukov2022a, leaving $n = 5001$ households, each of which has an outcome -- log of the petrol consumption, covariates -- log age, log income, log distance, and other time, geographical, household indicators, and treatment -- log of petrol price per liter. In this example, the ADE measures the average price elasticity of petrol demand.
Instead of constructing pointwise confidence intervals, we concatenated the point estimates and standard errors for $E[-s(A \mid X)Y]$ as estimated in Chernozhukov2022a with estimates and standard errors for the correction term from Equation (ref) to produce simultaneous confidence bands. This exercise demonstrates the ease in which the sensitivity analysis can be conducted after (and completely separate from) a primary analysis assuming no unmeasured confounding has been completed. For estimation of the correction term, the only nuisance function is the conditional median. As in the simulation, we used the xgboost package with absolute error loss to fit the conditional median $M(a, x)$, in conjunction with 5-fold cross-fitting. The results are displayed in Figure (ref) when the generalized Dantzig selector (GDS) and Lasso estimators are used for estimating $E[-s(A \mid X)Y]$ as described in Chernozhukov2022a. Each black dot represents a point estimate of an upper or lower bound at some value of $\gamma$. The shaded region depicts the 95% simultaneous confidence intervals. Assuming no unmeasured confounding, the point estimates for the ADE are $-0.28$ and $-0.16$ for GDS and the Lasso, respectively. As $\gamma$ increases, upper point estimates and confidence bounds increase. The point estimates cross 0 at $\gamma = 0.924$ and $\gamma = 0.528$ and the 95% simultaneous confidence bounds cross 0 at $\gamma = 0.524$ and $\gamma = 0.2$ for GDS and Lasso, respectively.
In this paper, we have proposed a new sensitivity model for the ADE estimand, along with valid closed-form bounds, and estimators and confidence intervals for said bounds. The form of the bounds differ for continuous and binary outcomes, and are particularly convenient, allowing easy construction of uniform confidence intervals. The extent to which the bounds introduced in this paper are sharp, in the sense of Dorn2023, is unclear and is a promising direction for future research. To the best of our knowledge, sharp bounds for causal estimands under sensitivity models resembling (ref) do not exist beyond the binary treatment case. Another promising direction for further inquiry might be a calibration procedure, in the vein of Hsu2013 and McClean2024a. Calibrating the sensitivity analysis to observed confounders, for example, could potentially help researchers gauge whether a certain magnitude of $\gamma$ is plausible, though such a practice has limitations. Finally, it may be of interest to study ordinal rather than continuous exposures. Such exposures may arise, for example, when doses of drugs are prescribed at a finite number of ordered levels.
We thank Abhinandan Dalal, Zhihan Huang, Ziang Niu, Zhimei Ren, Dylan Small, Eric Tchetgen Tchetgen, and participants at ACIC 2025 for helpful discussions and comments.