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.
87,263 characters · 22 sections · 45 citation commands
Long Story Short: Omitted Variable Bias in Causal Machine Learning
\setcounter{page}{0}
\thispagestyle{empty}
Unmeasured confounding is a pervasive issue in studies that aim to draw causal inferences from observational data. Such studies typically rely on a conditional ignorability (also known as unconfoundedness) assumption, which states that the treatment assignment is independent of potential outcomes given a set of observed covariates Rosenbaum1983, pearl:causality,AngristBook, imbens2015causal. This assumption, however, requires that there are no unobserved confounders influencing both the treatment and the outcome. When such variables are omitted from the analysis, empirical estimates may differ from the true causal effect of interest, giving rise to what is now commonly known as “omitted variable bias.”
The omitted variable bias (OVB) problem is one the most significant threats to the identification of causal effects. In the context of linear models, this bias amounts to the difference between the coefficients of the treatment variable from two distinct outcome regressions: one that controls only for observed covariates (the “short” regression) and another that would additionally control for unobserved variables (the “long” regression). Formulas characterizing this difference play a foundational role in statistics, econometrics, and related fields (see, e.g., discussions in classical and modern textbooks, such as goldberger1991course, AngristBook and wooldridge:text). Such results allow empirical analysts to understand and bound the maximum size of the bias, by making plausibility judgments on the magnitude of parameters that comprise the OVB formula.
But while linear models are widely used in applied work, they are often overly restrictive. For example, in the binary treatment case, using linear models when treatment effects are heterogeneous may yield unintuitive or even misleading estimates of the causal effects of interest aronow2016does,sloczynski2022interpreting. To address these limitations, many empirical analysts have turned their attention to more flexible nonlinear or nonparametric models, often leveraging modern machine learning techniques for estimation and inference vanderlaan:book,BCFH2013,DML,athey:grf. These tools offer the flexibility to capture complex relationships between variables, avoiding stringent functional form assumptions in causal effect estimation. Yet, we currently lack general OVB results for nonlinear models (whether parametric or nonparametric), as we have for the linear case. Our work provides such results.
In this paper we develop a general theory of omitted variable bias for a wide range of common causal parameters that can be identified as linear functionals of the conditional expectation function (CEF) of the outcome. Such functionals encompass many (if not most) of the traditional targets of investigation in causal inference studies, such as averages of potential outcomes, average treatment effects, average causal derivatives, and policy effects from covariate shifts. We allow for arbitrary treatment (e.g continuous or binary) and outcome variables. Our theory applies to general nonparametric models, while naturally allowing for (semi-)parametric restrictions (such as partial linearity) when such assumptions are made. Our formulation recovers well-known and familiar OVB results for linear models as a special case, and it can be seen as its natural generalization to nonlinear models. Importantly, we show that the general nonparametric bounds on the bias still have a simple and interpretable form.
More specifically, we first formalize the OVB problem in the nonparametric setting. Paralleling the linear case, we define the OVB as the difference between the “short” and “long” functionals of the outcome regression, where the former omits and the latter includes the latent variables. To derive the OVB, our construction then leverages the Riesz-Frechet representation of the target functionals, which allows us to rewrite the parameters of interest as weighted averages of the outcome regression, with weights given by the Riesz representers (RRs). We show that the OVB arises as a by-product of confounders introducing systematic errors in both the outcome regression and in the RRs for the parameter of interest. Furthermore, the bound on the bias has a simple characterization, depending only on the additional variation that latent variables create in the outcome regression and in the RRs. As a result, plausibility judgments on the maximum explanatory power of latent variables suffice to place overall bounds on the bias, simplifying the task of sensitivity analysis even when using nonparametric or otherwise complex models.
Although these general results may initially seem abstract to those not familiar with Riesz representation theory, in many leading examples the RRs in fact correspond to quantities that are well-known to empirical researchers. For instance, when estimating the average treatment effect in a partially linear model, the RR is the (variance scaled) residualized treatment, after “partialling out” the control covariates. Or, when estimating an average treatment effect in a general nonparametric model with a binary treatment, the RR is now given by another familiar quantity---the inverse probability of treatment weights (IPTW). In such cases, we show that the bounds on the bias can be reparameterized in terms of simple percentage gains in variance explained (or precision) in the treatment and the outcome regression due to unmeasured confounders, again facilitating the interpretation and use of the OVB formulas in practice. We further help analysts make plausibility judgments on the magnitude of sensitivity parameters by means of comparison of the relative strength of unobserved confounders against the strength of observed covariates.
Finally, we provide statistical inference for these bounds using debiased machine learning (DML) and auto-DML DML, RR-local, Adverse, AutoDML-Lasso. Our construction makes it possible to use modern machine learning methods for estimating the identifiable components of the bounds, including regression functions, Riesz representers, the norm of regression residuals, and the norm of Riesz representers. These results enables flexible and efficient statistical inference on the bounds, allowing researchers to perform sensitivity analyses against unmeasured confounding in a flexible class of machine-learned causal models using simple and interpretable tools.\footnote{ Here we provide DML-based statistical inference on the bounds, but we note that our approach can also be used with classical parametric and nonparametric estimation methods.}
Our work is most closely related to the literature that derives OVB formulas for linear models, such as those found in traditional textbooks and recent extensions goldberger1991course,AngristBook,wooldridge:text,frank:smr2000,oster2019unobservable,CH2020. We advance this literature by providing analogous, easily explainable OVB formulas for a broad and rich class of causal parameters, all for general nonlinear models, with or without further parametric restrictions. Importantly, we provide a single unifying framework that covers all these cases, and that can be easily specialized depending on the target parameter and on whether additional parametric assumptions (if any) are made. We further advance the OVB literature by providing flexible and efficient statistical inference methods, leveraging modern machine learning algorithms with debiased machine learning.
More broadly, our work is related to the extensive literature on sensitivity analysis against unmeasured confounders. Here we highlight the key differences between our approach and existing methods, while relegating a more detailed review to the Appendix, Section (ref). First, many prior works on sensitivity analysis either focus exclusively on binary treatments (e.g., rosenbaum2002gamma,tan:jasa2006,masten:e2018,kallus2019interval,zhao:jrssb2019,bonvini2021sensitivity), target a single estimand of interest, such as a causal risk ratio ding:ecm2016,evalue, or impose parametric assumptions on the observed data or on the nature of unobserved confounding rosenbaum1983assessing,imbens2003sensitivity,dorie2016flexible,cinelli:icml2019. Our approach differs from these in that (i) it is not limited to binary treatments, (ii) it covers a broader range of target parameters, such as average causal derivatives and average policy effects from covariate shifts, and (iii) it does not require parametric assumptions on the observed data nor on the nature of confounding.
Even if we focus solely on the important special case of estimating an average treatment effect (ATE) with a binary treatment, our OVB results usefully complements other seminal approaches on this problem such as those of rosenbaum2002gamma or the marginal sensitivity models of tan:jasa2006. Whereas such approaches limit the strength of confounding through its impact on the worst case change that confounders could cause in the odds ratio of treatment assignment---a quantity economists rarely focus on---our approach limits the strength of confounding through its impact on the gains in precision in the treatment regression, a measure of explanatory power similar in nature to a simple $R^2$ in a linear model. Moreover, even in stylized models of treatment assignment (e.g, a logistic model with a Gaussian latent confounder), worst-case approaches such as the ones in rosenbaum2002gamma and tan:jasa2006 have a naturally unbounded sensitivity parameter, no matter how small the actual degree of confounding is, whereas our approach does not suffer from this problem (see Section (ref) of the Appendix for an example).
Our OVB-based approach also differs from traditional sensitivity analyses in that it derives the exact OVB formula for the target parameters we cover. For example, our results show that the bias of the ATE in the binary treatment case is not determined by deviations on the odds of treatment; rather, it is determined by three quantities: (i) the maximum explanatory power of confounders in the treatment regression, as given by gains in precision, (ii) the maximum explanatory power of confounders in the outcome regression, as given by gains in variance explained, and (iii) by the correlation of errors in the regression function and the IPTW. Therefore, beyond being a tool for sensitivity analysis, OVB results such as ours provide a precise characterization of the bias, and reveal that any alternative approach that parameterize deviations from unconfoundedness in a different way can only affect the bias insofar as it constraints these three quantities.
Section (ref) presents our method in the simpler context of partially linear models. The results in that section serve not only as an accessible introduction to the main ideas of our general framework, but are also important in their own right, since partially linear models are widely used in applied work. Section (ref) derives the main results of the paper---we characterize and bound the omitted variable bias for continuous linear functionals of the conditional expectation function of the outcome, based on their Riesz representations, all for general, nonparametric causal models. In Section (ref) we construct high-quality inference methods for the bounds on the target parameters by leveraging recent advances in debiased machine learning with Riesz representers. Section (ref) demonstrates the use of our tools to assess the robustness of causal claims in a detailed empirical example that estimates the average treatment effect of 401(k) eligibility on net financial assets. Section (ref) concludes with suggestions for possible extensions. The Appendix contains all proofs, provides a more extensive literature review, as well as an additional empirical example that illustrates sensitivity analyses for average causal derivatives with continuous treatments.
All random vectors are defined on the probability space with probability measure ${\mathrm{P}}$. We consider a random vector $Z=(Y,W)$ with distribution $P$ taking values $z$ in its support $\mathcal{Z}$; we use $P_{V}$ to denote the probability law of any subvector $V$ and $\mathcal{V}$ denote its support. We use $\| f \|_{P,q} = \| f(Z)\|_{P,q}$ to denote the $L^q(P)$ norm of a measurable function $f: \mathcal{Z} \to \mathbb{R}$ and also the $L^q(P)$ norm of random variable $f(Z)$. For a differentiable map $x \mapsto g(x)$, from $\mathbb{R}^d$ to $\mathbb{R}^k$, $\partial_{x'} g$ abbreviates the partial derivatives $(\partial/\partial x') g(x)$, and $\partial_{x'} g(x_0)$ means $\partial_{x'} g (x) \mid_{x = x_0}$. We use $x'$ to denote the transpose of a column vector $x$; we use $R^2_{U\sim V}$ to denote the $R^2$ from the orthogonal linear projection of a scalar random variable $U$ on a random vector $V$. We use the conventional notation $dL/dP$ to denote the Radon-Nykodym derivative of measure $L$ with respect to $P$.
To fix ideas, we begin our discussion in the context of partially linear models (PLM). These results not only provide the key intuitions and the building blocks for the general case of nonseparable, nonparametric models of Section (ref), but they are also important in their own right, as these models are widely used in applied work.
Consider the partially linear regression model of the form
Here $Y$ denotes a real-valued outcome, $D$ a real-valued treatment, $X$ an observed {vector} of covariates, and $A$ an unobserved vector of covariates. We refer to $W: = (D,X,A)$ as the “long” list of regressors, and to equation ((ref)) as the “long” regression. For exposition purposes, we assume the error term $\epsilon$ obeys ${\mathrm{E}}[\epsilon |D, X,A] =0$ and thus ${\mathrm{E}}[Y|D, X, A] = \theta D + f(X, A)$, though we note this assumption is not necessary.\footnote{We can also consider the case where $\theta D + f(X, A)$ is the projection of the CEF on the space of functions that are partially linear in $D$.}
Under the traditional assumption of conditional ignorability,\footnote{Along with consistency and usual regularity conditions.} we have that the regression coefficient $\theta$ identifies the average treatment effect of a unit increase of $D$ on the outcome $Y$, $$ {\mathrm{E}}[Y(d+1) - Y(d)] = {\mathrm{E}}[{\mathrm{E}}[Y|D = d+1, X, A] - {\mathrm{E}}[Y|D=d, X, A]] = \theta, $$ where $Y(d)$ denotes the potential outcome of $Y$ when the treatment $D$ is experimentally set to $d$. The problem, however, is that $A$ is not observed, and thus both the long regression, and the regression coefficient $\theta$ cannot be computed from the available data.
Since the latent variables $A$ are not measured, an alternative route to obtain an approximate estimate of $\theta$ is to consider the partially linear projection of $Y$ on the “short” list of observed regressors $W^s:= (D,X)$, as in,
where here we do not make the assumption that the regression is correctly specified, and thus the error term simply obeys the orthogonality condition ${\mathrm{E}}[ \epsilon_s (D- {\mathrm{E}}[D\mid X])]=0$. Following convention, we call equation ((ref)) the “short regression.” We can then use the “short” regression parameter $\theta_s$ as a proxy for $\theta$.
Evidently, in general $\theta_s$ is not equal to $\theta$, and this naturally leads to the question of how far our “proxy” $\theta_s$ can deviate from the true inferential target $\theta$. Our goal is, thus, to analyze the difference between the short and long parameters---the omitted variable bias (OVB): $$ \theta_s - \theta, $$ and perform inference on this bias under various hypotheses on the strength of the latent confounders $A$.
Recall that, using a Frisch-Waugh-Lovell partialling out argument, one can express the long and short regression parameters, $\theta$ and $\theta_s$, as the linear projection coefficients of $Y$ on the residuals $D- {\mathrm{E}}[D\mid X,A]$ and $D- {\mathrm{E}}[D\mid X]$, respectively. That is,
where here we define $$ \alpha(W) := \frac{D- {\mathrm{E}}[D\mid X,A]}{{\mathrm{E}} (D- {\mathrm{E}}[D\mid X,A ])^2 }, \quad \alpha_s(W^s) := \frac{D- {\mathrm{E}}[D\mid X]}{{\mathrm{E}} (D- {\mathrm{E}}[D\mid X])^2 }. $$ For reasons that will become clear in the next section, we can refer to $\alpha(W)$ and $\alpha_s(W^s)$ as the “long” and “short” Riesz representers (RR).\footnote{We are deliberately introducing Riesz representers in this section to smooth the transition to the general case. The formulation in terms of Riesz representers is a key innovation of this paper and it has not appeared in previous works on omitted variable bias. }
Now let $g(W) := {\mathrm{E}}[Y\mid D, X,A]$ and $g_{s}(W^s) := \theta_s D + f_s(X)$ denote the long and short regressions, respectively. Using the orthogonality conditions in ((ref)) and ((ref)), we can further express $\theta$ and $\theta_s$ as
Our first characterization of the OVB is thus as follows, where we use the shorthand notation: $g =g(W)$, $g_s = g_s(W^s)$, $\alpha= \alpha(W)$, and $\alpha_s = \alpha_s(W^s)$.
This result for partially linear models is new and it naturally generalizes the traditional OVB formula for linear models. It is worth noting that the proof of Theorem (ref) does not rely on the assumption that the long regression is partially linear, even though this assumption was made for expository purposes. In general, if we define both $g$ and $g_s$ to be projections of $Y$ onto the space of functions that are partially linear on $D$, the results of the theorem still hold.
Sensitivity analysis requires making plausibility judgments on the values of the sensitivity parameters. Therefore, it is important that such parameters be well-understood, and easily interpretable in applied settings. Here we show how the bias of Theorem (ref) can be reparameterized in terms of conventional $R^2$s.
Recall that, when the CEF is not linear, a natural measure of the strength of relationship between some variable $W$ and another variable $V$ is the nonparametric $R^2$---also known as Pearson's correlation ratio pearson1905general,doksum1995nonparametric: $$ \eta^2_{V\sim W} := R^2_{V\sim {\mathrm{E}}[V|W]} = \operatorname{Var}({\mathrm{E}}[V|W])/\operatorname{Var}(V) = \frac{\operatorname{Var}(V)-{\mathrm{E}}[\operatorname{Var}(V|W)]}{\operatorname{Var}(V)}. $$ Further, the nonparametric partial $R^2$ of a variable $V$ with another variable $A$ given $X$ measures the additional gain in the explanatory power that $A$ provides, beyond what is already is explained by $X$. This also equals the relative decrease in the average residual variance:
We are now ready to rewrite the bound of Theorem (ref).
The bound is the product of the term $S^2$, which is directly identifiable (and thus estimable) from the observed distribution of $(Y,D,X)$, and the term $C^2_Y C^2_D$, which is not identifiable from the data, and needs to be restricted through hypotheses that limit the strength of confounding. The factors $C^2_Y$ and $C^2_D$ measure the strength of confounding that the omitted variables generate in the outcome and treatment regressions. More precisely,
Note how this parameterization simplifies the complexity of plausibility judgments. Researchers now need only to reason about the maximum explanatory power that unobserved confounders have in explaining treatment and outcome variation, as given by familiar $R^2$ measures, in order to place bounds on the size of the bias. Finally, in practice, both $\theta_s$ and $S^2$ need to be estimated from finite samples. This can be readily done using debiased machine learning, as we discuss in Section (ref).
We now derive the main results of the paper, and construct sharp bounds on the size of the omitted variable bias for a broad class of causal parameters that can be identified as linear functionals of the conditional expectation function of the outcome. Although more abstract, the presentation of this section largely parallels the special case of partially linear models given in Section (ref).
As a motivating example, consider the following nonparametric structural equation model (SEM):
where $Y$ is an outcome variable, $D$ is a treatment variable, $X$ is a vector-valued observed confounder variable, $A$ is a vector-valued latent confounder variable, and $\epsilon _Y, \epsilon_D, \epsilon_A$ are vector-valued structural disturbances that are mutually independent. This model has an associated Directed Acyclic Graph (DAG) pearl:causality as shown in Figure (ref).
The SEM above induces the potential outcome $Y(d)$ under the intervention that sets $D$ experimentally to $d$, $$ Y(d) := g_Y(d,X,A, \epsilon_Y). $$ The structural model also encodes a consistency assumption between observed and potential outcomes, $Y=Y(D)$. Additionally, the independence of the structural disturbances implies the following conditional ignorability condition:
which states that the realized treatment $D$ is independent of the potential outcomes, conditionally on $X$ and $A$. More generally, we can work with any causal inference framework that implies the existence of potential outcomes, the consistency of observed and potential outcomes, and such that the conditional ignorability assumption ((ref)) holds AngristBook,pearl:causality,imbens2015causal.\footnote{There are many structural models that satisfy the conditional ignorability assumption ((ref)); see e.g. pearl:causality and Figure (ref) for concrete examples.}
Under this set-up and when $d$ is in the support of $D$ given $X$, $A$, we then have the following (well-known) identification result $$ {\mathrm{E}}[Y (d) \mid X, A]= {\mathrm{E}} [Y(d) \mid D=d, X, A] = {\mathrm{E}} [Y \mid D=d, X, A] =: g(d,X,A), $$ that is, the conditional average potential outcome coincides with the “long” regression function of $Y$ on $D$, $X$, and $A$. Therefore, we can identify various causal parameters---functionals of the average potential outcome---from the regression function. Important examples include: (i) the average treatment effect (ATE) $$ \theta = {\mathrm{E}}[Y(1) - Y(0)] = {\mathrm{E}} [g(1,X,A) - g(0, X, A)], $$ for the case of a binary treatment $D$; and, (ii) the average causal derivative (ACD) $$ \theta = {\mathrm{E}}\left [\partial_d {\mathrm{E}} [Y(D)\mid X, A] \right] = {\mathrm{E}}[\partial_d g(D,X,A)], $$ for the case of a continuous treatment $D$.
In fact, our framework is considerably more general, and it covers any target parameter of the following form.
This formulation covers the two previous examples with scores $m(W, g) = g(1, X, A) - g(0, X, A)$ for the ATE and $m(W, g) = \partial_d g(D, X, A)$ for the ACD. The continuity condition holds under the regularity conditions provided in the remark below. We discuss many other examples of this form later in Section (ref).
The key problem is that we do not observe $A$. Therefore we can only identify the “short” conditional expectation of $Y$ given $D$ and $X$, i.e. $$ g_s(D, X) := {\mathrm{E}} [Y \mid D, X]. $$ With the short regression in hand, we can compute proxies (or approximations) $\theta_s$ for $\theta$. In particular, for the ATE, the short parameter consists of $$\theta_s = {\mathrm{E}} [g_s(1,X) - g_s(0,X)],$$ and for the ACD, $$\theta_s = {\mathrm{E}} [\partial_d g_s(D,X)].$$
In this general framework, the proxy parameter can also be expressed as the same linear functional applied to the short regression, $g_s(W^s)$.
In the two working examples this assumption is satisfied, since $m(W, g_s) = m(W^s, g_s)= g_s(1, X) - g_s(0, X)$ for the ATE and $m(W, g_s) = m(W^s, g_s) = \partial_d g_s(D, X)$ for the ACD. Section (ref) verifies this assumption for other examples.
Our goal is to characterize and provide bounds on the omitted variable bias (OVB), ie., the difference between the “short” and “long” functionals, $$ \theta_s - \theta, $$ under assumptions that limit the strength of confounding, and perform statistical inference on its size.
The key to bounding the bias is the following lemma that characterizes the target parameters and their proxies as inner products of regression functions with terms called Riesz representers (RR).
In the case of the ATE with a binary treatment, the representers are just the classical inverse probability of treatment (Horvitz-Thompson) weights: $$ \alpha(W) = \frac{D}{P(D=1\mid X,A)}-\frac{1-D}{P(D=0\mid X,A)}, \quad \alpha_s(W) = \frac{D}{P(D=1\mid X)}-\frac{1-D}{P(D=0\mid X)}. $$ This follows from change of measure arguments. While it may not be immediately obvious that $\alpha_s = E[\alpha|D, X]$, one can easily show that by applying Bayes' rule.
In the case of the ACD with a continuous treatment, using integration by parts we can verify that the representers are logarithmic derivatives of the conditional densities: $$ \alpha(W) = -\partial_d \log f(D \mid X,A ), \quad \alpha_s(W^s) = -\partial_d \log f(D \mid X). $$ We give more involved examples in the next section.
Using this lemma, we obtain the following characterization of the OVB and bounds on its size.
This is the main conceptual result of the paper, and it is new. It covers a rich variety of causal estimands of interest, as long as they can be written as linear functionals of the long regression. We analyze further examples of this class of estimands in Section (ref).
In the same spirit of Section (ref), we can further derive useful characterizations of the bounds.
This generalizes the result of Corollary (ref) to fully nonlinear models, and general target parameters defined as linear functionals of the long regression. As before, the bound is the product of the term $S^2$, which is directly identifiable from the observed distribution of $(Y,D,X)$, and the term $C^2_Y C^2_D$, which is not identifiable, and needs to be restricted through hypotheses that limit strength of confounding.
Here, again, the terms $C^2_Y$ and $C^2_D$ generally measure the strength of confounding that the omitted variables generate in the outcome regression and in the treatment:
Likewise, we have the same useful interpretation of $C^2_Y$ as the nonparametric partial $R^2$ of $A$ with $Y$, given $D$ and $X$, namely, $C^2_Y = \eta^2_{Y \sim A\mid D, X}$. The interpretation of $1-R^2_{\alpha \sim \alpha_s}$ can be further specialized for different cases, as follows.
And an analogous interpretation applies for average causal derivatives.
Beyond making direct plausibility judgments on the strength of confounding using the above quantities, analysts can also leverage judgments of relative importance of variables to bound the size of the bias (see, e.g. imbens2003sensitivity,CH2020). For instance, if one has reasons to believe that $A$ would not generate as much gains in explanatory power as certain key observed covariates $X_j$, this can be used to formally place bounds on the strength of confounding due to $A$. This allows one to assess, for instance, whether confounders as strong or stronger then observed covariates would be sufficient to overturn an empirical result. We elaborate the benchmarking procedure formally in Section (ref) of the appendix and illustrate its use in the empirical example. These results extend previous benchmarking ideas for linear regression models to the general case.
We now provide theoretical details for a wide variety of interesting and important causal estimands. Recall that we use $W = (D,X, A)$ to denote the “long” set of regressors and $W^s = (D,X)$ to denote the “short” list of regressors.
Let us start with examples for the binary treatment case, with the understanding that finitely discrete treatments can be analyzed similarly.
The long parameter is a weighted average potential outcome (PO) when we set the treatment to $\bar d$, under the standard conditional ignorability assumption ((ref)). The short parameter is a statistical approximation based on the short regression. In this example, setting
Above we can consider $\mathcal{N}$ as small regions shrinking in volume with the sample size, to make the averages local, as in RR-local, but for simplicity we take them as fixed in this paper.
The long parameter is a weighted average treatment effect under the standard conditional ignorability assumption. In this example, setting
where the policy $\pi$ assigns a fraction $0\leq \pi(x) \leq 1$ of the subpopulation with observed covariate value $x$ to receive the treatment.
In what follows $D$ does not need to be binary. We next consider a weighted average effect of changing observed covariates $W^s$ according to a transport map $w^s \mapsto T(w^s)$, where $T$ is deterministic measurable map from $\mathcal{W}^s$ to $\mathcal{W}^s$. For example, the policy $$ (D,X,A) \mapsto (D + 1, X, A)$$ adds a unit to the treatment $D$, that is $T(W^s) = (D+1,X)$. This has a causal interpretation if the policy induces the equivariant change in the regression function, namely the counterfactual outcome $\tilde Y$ under the policy obeys ${\mathrm{E}} [\tilde Y |X,A] = g(T(W^s),A)$, and the counterfactual covariates are given by $\tilde W = (T(W^s),A)$.
We now turn to examples with continuous treatments $D$ taking values in $\mathbb{R}^k$. Consider the average causal effect of the policy that shifts the distribution of covariates via the map $W=(D,X,A) \mapsto (T(W^s),A) = (D + r t(W^s), X,A )$ weighted by $\ell(W^s)$, keeping the long regression function invariant. The following long parameter $\theta$ is an approximation to $1/r$ times this average causal effect for small values of $r$. This example is a differential version of the previous example.
Another example is that of a policy that shifts the entire distribution of observed covariates, independently of $A$. The following long parameter corresponds to the average causal contrast of two policies that set the distribution of observed covariates $W^s$ to $F_0$ and $F_1$, independently of $A$. Note that this example is different from the transport example, since here the dependence between $A$ and $W^s$ is eliminated under the interventions.
The following result establishes the validity of the OVB formulas and bounds for all examples.
As we have seen in Remarks (ref) and (ref), it may be useful to further specialize the interpretation of the sensitivity parameters $1-R^2_{\alpha\sim\alpha_s}$ for the many cases encompassed by the examples of Theorem (ref). As this would be an extensive task, we leave such specializations to future work.
The bounds for the target parameter $\theta$ take the form $$ \theta_{\pm} =\theta_s \pm |\rho| C_Y C_D S, \quad S^2 = {\mathrm{E}} (Y-g_s)^2{\mathrm{E}} \alpha^2_s. $$ The components $C_Y$, $C_D$ are set through hypotheses on the maximum explanatory power of omitted variables. Without further assumptions on the data generating process, $|\rho|$ is set to its upper bound of $|\rho|=1$, which is the most conservative scenario. Researchers may also investigate less conservative scenarios for $|\rho|$ based on, for example, empirical benchmarking as we illustrate in the empirical example. The estimable components of the bounds are $S$ and $\theta_s$. We can estimate these components via debiased machine learning (DML), which is a form of the classical “one-step” semi-parametric correction levit1975efficiency,hasminskii1978,pfanzagl:wefelmeyer, bickel1993efficient,newey1994,DML,DML-LR based on Neyman orthogonal scores we give for the these components, combined with cross-fitting, an efficient form of data-splitting.
For debiased machine learning of $\theta_s$, we exploit the representation $$ \theta_s = {\mathrm{E}} [m(W^s, g_s) + (Y- g_s) \alpha_s], $$ as in AutoDML-Lasso, AutoDML. This representation is Neyman orthogonal with respect to perturbations of $(g_s, \alpha_s)$, which is a key property required for DML. Another component to be estimated is $$ {\mathrm{E}} (Y-g_s)^2 =: \sigma^2_{s}, $$ which is also Neyman-orthogonal with respect to $g_s$. The final component to be estimated is ${\mathrm{E}} \alpha^2_s$. For this we explore the following formulation: $$ {\mathrm{E}} \alpha^2_s = 2 {\mathrm{E}} m(W^s, \alpha_s) - {\mathrm{E}} \alpha^2_s =: \nu^2_s, $$ where the latter parameterization is Neyman-orthogonal. Specifically Neyman orthogonality refers to the property:
where $\partial$ is the Gateaux (pathwise derivative) operator over directions $h \in L^2(P_{W^s})$.
Application of DML theory in DML and the delta-method gives the statistical properties of the estimated bounds under the condition that machine learning of $g_s$ and $\alpha_s$ is of sufficiently high quality, with learning rate faster than $n^{-1/4}$. The estimation relies on the following generic algorithm.
Therefore the estimators are defined as $$ \widehat \theta_s := \mathrm{DML} (\psi_{\theta}); \quad \widehat \sigma^2_s := \mathrm{DML} ( \psi_{\sigma^2}); \quad \widehat \nu^2_s := \mathrm{DML} (\psi_{\nu^2}); $$ for the scores
We say that an estimator $\hat \beta$ of $\beta$ is asymptotically linear and Gaussian with the centered influence function $\psi^o(Z)$ if $$ \sqrt{n} (\hat \beta - \beta) = \frac{1}{\sqrt{n}} \sum_{i=1}^n \psi^o(Z_i) + o_{{\mathrm{P}}} (1) \leadsto N(0, {\mathrm{E}} \psi^{o2}(Z)). $$
The application of the results in DML for linear score functions yields the following result.
The resulting plug-in estimator for the bounds is then: $$ \widehat \theta_{\pm} = \widehat \theta_s \pm |\rho| C_Y C_D \widehat S, \quad \widehat S^2 = \widehat \sigma^2_s \widehat \nu^2_s. $$ Confidence bounds for the bounds can be constructed using the following result.
We focus on the one-sided covering property stated in the theorem, since in applications the relevant hypotheses are typically one-sided. We can use further adjustments of stoye:CI to construct uniformly valid two-sided intervals.
The following remark discusses learning the regression function $g_s$ and the Riesz representer $\alpha_s$.
In this section we demonstrate the utility of our approach in an empirical example that estimates the average treatment effect of 401(k) eligibility on net financial assets pvw:94, pvw:95,DML. Our goal is to determine whether prior conclusions, reached under the assumption of conditional ignorability, are robust to plausible scenarios of unmeasured confounding. This example illustrates our bounding approach for the ATE in a partially linear model and in a nonparametric model with a binary treatment. In the Appendix we provide an additional example that estimates the price elasticity of gasoline demand blundell2012measuring, blundell2017nonparametric, chetverikov2017nonparametric and illustrates bounds for the average causal derivative with a continuous treatment.
A 401(k) plan is an employed sponsored tax-deferred savings option that allows individuals to deduct contributions from their taxable income, and accrue tax-free interest on investments within the plan. Introduced in the early 1980s as an incentive to increase individual savings for retirement, an important question in the savings literature is precisely to quantify the causal impact of 401(k) eligibility on net financial assets. Indeed, a naive comparison of net financial assets between those individuals with and without 401(k) eligibility suggests a positive and large impact: using data from the 1991 Survey of Income and Program Participation (SIPP), this difference amounts to \$19,559.
The problem of this naive comparison, however, is that 401(k) plans can be obtained only by those individuals that work for a firm that offers such savings option---and employment decisions are far from randomized. As an attempt to overcome this lack of random assignment, pvw:94, pvw:95, and more recently DML, leveraged the 1991 SIPP data to adjust for potential confounding factors between 401(k) eligibility and the financial assets of an individual. As explained in pvw:94, at least around the time 401(k) plans initially became available, people were unlikely to make employment decisions based on whether an employer offered a 401(k) plan; instead, their main focus were on salary and other aspects of the job. Thus, as a first approximation, whether one is eligible for a 401(k) plan could be taken as ignorable once we condition on income and other covariates related to job choice.
It is useful to think about causal diagrams pearl:causality that represent this identification strategy. One possible model is shown Figure (ref). Here the outcome variable, $Y$, consists of net financial assets;\footnote{Defined as the sum of IRA balances, 401(k) balances, checking accounts, U.S. saving bonds, other interest-earning accounts in banks and other financial institutions, other interest-earning assets (such as bonds held personally), stocks, and mutual funds less non-mortgage debt.} the treatment variable, $D$, is an indicator for being eligible to enroll in a 401(k) plan; finally, the vector of observed covariates, $X$, consists of: (i) age; (ii) income; (iii) family size; (iv) years of education; (iv) a binary variable indicating marital status; (v) a “two-earner” status indicator; (vi) an IRA participation indicator; and, (vii) a home ownership indicator. We consider that the decision to work for a firm that offers a 401(k) plan depends both on the observed covariates $X$, but also on latent firm characteristics, denoted by $A$; moreover, $X$, $A$, and $D$ are jointly affected by a set of latent factors $U$. Most importantly, note the assumption of absence of direct arrows, both from $A$ and $U$, to $Y$. Under such assumption, conditional ignorability holds adjusting for $X$ only. The story represented by the DAG of Figure (ref) is one way of rationalizing the identification strategy used in earlier papers.
The first three columns of Table (ref) shows the estimates for the average treatment effect (ATE) of 401(k) eligibility on net financial assets under this conditional ignorability assumption. For these estimates, we follow the same strategy used in DML, and we estimate the ATE using DML with Random Forests, considering both a partially linear model (PLM), and a nonparametric model (NPM).\footnote{We use Random Forest both for the outcome and treatment regression and estimate the parameters using DML with 5-fold cross-fitting. In order to reduce the variance that stems from sample splitting, we repeat the procedure 5 times. Estimates are then combined using the median as the final estimate, incorporating variation across experiments into the standard error as described in DML. } As we can see, after flexibly taking into account observed confounding factors, although the estimates of the effect of 401(k) eligibility on net financial assets are substantially attenuated, they are still large, positive and statistically significant (approximately \$9,000 for the PLM and \$8,000 for the NPM). With the nonparametric model, we further explore heterogeneous treatment effects, by analyzing the ATE within income quartile groups. The results are shown in Figure (ref). We see that the ATE varies substantially across groups, with effects ranging from approximately \$4,000 (first quartile) to almost \$18,000 (last quartile).
It is now useful to consider scenarios in which conditional ignorability fails. Figure (ref) presents one such scenario, where a violation of conditional ignorability is credible.\footnote{We note that Figure (ref) is just one example, and our sensitivity analysis results hold for {any} model in which conditional ignorability holds given observed variables and latent confounders.} Employers often offer a benefit in which they “match” a proportion of an employee's contribution to their 401(k) up to 5% of the employee's salaries. The model in Figure (ref) allows this “matched amount,” denoted by $M$, to be determined by unobserved firm characteristics $A$, observed worker characteristics $X$, and by 401(k) eligibility $D$. In this model, adjustment for $X$ alone is not sufficient for control of confounding. Instead, we now need to condition both on observed covariates $X$ and latent confounders $A$ for ignorability to hold.\footnote{Note that in this case the average treatment effect is still defined as ${\mathrm{E}}[ Y(1) - Y(0)]$. The relevant counterfactuals $Y(d)$ are obtained by setting $D=d$ for all descendants of $D$, that is $Y(d) = g_Y(d, M(d), X, \epsilon_Y)$, where $M(d) = g_M(d, F, X, \epsilon_M).$ } How strong would the omitted firm characteristics $A$ have to be in order to overturn our previous conclusions? And how plausible are the strengths revealed to be problematic? In what follows, we use our sensitivity analysis results to address these questions.
In reporting empirical results, the following definition will be useful.
Whereas standard errors, t-values or p-values communicate how robust the short estimate is to sampling errors, the idea of robustness values is to quickly communicate how robust the short estimate is to systematic errors due to residual confounding. For example, $\text{RV}_{\theta = 0,a=.05}$ measures the minimal strength on both confounding factors such that the estimated confidence bound for the ATE would include zero, at the 5% significance level.
Table (ref) illustrates our proposal for a minimal sensitivity reporting of causal effect estimates. Beyond the usual estimates under the assumption of conditional ignorability, it reports the robustness values of the short estimate. Starting with the PLM, the $\text{RV}_{\theta = 0, a = 0.05}= 5.4\%$ means that unobserved confounders that explain less than 5.4% of the residual variation, both of the treatment, and of the outcome, are not sufficiently strong to bring the lower limit of the confidence bound to zero, at the 5% significance level. Moving to the nonparametric model, we obtain a similar, but somewhat lower value of $\text{RV}_{\theta = 0, a = 0.05} = 4.5\%$. The RV thus provides a quick and meaningful reference point that summarizes the robustness of the short estimate against unobserved confounding---any postulated confounding scenario that does not meet this minimal criterion of strength cannot overturn the results of the original study.
We now proceed to construct a particular confounding scenario, based on the contextual details of the problem. We start with the assumption that $A$ explains as much variation in net financial assets as the total variation of the maximal matched amount of income (5%) over the period of three years (roughly the period over which the effect is measured).\footnote{This strategy is based on a suggestion by James Poterba.} In the worst case scenario, this would lead to an additional $3\%$ of total variation explained, resulting in a partial $R^2$ of outcome with omitted firm characteristics $A$ of $C_Y^2 = \eta^2_{Y \sim F | DX} = 4\%$.\footnote{ $ \eta^2_{Y \sim F | DX} = \frac{\eta^2_{Y \sim FDX} - \eta^2_{Y \sim DX}}{1 - \eta^2_{Y \sim DX}} = \frac{0.28 + 0.03 - 0.28}{1-0.28} \approx 4\%, $ } This amounts to a relative increase of approximately 10% in the baseline $R^2$ of the outcome regression of 28%. Following similar reasoning, and more conservatively, we posit that omitted firm characteristics can explain an additional $2.5\%$ of the variation in 401(k) eligibility, corresponding to a $22\%$ relative increase in the baseline $R^2$ of the treatment regression of 11.4%. For the partially linear model, this results in $1-R^2_{\alpha \sim \alpha_s}= \eta^2_{D\sim F\mid X}\approx 3\%$ (and also $C^2_D \approx 3\%$).\footnote{$ 1-R^2_{\alpha \sim \alpha_s} = \eta^2_{D \sim F |X} = \frac{\eta^2_{D \sim FX} - \eta^2_{D \sim X} } { 1- \eta^2_{D \sim X}} = \frac{0.114 + .025 - 0.114}{1-0.114} \approx 3\%. $} We adopt the same scenario for the nonparametric model, with the understanding that now this would correspond to gains in precision (see Remark (ref)). Since both $\eta^2_{Y \sim F | DX} \approx 4\%$ and $1-R^2_{\alpha \sim \alpha_s} \approx 3\%$ are below the robustness value of 5.4% (or 4.5%), we immediately conclude that such confounding scenario is not capable of bringing the lower limit of the confidence bound of the ATE to zero.
The exact bias, bounds, and confidence bounds on the ATE implied by the posited scenario are shown in Table (ref).\footnote{We use the same estimation procedure as described in footnote (ref).} Starting with the partially linear model, the confounding scenario has an estimated absolute value of the bias of $\$4,196$. Accounting for statistical uncertainty, we obtain a lower limit for the confidence bound of $\$2,497$. The results for the nonparametric model are qualitatively similar, with a bias of similar magnitude, and point estimates, bounds, and confidence bounds for the ATE shifted down by roughly one thousand dollars. Confidence bounds for group-wise ATEs can also be computed, and are shown in Figure (ref). Note how the bounds are still largely positive, with only a small excursion into the negative side in the case of the second quartile group. These results suggest that the main qualitative findings reported in earlier studies are relatively robust to plausible violations of unconfoundedness, such as the one specified by our confounding scenario.
A useful tool for visualizing the whole sensitivity range of the target parameter, under different assumptions regarding the strength of confounding, is a bivariate contour plot showing the collection of curves in the space of $R^2$ values along which the confidence bounds are constant imbens2003sensitivity, CH2020. These plots allow investigators to quickly and easily assess the robustness of their findings against any postulated confounding scenario. Here we focus on contour plots for the lower limit of the confidence bounds, as this is the direction of the bias that threatens the preferred hypothesis in this empirical example. Analogous contours can be constructed for the upper limit of the confidence bounds, and are omitted.
Starting with the partially linear model, the results are shown in Figure (ref). The horizontal axis describes the fraction of residual variation of the treatment explained by unobserved confounders, whereas the vertical axis describes the share of residual variation of the outcome explained by unobserved confounders. The contour lines show the lower limit of the confidence bounds $[l, u]$ for the ATE (see Theorem (ref)), given a pair of hypothesized values of partial $R^2$. Note $\text{RV}_{\theta=0, a=0.05}$ of Table (ref) is simply the point where the 45-degree line crosses the contour line of zero (red dashed line), offering a convenient summary of the critical contour. We can further place reference points on the contour plots, indicating plausible bounds on the strength of confounding, under alternative assumptions about the maximum explanatory power of omitted variables. The red triangle point on the plot---Max Match---shows the bounds on the partial $R^2$ as previously discussed, resulting in a lower limit of the confidence bound for the ATE of \$2,497, in accordance with Table (ref). Note here the correlation $|\rho|$ is set to its upper bound of 1.
Another approach to construct confounding scenarios is to use observed covariates to bound the plausible strength of unobserved covariates. For instance, in our empirical example, we know that employment decisions are largely driven by salary considerations. Similarly, salary is clearly an important determinant of net financial assets. One could therefore argue that it is implausible to imagine other latent firm characteristics that would be even a fraction as strong as the observed income of individuals, in terms of explanatory power in predicting 401(k) eligibility and net financial assets. Whenever such claims of relative importance can be made, they can be used to set plausible bounds on the strength of unmeasured confounding. Formal details of this benchmarking procedure are provided in Section (ref) of the Appendix.
The red diamonds of Figure (ref) shows the bounds on the strength of the latent variable $A$ if it were as strong as (i) income (1 x Income), (ii) whether a worker has an individual retirement account (1 x Part. in IRA), and (iii) whether the worker's family has a two-earner status (1 x Two Earners). Note that, apart from income, latent variables as strong as these covariates would result in a weaker confounding scenario than the one we have previously considered (Max Match). As for income, the worst-case bound indicates that omitted firm characteristics as important as income would indeed be sufficient to overturn the original results. However, one could argue such scenario to be implausible, as it is hard to imagine latent firm characteristics that would explain more variation in job choice than income itself. A more realistic, but still conservative, scenario is thus provided by the benchmark point 1/4 x Income, which shows the bound on the strength of $A$ if it were 25% as strong as income in predicting treatment and outcome variation. Note this scenario is comparable to the \emph{Max Match} scenario, and not enough to bring the lower limit of the confidence bound to zero.
All results of Figure (ref) were computed under the very conservative assumption that, given a pair of partial $R^2$ values for the latent variable $A$, the confounders enter both the outcome and treatment equations in a way that maximizes the bias, resulting in $|\rho|=1$. Although we can always construct such a confounder (absent further assumptions on the data generating process), it may be an unnatural scenario in practice, especially in nonlinear models.\footnote{ For an extreme example, consider the model $D = A^2$, $Y = \theta D + A$, with $A\sim N(0,1)$. Even though the latent variable $A$ nonparametrically explains 100% of the residual variation in both the treatment and the outcome equations, the nonlinearity of the confounding model attenuates this bias, making it effectively zero ($A^2$ is uncorrelated with $A$).} Thus, similar benchmarking procedures used for assessing the plausibility of the $R^2$ values can also be employed to calibrate judgments on the magnitude of $\rho$. Section (ref) of the appendix shows that in fact none of the observed covariates result in $|\rho|$ values exceeding 1/2. With this in mind, Figure (ref) presents the same contour plots as before, but now with $|\rho|$ set to a less conservative value of 1/2. Note how this substantially attenuates the bias, with the lower limits of the confidence bounds reaching approximately $\$4,600$ and $\$5,400$ for the Max Match and 1/4 x Income, respectively.
Sensitivity contour plots for the nonparametric model are similar but slightly more conservative, and are provided in Figure (ref). The interpretation of the contours is the same as before, with the main difference being that the horizontal axis now describes gains in precision instead of gains in variance explained (see, e.g, Remark (ref)).
In this paper we provide a general theory of omitted variable bias for continuous linear functionals of the conditional expectation function of the outcome---all for general, nonparametric, causal models, while naturally allowing for (semi-)parametric restrictions (such as partial linearity), when such assumptions are made. We allow for arbitrary (e.g., binary or continuous) treatment and outcome variables, and we show that the bounds on the bias depends only on the maximum explanatory power of latent variables. We provide theoretical details of many leading causal estimands, and, in particular, we derive novel bounds for the important special cases of average treatment effects in partially linear models, in nonparametric models with a binary treatment, as well as for average causal derivatives. Finally, we leverage the Riesz representation of our bounds to offer flexible statistical inference through (debiased) machine learning, with rigorous coverage guarantees. Therefore, we provide a concise and complete solution to the OVB problem and the bounding of its size, as well performing statistical inference on these bounds, for a rich and important class of causal parameters.
Our results can potentially be extended to nonlinear functionals, such as those arising in instrumental variable (IV) methods. For instance, consider a variant of the IV problem imbens:angrist:94, where the instrumental variable $Z$ is valid only when conditioning both on observed covariates $X$, and latent variables $A$. In this case, the IV estimand is given by the ratio of two average treatment effects, $$ \mathrm{IV}= \frac{\text{ATE}(Z \to Y)}{\text{ATE}(Z \to D)}. $$ Both the numerator and denominator can be bounded using the methods for the ATE proposed in this paper. Another interesting direction for future work is to consider causal estimands that are functionals of the long quantile regression, or causal estimands that are values of a policy in dynamic stochastic programming. When the degree of confounding is small, it seems possible to use the results in DML-LR to derive approximate bounds on the bias that can be estimated using debiased ML approaches.
\paragraph{Data availability.} All data is is publicly available in our GitHub repository: \url{https://github.com/carloscinelli/dml.sensemakr}. \paragraph{Conflict of interest.} There are no relevant financial or nonfinancial competing interests to report. \paragraph{Funding.} This research was partially funded by the Royalty Research Fund at the University of Washington.
{
}
\setcounter{page}{1}