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.
62,481 characters · 15 sections · 60 citation commands
Doubly-Valid/Doubly-Sharp Sensitivity Analysis for Causal Inference with Unmeasured Confounding
Investigating causal relationships using only observational data is often necessary when experimentation is infeasible. In the absence of a natural experiment, estimation often proceeds under the “unconfoundedness" assumption that all relevant confounders have been measured. Since this assumption is fundamentally untestable, it is imperative to conduct sensitivity analyses that explore how unobserved confounders might affect our inferences.
We consider sensitivity analysis under a nonparametric relaxation of unconfoundedness known as the marginal sensitivity model (MSM). This assumption allows for the existence of arbitrarily many unmeasured confounders, but posits that within each stratum of the observed covariates, measuring these confounders would not change the odds of receiving treatment by more than some factor. The posited factor can be varied by the data analyst to assess the robustness of conclusions to different magnitudes of unobserved confounding. The MSM was introduced by tan2006 and has been applied in kallus2018interval, zsb2019, kallus2020confoundingrobust, kallus_zhou2020, lee2021causal, rosenman2020combining, dornGuo2021sharp, rosenman2021designing,soriano2021interpretable, jin2021sensitivity, yin2021conformal, nie2021covariate, among others.
This paper draws a connection between sensitivity analysis under the MSM and classical results on distributionally robust optimization. In particular, we show that the problem of computing sharp upper and lower bounds on counterfactual means and treatment effects in the MSM reduces to the dual problem defining the conditional value at risk of a probability distribution rockafellar2000. Using this observation, we derive new “outcome regression" formulae for the optimal bounds on the average treatment effect (ATE) under the MSM. These formulae complement the “propensity weighting" formulae given in dornGuo2021sharp, jin2021sensitivity.
Based on this characterization, we propose new estimators for the sharp ATE bounds which we call the Doubly-Valid/Doubly-Sharp (DVDS) estimators. These estimators combine three nuisances --- a quantile regression, a propensity model, and a transformed-outcome regression --- to yield a number of desirable robustness properties:
No efficient or doubly-robust estimators were previously known under the MSM. Moreover, all previously-proposed confidence intervals under the MSM tan2006, zsb2019, dornGuo2021sharp, soriano2021interpretable have coverage guarantees only when first-stage nuisance functions are estimated by parametric models. In contrast, our method's estimation and inference can be valid even when nuisance functions are estimated by black-box machine learning techniques, possibly even inconsistently.
Double validity is --- to our knowledge --- a never-before-seen robustness property that uniquely enables credible sensitivity analysis. It is a refinement of the single validity guarantee in dornGuo2021sharp. Their sensitivity analysis estimates sharp bounds when both the propensity score and quantile regression are well-specified and estimates valid bounds when only the propensity score is well-specified. Our estimator offers additional protection by having an honest shot at valid inference when the quantile regression estimator is inconsistent and at valid estimates when the propensity score is inconsistent as well. These validity guarantees are particularly relevant as quantile regression may be more challenging than the usual propensity score and outcome regression problems.
Sensitivity analysis using DVDS estimators is especially convenient when the outcome variable is binary. The binary-outcome DVDS estimators have closed-form expressions in terms of the outcome regression and the propensity score, meaning they can be computed “for free" alongside the usual Augmented Inverse Propensity Weighting (AIPW) estimator of the ATE RRZDoubleRobust. Moreover, the validity of the DVDS confidence intervals holds under the same product-of-rates conditions typically used in the analysis of the AIPW estimator.
The rest of this paper is organized as follows: (ref) introduces the sensitivity assumption; (ref) characterizes the sharp partial identification region; (ref) defines the DVDS estimators; Sections (ref), (ref) present our statistical guarantees for point estimates and confidence intervals, respectively; (ref) illustrates the DVDS estimators empirically; and (ref) concludes. All proofs are deferred to the appendix.
In this section, we describe the mathematical setting of this paper and introduce the MSM.
We consider the Neyman-Rubin potential outcomes model with a binary treatment neyman, rubin1974. We assume $n$ units $(X_i, Y_i(0), Y_i(1), Z_i, U_i)$ are sampled independently from a distribution $P_{\textup{full}}$. Here, $(X_i, U_i) \in \mathbb{R}^d \times \mathbb{R}^k$ is a vector of confounders, $( Y_i(0), Y_i(1) )$ are real-valued potential outcomes, and $Z_i \in \{ 0, 1 \}$ is a binary treatment. However, the data analyst only observes $(X_i, Y_i, Z_i)$ where $Y_i = Y_i(Z_i)$. The confounders $U_i$ are also never observed, and their dimension $k$ may be unknown. We use $P$ to denote the distribution of the observables, so that $(X_i,Z_i,Y_i)\sim P$ are independent and identically distributed. We omit the subscript $i$ to denote a generic draw from $P_{\textup{full}}$ or $P$.
Our goal is to use the observed data to draw inferences about the average treatment effect $\psi_{\textup{ATE}}(P_{\textup{full}}) = \mathbb{E}_{P_{\textup{full}}}[Y(1) - Y(0)]$ and the counterfactual means $\psi_z(P_{\textup{full}}) = \mathbb{E}_{P_{\textup{full}}}[Y(z)], z \in \{ 0, 1 \}$. Extensions to the average treatment effect on the treated are presented in the Appendix.
Since we allow for unmeasured confounders, these quantities are not point-identified (meaning, not functions of $P$ alone) and cannot be consistently estimated. Nevertheless, they may be bounded if we impose some assumptions on the unmeasured confounders $U_i$. In this paper, we assume that measuring unobserved confounders could not change the odds of treatment by more than some factor $\Lambda \geq 1$. Formally, this restricts the distribution $P_{\textup{full}}$ to the marginal sensitivity model tan2006.
The identifying power of this assumption come from the bounded odds ratio condition (ref). When $\Lambda = 1$, this condition implies that measuring $U_i$ would not change the odds of treatment at all, i.e., the observed data is already unconfounded. As $\Lambda \rightarrow \infty$, stronger and stronger confounding is allowed. Condition (ref) is innocuous as it may always be satisfied by taking $U = ( Y(0), Y(1) )$. We frame the assumption in terms of an unobserved variable $U$ since it may be helpful to have some specific confounders in mind when choosing $\Lambda$. Meanwhile, any distribution not satisfying (ref) evidently cannot be the true complete-data distribution.
For any finite $\Lambda$, the assumption $P_{\textup{full}} \in \mathsf{MSM}(P, \Lambda)$ restricts the estimands $\psi_0, \psi_1, \psi_{\textup{ATE}}$ to a closed interval known as the partially identified set dornGuo2021sharp:
The quantities $\psi_{\diamond}^-, \psi_{\diamond}^+$ are called the sharp upper and lower bounds for $\psi_{\diamond}$, as they are the smallest upper bound and largest lower bound on $\psi_{\diamond}$ that can be derived from observations from $P$, assuming the MSM holds. The main goals of this paper are to develop robust estimators of these bounds and reliable confidence intervals for the partially identified set. For the rest of the paper, excluding (ref), we will fix some $\Lambda\geq1$ and study the corresponding sharp bounds and partially identified set.
Notations. We set $\tau=\frac{\Lambda}{\Lambda + 1}\in[0.5,1)$ throughout the paper. For a real number $t$, we set $\{ t \}_+ = \max \{ t, 0 \}$, $\{ t \}_- = \min \{ t, 0 \}$, and $\operatorname{sign}(t) = 1$ if $t \geq 0$ and $-1$ otherwise. For a real-valued function $f(x, y, z)$, we set $|| f ||_q = (\int |f|^q \, \textup{d} P)^{1/q}$ if $q \in [1, \infty)$ and $|| f ||_{\infty} = \inf \{ t \in \mathbb{R} \, : \, P( | f(X,Y,Z)| \leq t ) = 1 \}$. We also define the following functions:
We assume throughout that $Y$ has a finite mean and that the propensity score is bounded away from zero and one. These conditions ensure that the above functions are well-defined. Finally, we adopt the following convention regarding the symbols $\pm$ and $\mp$: any equation containing these symbols should be read twice, first with $\pm = +$ and then with $\pm = -$. Meanwhile, $\mp$ should always read with the opposite sign as $\pm$. This convention helps to unify formulae for upper and lower bounds.
In this section, we present two characterizations of the sharp partial identification bounds. The first is based on inverse propensity weighting (IPW) and extends the results of dornGuo2021sharp, jin2021sensitivity. The second is based on regression adjustment and is derived using classical results from distributionally robust optimization rockafellar2000. These characterizations motivate certain plug-in estimators of the identified set, which we combine in (ref) to form the DVDS estimators.
For simplicity, this section focuses on $\psi_1^+$, the sharp upper bound for $\mathbb{E}_{P_{\textup{full}}}[Y(1)]$. The same arguments apply to characterize $\psi_1^-, \psi_0^+$ or $\psi_0^-$, after exchanging $Y$ with $-Y$ and/or $Z$ with $1 - Z$. These bounds can then be subtracted to obtain sharp bounds on the ATE dornGuo2021sharp. Note that, in contrast, in other sensitivity models, subtracting sharp bounds on counterfactual means does not always give sharp bounds for the ATE yadlowsky2018bounds. This makes the MSM especially convenient for sensitivity analysis.
In the absence of unobserved confounding, counterfactual averages can be identified using the IPW formula $\mathbb{E}_{P_{\textup{full}}}[Y(1)] = \mathbb{E}[ YZ / e(X)]$. dornGuo2021sharp shows that a similar formula identifies $\psi_1^+$ once the propensity score $e(X)$ is replaced by an “adversarial" propensity score.
The interpretation of Proposition (ref) is the following. To maximize $\mathbb{E}_{P_{\textup{full}}}[Y(1)]$, the worst-case confounding structure in $\mathsf{MSM}(P, \Lambda)$ treats all “small" values of $Y(1)$ with the largest allowed probability and treats all “large" values of $Y(1)$ with the smallest allowed probability. The cutoff between large and small is determined separately for each covariate level $x$. Similar cutoff phenomena are observed in many other partial identification problems AronowLeeInterpretable, manski2003partial, lee2009training, jin2021sensitivity.
Despite its interpretability, it is not obvious how the identification formula from (ref) might be used for estimation. One issue is that (ref) does not specify the value of the worst-case propensity score $e_+(x, y)$ on the boundary event $y = Q_{\tau}(x, 1)$.
We now derive a new “propensity weighting" identification formula from (ref) that is more useful for estimation purposes. Since the worst-case propensity score $e_+$ satisfies $\mathbb{E}[ Z / e_+(X, Y) \mid X] = 1$, we may add $0 = \mathbb{E}[ Q_{\tau}(X, 1) Z / e(X) - Q_{\tau}(X, 1) Z / e_+(X, Y)]$ in (ref) to obtain:
In the final step, we plugged in the value of $e_+(X,Y)$ specified by (ref). This was possible even when $Y = Q_{\tau}(X, 1)$ because, on that event, the second term in (ref) is zero anyway.
The identification formula in (ref) is fully explicit and suggests a natural strategy for estimating $\psi_1^+$: first, estimate the propensity score $e$ and the conditional quantile $Q_{\tau}$ from data, then plug these estimated quantities into (ref).
The following Proposition shows that the natural strategy is singly valid: if the analyst's estimate of the propensity score is correct, then the plug-in estimator based on (ref) yields a valid (but possibly conservative) upper bound for $\psi_1^+$ even if the analyst estimates the quantile function $Q_{\tau}$ incorrectly. This property will play a crucial role in our estimation results to come.
The source of this robustness is a simple sign-matching argument: $( Y - \hat{Q}_{\tau}(X, 1) ) \Lambda^{\operatorname{sign} ( Y - \hat{Q}_{\tau}(X, 1) )}$ is always at least as large as $( Y - \hat{Q}_{\tau}(X, 1) ) \Lambda^{\operatorname{sign} ( Y - Q_{\tau}(X, 1) )}$, since $\Lambda^{\operatorname{sign} ( Y - \hat{Q}_{\tau}(X, 1) )}$ takes on its largest (smallest) value whenever $Y - \hat{Q}_{\tau}(X, 1)$ is positive (negative). Thus, using a misspecified cutoff is at least as conservative as using the true cutoff.
A second way of identifying counterfactual averages under unconfoundedness is via the regression adjustment formula $\mathbb{E}_{P_{\textup{full}}}[Y(1)] = \mathbb{E}[ZY + (1 - Z) \mu(X, 1)]$. In this section, we show that a similar formula will identify $\psi_1^+$ in the presence of unobserved confounding, after replacing the outcome regression $\mu(X, 1)$ by an “adversarial" outcome regression. The formula motivates another singly-valid plug-in estimator.
We start with the identity $\mathbb{E}_{P_{\textup{full}}}[Y(1)] = \mathbb{E}[ZY + (1 - Z) \mathbb{E}_{P_{\textup{full}}}[Y(1) \mid X, Z = 0]]$, which holds even in the presence of unmeasured confounding. The only unknown quantity here is the counterfactual regression $\mathbb{E}_{P_{\textup{full}}}[ Y(1) \mid X = x, Z = 0]$, which can be written as an integral relative to the distribution $F_{\textup{full}}(y \mid x, 0) = P_{\textup{full}}(Y(1) \leq y \mid X = x, Z = 0)$.
While the counterfactual distribution $F_{\textup{full}}(y \mid x, 0)$ is unobservable, the marginal sensitivity model allows us to bound it in terms of the distribution of the observed outcome, $F(y \mid x, 1)$. In particular, the bounded odds ratio condition (ref) and Bayes' theorem imply that their likelihood ratio must satisfy:
for almost every $x$. Moreover, these bounds are achievable when $P_{\textup{full}}$ makes the potential outcomes conditionally independent of one another. Therefore, the largest allowed counterfactual regression can be computed by maximizing the integral ((ref)) over measures $F_{\textup{full}}$ satisfying the likelihood ratio constraint ((ref)). This optimization problem was considered but not solved in tan2006.
After a simple reparameterization, the counterfactual regression maximization problem can be solved using results from distributionally robust optimization. The likelihood ratio constraint on $F_{\textup{full}}$ is equivalent to the condition that $F_{\textup{full}}(y \mid x, 0) = \Lambda^{-1} F(y \mid x, 1) + (1 - \Lambda^{-1}) G(y)$ for some distribution $G$ with $\textup{d} G(y)/\textup{d} F(y \mid x, 1) \leq 1/(1- \tau)$. Thus, the largest value of the counterfactual regression allowed under the MSM is:
In robust optimization, the problem of maximizing $\int y \, \textup{d} G(y)$ subject to the likelihood ratio constraint $\textup{d} G / \textup{d} F \leq 1/(1 - \tau)$ is known as the dual problem defining the level-$\tau$ conditional value at risk (CVaR) of the distribution $F$ ang_etal2017, rockafellar2000, Rockafellar2006Generalized, ruszczynski2006optimization. It can be solved explicitly in terms of the conditional quantiles of $F$:
Plugging the CVaR into (ref) yields our adversarial “regression adjustment" formula for $\psi_1^+$.
The interpretation of (ref) is the following. The adversarial counterfactual regression $\rho_+$ is a mixture between the average of the distribution $F(y\mid x,z)$ and the average of only the largest $100(1-\tau)\%$ of values under that distribution. As $\Lambda$ grows, the CVaR term both gets more weight and uses a more extreme tail average. In the limit, the CVaR converges to the upper bound of the support of $F(y \mid x, z)$. Because of its interpretation as an extremal average, the CVaR appears naturally in partial identification problems lee2009training, semenova2021generalized.
We now discuss how this identification formula might be used to form a plug-in estimator for $\psi_1^+$. Observe from (ref) that $\rho_+$ is the conditional mean of a certain transformed outcome:
Therefore, to estimate $\rho_+$, one can first learn $Q_{\tau}$ by quantile regression and then learn the conditional mean of the estimated transformed outcome $\Lambda^{-1} Y_i + (1 - \Lambda^{-1}) [ \hat{Q}_{\tau}(X_i, Z_i) + \tfrac{1}{1 - \tau} \{ Y_i - \hat{Q}_{\tau}(X_i, Z_i)]$. The resulting estimate of $\rho_+$ may then be plugged into the identification formula in (ref).
The following Proposition shows that this plug-in estimator for $\psi_1^+$ also has a single validity property: even if the quantile regression is misspecified, we still obtain a valid (but possibly conservative) upper bound for $\psi_1^+$ as long as the transformed outcome regression correctly learns the conditional mean of the estimated transformed outcome.
The explanation for this robustness is that the CVaR is the optimal value of the following minimization problem solved by the true quantile $Q_{\tau}$:
See koenker_bassett_1978 or rockafellar2000 for a proof. Therefore, using any function other than $Q_\tau$ yields a bound that is equal or larger.
As an aside, this derivation also directly motivates the consideration of more general sensitivity models. (ref) may be replaced by other ambiguity sets for $F_{\textup{full}}(y \mid x, 0)$, such as $f$-divergence or Wasserstein balls around $F(y \mid x, 1)$. The Rosenbaum model considered in yadlowsky2018bounds also falls in this class. The duality between robust optimization and convex risk measures ruszczynski2006optimization then implies the resulting adversarial regression would be the corresponding risk measure on the conditional outcome distribution. For example, had we used the Kullback-Leibler divergence, we would have obtained the entropic value at risk ahmadi2012entropic in place of CVaR. Moreover, like CVaR, the resulting risk measure will be expressible as the value of a minimization problem. Thus, identification formulae obtained in this way are always “singly valid."
In this section, we introduce the DVDS estimators and give some intuition for their properties, proving them formally in subsequent sections. These estimators combine the “propensity weighting" strategy from Section (ref) with the “regression adjustment" strategy from Section (ref) to simultaneously enjoy the robustness properties of both Propositions (ref) and (ref). This is analogous to how the AIPW estimator combines propensity weighting and regression adjustment to achieve double robustness under unconfoundedness.
The first step in computing the DVDS estimators is to learn the nuisance functions from Section (ref). We re-define these below in a generalized notation that accommodates either upper or lower bounds:
Let $\eta = (e, Q_+, Q_-, \rho_+, \rho_-)$ denote the full vector of DVDS nuisance functions.
While estimating the propensity score $e$ is a standard task in causal inference, estimating the other nuisances merits further discussion.
For continuous outcomes, we recommend estimating $Q_{\pm}, \rho_{\pm}$ in the manner described in (ref). First, $Q_{\pm}$ is estimated by any quantile regression method quantile_neural_networks,quantile_random_forest,generalized_random_forests,belloni2011quantile,koenker_bassett_1978. Then, $\rho_{\pm}$ is estimated using regression to learn the conditional mean of the estimated transformed outcome using the estimated quantile. In the Appendix, we show that this two-step approach can yield accurate estimates of $\rho_{\pm}$ despite first-step estimation error. If this second-step regression is performed by combining separate regression estimates of the $\mu$ and CVaR components, then, as $\Lambda$ tends to one, the DVDS estimator will tend to the AIPW estimator that uses the same $\mu$ and $e$-estimates.
For binary outcomes, estimating the additional nuisances is considerably simpler. The quantiles and CVaR of a Bernoulli($\mu(x, z)$)-distributed random variable are explicit functions of $\mu(x, z)$, so the nuisances $Q_{\pm}, \rho_{\pm}$ can be computed by directly plugging an estimate $\hat{\mu}$ of the outcome regression into the following formulae:
When these formulae are used, the DVDS estimators tend to the AIPW estimator as $\Lambda$ tends to one and they tend to Manski1990's “no assumptions" bounds as $\Lambda$ grows large.
For either continuous or binary outcomes, we assume that the nuisances are estimated using a standard sample-splitting technique known as cross-fitting. This technique is described in (ref) and has been applied in schick1986, robins_etal2008, zheng2011cross, doubleML (among others). Cross-fitting allows black-box machine learning methods to be used for nuisance estimation.
The second step in computing the DVDS estimators is to aggregate the estimated nuisances using the following recentered influence functions:
(ref) gives a complete description of this aggregation.
To provide some intuition for these recentered influence functions, we present two decompositions of $\phi_1^+$. For conceptual clarity, we assume in the discussion below that the nuisances $\hat{\eta}$ are fixed rather than estimated functions.
In summary, we expect the estimator $\hat{\psi}_1^+$ to have two chances at estimating sharp bounds and, if quantiles are misspecified, two chances at estimating valid bounds.
In this section, we present the formal asymptotic properties of the DVDS point estimators. For simplicity, we state our results only for the ATE bounds $(\hat{\psi}_{\textup{ATE}}^-, \hat{\psi}_{\textup{ATE}}^+)$, although the same results hold for the counterfactual mean bounds as well.
Our first main result on estimation shows that the robustness properties intuitively derived in (ref) hold even when the nuisances are estimated from data. This requires only mild regularity conditions, which are standard in semiparametric causal inference Kennedy2016, doubleML, TargetedLearningVDLRose.
We suppress cross-fitting notation from this assumption and the ones that follow. It should be understood that conditions stated in terms of “$\hat \eta$” are imposed on each of $\hat\eta^{(-k)}$, for $k=1,\dots,K$.
(ref) below states that correct specification of either $e$ or $\rho_{\pm}$ suffices for the DVDS estimators to estimate the sharp ATE bounds when quantiles are correctly specified. However, if quantiles are misspecified, then the DVDS estimators may still estimate valid bounds on the ATE as long as the propensity score is well-specified or the transformed outcome regression accurately learns the conditional mean of the estimated transformed outcome.
While the “double sharpness" result (ref) is a relatively standard multiple robustness property, the “double validity" result (ref) appears to be new. In our view, it enables uniquely credible sensitivity analysis by offering two chances at valid bounds even when consistency fails.
For binary outcomes, the condition $|| \hat{\rho}_+ - \varrho_{\pm}(\cdot,\cdot;\hat Q_\pm) ||_2 = o_P(1)$ will be satisfied whenever the outcome regression $\mu$ is consistently estimated. As a result, Theorem (ref) for binary outcomes may be phrased entirely in terms of the outcome regression and propensity score.
Our second main result on estimation shows that the DVDS estimators are asymptotically efficient when all first-stage nuisance parameters are estimated consistently with certain rates of convergence. These convergence rates accommodate nonparametric nuisance learners based on machine learning.
This result requires certain density assumptions. In the continuous case, we assume that $Y$ has a bounded conditional density, as required by many standard quantile regression methods (e.g., generalized_random_forests, or belloni2011quantile). In the binary case, we assume $\mu(X, Z)$ has a bounded density. This ensures that thresholding an accurate estimate of $\mu$ yields accurate estimates of $Q_{\pm}$.
The semiparametric efficiency bound $\mathbf{\Sigma}$ is the smallest asymptotic variance attainable by any estimator which converges locally uniformly to its asymptotic distribution. Moreover, $\textup{trace}(\mathbf{\Sigma})$ is an asymptotic lower bound on the local minimax mean squared error attainable by any estimator. If the rate conditions and bounds assumed in Theorem (ref) hold uniformly in a nonparametric model $\mathcal{P}$, the convergence in (ref) can be improved from locally uniform to globally uniform as in doubleML or kallus2020localized. We skip this detail in the analysis for brevity.
The rate conditions in (ref) will be satisfied if all first-stage nuisance parameters are estimated at rates faster than $n^{-1/4}$. This is a standard requirement in semiparametric inference. For the outcome regression, propensity score, and quantile regression nuisances, primitive conditions that imply this rate can be found in the rich literature on nonparametric function estimation (see generalized_random_forests, belloni2011quantile, gyorfi_etal2002, wainwright_2019, and references therein). Even the $L^{\infty}$ rate requirement on $\hat\mu$ needed to estimate (discrete) quantiles in the binary case is achievable in many cases stone1982optimal.
Convergence rates for $\hat{\rho}_{\pm}$ in the continuous case are slightly nonstandard since the transformed outcome is itself estimated from the data. Fortunately, this regression problem satisfies a conditional Neyman orthogonality doubleML, foster2020orthogonal:
In the Appendix, we show that this property allows us to essentially ignore estimation error in the transformed outcome when deriving convergence rates for $\hat{\rho}_{\pm}$, meaning the rate $|| \hat{\rho}_{\pm} - \rho_{\pm} ||_2 = o_P(n^{-1/4})$ may be achieved even when $\hat{Q}_{\pm}$ is itself estimated by black-box methods. olma2021truncated obtains similar results for two-step CVaR estimation using local polynomial regression methods.
In this section, we explain how the DVDS estimators can be used to set confidence limits on the endpoints of the partially-identified set or the entire set itself. While all previous approaches for inference in the MSM have required computationally-intensive bootstrap procedures dornGuo2021sharp, soriano2021interpretable, zsb2019, we propose using simple Wald-type intervals based on the following standard errors:
Here, we recall that $\phi_{\textup{ATE}}^{\pm}$ is the estimated recentered influence function defined in (ref).
The following Theorem shows that confidence limits based on these standard errors are asymptotically valid and optimal under the conditions of Theorem (ref). This is a straightforward consequence of asymptotic normality and semiparametric efficiency, so we view it as the inferential analog of double sharpness.
The theoretical guarantees for lower confidence bounds of the form $\hat{\psi}_{\textup{ATE}}^- - z_{1 - \alpha} \hat{\sigma}_-$ are analogous. Optimal one-sided confidence limits for the partially identified also provide optimal one-sided confidence limits for the true (unidentified) ATE. To obtain an asymptotic $100(1 -\alpha)\%$ two-sided confidence region for the entire identified set, one can intersect the $100(1 - \alpha/2)\%$ upper and lower confidence regions. By the union bound, this will cover the identified set with asymptotic probability at least $1 - \alpha$. It is possible to construct slightly refined intervals with asymptotic coverage exactly $1 - \alpha$ by accounting for the correlations between $\hat{\psi}_{\textup{ATE}}^+, \hat{\psi}_{\textup{ATE}}^-$ kallus2021assessing. However, we expect the scope of this refinement to be limited in practice as we typically see strong positive correlations between the bounds in our simulations.
One dissatisfying feature of (ref) is that it requires stronger rate conditions than those typically used for inference under unconfoundedness. For example, inference based on the usual AIPW estimator does not require any quantile regression or $L^{\infty}$ rate conditions.
Thus, our second inference result considers what happens when these extra rate conditions are not satisfied. In other words, we allow for the quantiles $Q_{\pm}$ to be estimated at a rate slower than $n^{-1/4}$, or even to be completely misspecified. In the binary case, this amounts to assuming only the standard product-of-rates condition used in the analysis of the AIPW estimator, with no $L^{\infty}$ rate condition. Because this result does not assume quantile consistency, we view it as the inferential analogue of double validity.
In the continuous case, (ref) assumes that $\hat{\rho}_{\pm}$ is not too far from $\varrho_{\pm}(\cdot, \cdot, \hat{Q}_{\tau})$, the conditional mean of the estimated transformed outcome. We show in the appendix that this is achievable even when the quantile regression is misspecified or converges at a much slower rate.
In our view, (ref) is surprising since Wald-type confidence intervals are typically invalid when asymptotic normality fails. We prove this result by showing that the Wald upper confidence bound based on $\hat{\psi}_{\textup{ATE}}^+$ is first-order equivalent to the confidence bound obtained from a certain bootstrap scheme. Then, we couple the bootstrap distribution of the DVDS estimator with that of an infeasible but asymptotically normal estimator $\bar{\psi}_{\textup{ATE}}^+$ to show that the quantiles of the former bootstrap distribution are always larger than the quantiles of the latter. Since $\bar{\psi}_{\textup{ATE}}^+$ is a well-behaved estimator, its bootstrap quantiles are asymptotically valid confidence limits. Therefore, the bootstrap quantiles based on the DVDS estimator are as well.
We make a technical remark on uniform validity. Under the stronger assumptions of (ref), the DVDS estimators are regular so the convergence in (ref) is automatically uniform over local (contiguous) neighborhoods of $P$. Under the weaker assumptions of (ref), the DVDS bounds may no longer be regular, so one may expect that the convergence in (ref) is only pointwise. However, this is not the case. Even when the DVDS estimators are nonregular, Wald confidence bounds retain their aforementioned uniform validity. We sketch the argument in our proof of (ref).
We next demonstrate DVDS in simulated examples and a case study.
In this section, we illustrate the DVDS estimators in two simulated examples.
In both cases, in a slight departure from our formal results, we use out-of-bag predictions instead of cross-fitting. This approach is also taken in athey2017efficient. All forests are fit using R package grf generalized_random_forests.
These computations were performed $1,000$ times in each example. For each replication, we also computed 95% Wald confidence intervals based on the DVDS estimators. For comparison, we also obtained bounds using the AIPW sensitivity analysis for the MSM proposed in zsb2019. The estimated bounds are plotted in (ref).
In the binary case, DVDS estimators are approximately unbiased for all considered values of $\Lambda$. Moreover, nominal level $95\%$ two-sided Wald confidence bounds had coverage ranging from $92.2\%$ (when $\Lambda = 1$) to $96.5\%$ (when $\Lambda = 2.1$) with average coverage of $95.3\%$. Meanwhile, the zsb2019 bounds are quite conservative in the binary outcome simulation.
In the continuous case, the DVDS upper bounds perform well but the lower bounds are somewhat conservative. Here, it appears that the random forest methods used for quantile and regression estimation did not adapt well to the outcome model used. Nevertheless, nominal level 95% two-sided Wald confidence bounds had reasonable coverage ranging from $91.7\%$ (when $\Lambda = 1$) to $98.1\%$ (when $\Lambda = 3$) with average coverage of $96.1\%$. The DVDS bounds improve only slightly over the zsb2019 bounds in this example. The difference may be more pronounced in examples with more heteroscedasticity.
In summary, these experiments largely validate the theoretical results presented in Sections (ref) and (ref). When nuisances are estimated well, DVDS point estimates and confidence intervals yield sharp inferences on the partially identified set even when machine learning methods are used. Meanwhile, when some nuisances are estimated poorly, the procedure errs on the side of conservatism.
We now apply the DVDS estimators in a real-data example. Specifically, we revisit ConnorsEtAl1996's influential finding that right heart catheterization (RHC) leads to lower survival rates among critically ill patients. Their data has been extensively analyzed in the causal inference literature (e.g. tan2020, ATTViaIPW, cui_tchetgen2019, lin1998), and was the first real-data application ever studied using the MSM tan2006.
We use a version of the data from BWSnooping. It consists of 5{,}735 observations on adult patients from five US hospital centers. The treatment $Z$ indicates whether a patient received RHC within 24 hours of admission, and the outcome $Y$ is 30-day survival rate. The data also contains measurements on a rich set of demographic and medical factors $X$ measured prior to treatment. For a complete list, see ATTViaIPW. Our estimand of interest is the ATE.
In our analysis, we estimate the propensity score $\hat{e}$ and outcome regression $\hat{\mu}$ using the Super Learner ensemble method SuperLearner. As base learners, we use logistic regression, random forest, gradient boosted trees, and support vector machine with Platt scaling.
The AIPW estimator based on our estimated nuisances yields an ATE estimate of $-5.4\%$ with a 95% confidence interval of $[-7.8\%, -3\%]$. Thus, assuming unconfoundedness, we find that RHC has a negative average effect on 30-day survival. cui_tchetgen2019, tan2020 also apply machine learning methods to this dataset and obtain similar estimates.
We compare the DVDS method to the AIPW-based sensitivity analysis from zsb2019. Since the outcome $Y$ in this example is binary, no nuisance functions beyond the propensity score and outcome regression are needed for either sensitivity analysis. We also compute (pointwise) 95% confidence intervals based on the DVDS estimator.
Estimated bounds for various values of $\Lambda$ are shown in (ref). Both sensitivity analyses find that fairly small amounts of unobserved confounding would suffice to explain away the negative point estimate ($\Lambda = 1.17$ for the zsb2019 method, $\Lambda = 1.23$ for the DVDS method). If we further account for statistical variability in the estimated bounds, the DVDS estimators would find statistical evidence for a negative treatment effect only under the assumption that unobserved confounders cannot change the odds of treatment by up to a factor of $\Lambda = 1.12$. Our bound estimates are comparable to those of tan2006 and suggest more caution is warranted than is conveyed by ConnorsEtAl1996's original sensitivity analysis.
This paper considered estimation of bounds on the ATE (and other causal quantities) under the assumption that measuring unobserved confounders does not change the odds of treatment by more than a factor of $\Lambda$. We characterized the sharp partial identification bounds in terms of CVaR regression and used this characterization to propose robust estimators of the partially identified set.
Perhaps the most novel phenomenon observed in our work is double validity and the associated inference-without-normality result of (ref). As far as we can tell, neither of these properties has been observed before, possibly because they do not arise in point-identified problems.
Since the literature on semiparametric estimation of partial identification bounds is relatively nascent tchetgen_tchetgen_shpitser2012, yadlowsky2018bounds, BonviniKennedy, semenova2021generalized, doubleMLSensitivity, kallus2021assessing, we expect that single and possibly double validity will be recurring phenomenon as the field grows. Since the first draft of this paper, double validity has even been observed in partial identification outside of sensitivity analysis kallus2022treatment,kallus2022s. For sensitivity analysis, as discussed in (ref), single validity should arise in any sensitivity model that bounds a convex distributional divergence between $P_{\textup{full}}(Y(1) \mid X = x, Z = 0)$ and $P(Y(1) \mid X = x, Z = 1)$. The MSM falls into this class, but numerous other choices have been popularized in the growing literature on distributionally robust optimization ben2013robust,bertsimas2018robust,esfahani2018data. This yields a rich family of sensitivity models that may be appropriate for different applications. Since the first draft of this paper, jin_etal2022 have explored such alternatives using $f$-divergences and indeed obtain single validity guarantees.
This material is based upon work supported by the National Science Foundation under Grant No. 1939704 and by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE-2039656. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. The authors thank Angela Zhou for helpful discussions while conceiving the idea for this paper and are grateful for comments from Bo Honor{\'e}, Michal Koles\'{a}r, and Mikkel Plagborg-Møller.