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.
47,683 characters · 11 sections · 18 citation commands
Orthogonal Machine Learning: Power and Limitations
\twocolumn[ \icmltitle{Orthogonal Machine Learning: Power and Limitations}
\icmlaffiliation{msr}{Microsoft Research New England, USA} \icmlaffiliation{mit}{Operations Research Center, MIT, USA}
\icmlcorrespondingauthor{Lester Mackey}{[email removed]} \icmlcorrespondingauthor{Vasilis Syrgkanis}{[email removed]} \icmlcorrespondingauthor{Ilias Zadik}{[email removed]}
\icmlkeywords{orthogonal machine learning, double machine learning, asymptotic normality, nuisance, semiparametric inference, partially linear regression, high-dimensional regression, Stein's lemma, treatment effect, causal inference}
\vskip 0.3in ]
\printAffiliationsAndNotice
The increased availability of large and complex observational datasets is driving an increasing demand to conduct accurate causal inference of treatment effects in the presence of high-dimensional confounding factors. We take as our running example demand estimation from pricing and purchase data in the digital economy where many features of the world that simultaneously affect pricing decisions and demand are available in large data stores. One often appeals to modern statistical machine learning (ML) techniques to model and fit the high-dimensional or nonparametric nuisance parameters introduced by these confounders. However, most such techniques introduce bias into their estimates (e.g., via regularization) and hence yield invalid or inaccurate inferences concerning the parameters of interest (the treatment effects).
Several recent lines of have begun address the problem of debiasing ML estimators to perform accurate inference on a low dimensional component of model parameters. Prominent examples include Lasso debiasing ZhangZh2013,vanDeGeer2013,Javanmard2015 and post-selection inference Belloni2013,Berk2013,Tibhsirani2014. The recent double / debiased ML work of Chernozhukov2016 describes a general-purpose strategy for extracting valid inferences for target parameters from somewhat arbitrary and relatively inaccurate estimates of nuisance parameters.
Specifically, Chernozhukov2016 analyze a two-stage process where in the first stage one estimates nuisance parameters using arbitrary statistical ML techniques on a first stage data sample and in the second stage estimates the low dimensional parameters of interest via the generalized method of moments (GMM). Crucially, the moments in the second stage are required to satisfy a Neyman orthogonality condition, granting them first-order robustness to errors in the nuisance parameter estimation. A main conclusion is that the second stage estimates are $\sqrt{n}$-consistent and asymptotically normal whenever the first stage estimates are consistently estimated at a $o(n^{-1/4})$ rate.
To illustrate this result, let us consider the partially linear regression (PLR) model, popular in causal inference. In the PLR model we observe data triplets $Z = (T,Y,X)$, where $T \in \ensuremath{\mathbb R}$ represents a treatment or policy applied, $Y \in \ensuremath{\mathbb R}$ represents an outcome of interest, and $X \in \ensuremath{\mathbb R}^p$ is a vector of associated covariates. These observations are related via the equations
where $\eta$ and $\epsilon$ represent unobserved disturbances with distributions independent of $(\theta_0, f_0, g_0)$. The first equation features the treatment effect $\theta_0$, our object of inference. The second equation describes the relation between the treatment $T$ and the associated covariates $X$. The covariates $X$ affect the outcome $Y$ through the nuisance function $f_0$ and the treatment $T$ through the nuisance function $g_0$. Using the Neyman-orthogonal moment of Chernozhukov2016, the authors show that it suffices to estimate the nuisance $(f_0,g_0)$ at an $o(n^{-1/4})$ rate to construct a $\sqrt{n}$-consistent and asymptotically normal estimator of $\theta_0$.
In this work, we provide a framework for achieving stronger robustness to first stage errors while maintaining second stage validity. In particular, we introduce a notion of higher-order orthogonality and show that if the moment is $k$-th order orthogonal then a first-stage estimation rate of $o(n^{-1/(2k+2)})$ suffices for $\sqrt{n}$-asymptotic normality of the second stage.
We then provide a concrete application of our approach to the case of estimating treatment effects in the PLR model. Interestingly, we show an impossibility result when the treatment residual follows a Gaussian distribution: no higher-order orthogonal moments with finite asymptotic variance exist, so first-order Neyman orthogonality appears to be the limit of robustness to first stage errors under Gaussian treatment residual. However, conversely, we also show how to construct appropriate second-order orthogonal moments whenever the treatment residual is not Gaussian. As a result, when the nuisance functions are linear in the high-dimensional confounders, our second-order orthogonal moments provide valid inferences whenever the number of relevant confounders is $o(\frac{n^{2/3}}{\log p})$; meanwhile the first-order orthogonality analyses of Chernozhukov2016 accommodate only $o(\frac{\sqrt{n}}{\log p})$ relevant confounders.
We apply these techniques in the setting of demand estimation from pricing and purchase data, where highly non-Gaussian treatment residuals are standard. In this setting, the treatment is the price of a product, and commonly, conditional on all observable covariates, the treatment follows a discrete distribution representing random discounts offered to customers over a baseline price linear in the observables. In Figure (ref) we portray the results of a synthetic demand estimation problem with dense dependence on observables. Here, the standard orthogonal moment estimation has large bias, comparable to variance, while our second-order orthogonal moments lead to nearly unbiased estimation.
\paragraph{Notational conventions} For each $n\in \mathbb{N}$, we introduce the shorthand $[n]$ for $\{1,\dots,n\}$. We let $\stackrel{p}{\to}$ and $\stackrel{d}{\to}$ represent convergence in probability and convergence in distribution respectively. When random variables $A$ and $B$ are independent, we use $\mathbb{E}_A[g(A,B)] \triangleq \mathbb{E}[g(A,B) \mid B]$ to represent expectation only over the variable $A$. For a sequence of random vectors $(X_n)_{n=1}^\infty$ and a deterministic sequence of scalars $(a_n)_{n=1}^\infty$, we write $X_n = O_P(a_n)$ to mean $X_n/a_n$ is stochastically bounded, i.e., for any $\epsilon > 0$ there is $R_{\epsilon},N_{\epsilon}>0$ with $\Pr(\|{X_n/a_n}\| > R_{\epsilon}) \leq \epsilon$ for all $n >N_{\epsilon}$. We let $N(\mu, \Sigma)$ represent a multivariate Gaussian distribution with mean $\mu$ and covariance $\Sigma$.
Our aim is to estimate an unknown target parameter $\theta_0\in \Theta \subseteq \ensuremath{\mathbb R}^d$ given access to independent replicates $(Z_t)_{t=1}^{2n}$ of a random data vector $Z\in \ensuremath{\mathbb R}^{\rho}$ drawn from a distribution satisfying $d$ moment conditions,
Here, $X\in \ensuremath{\mathbb R}^{p}$ is a sub-vector of the observed data vector $Z$, $h_0 \in \mathcal{H} \subseteq \{h:\ensuremath{\mathbb R}^{p}\rightarrow \ensuremath{\mathbb R}^\ell\}$ is a vector of $\ell$ unknown nuisance functions, and $m: \ensuremath{\mathbb R}^{\rho} \times \ensuremath{\mathbb R}^d\times \ensuremath{\mathbb R}^\ell\rightarrow \ensuremath{\mathbb R}^d$ is a vector of $d$ known moment functions. We assume that these moment conditions exactly identify the parameter $\theta_0$, and we allow for the data to be high-dimensional, with $\rho$ and $p$ potentially growing with the sample size $n$. However, the number of parameters of interest $d$ and the number of nuisance functions $\ell$ are assumed to be constant.
We will analyze a two-stage estimation process where we first estimate the nuisance parameters using half of our sample\footnote{Unequal divisions of the sample can also be used; we focus on an equal division for simplicity of presentation.} and then form a $Z$-estimate of the target parameter $\theta_0$ using the remainder of the sample and our first-stage estimates of the nuisance. This sample-splitting procedure proceeds as follows.
Relegating only half of the sample to each stage represents a statistically inefficient use of data and, in many applications, detrimentally impacts the quality of the first-stage estimate $\hat{h}$. A form of repeated sample splitting called $K$-fold cross-fitting Chernozhukov2016 addresses both of these concerns. $K$-fold cross-fitting partitions the index set of the datapoints $[2n]$ into $K$ subsets $I_1,\ldots,I_K$ of cardinality $\frac{2n}{K}$ (assuming for simplicity that $K$ divides $2n$) and produces the following two-stage estimate:
Throughout, we assume $K$ is a constant independent of all problem dimensions. As we will see in Theorem {(ref)}, a chief advantage of cross-fitting over sample splitting is improved relative efficiency with an asymptotic variance that reflects the use of the full dataset in estimating $\theta$.
Main Question. Our primary inferential goal is to establish conditions under which the estimators ${\hat{\theta}^{SS}}$ in (ref) and ${\hat{\theta}^{CF}}$ (ref) enjoy $\sqrt{n}$-asymptotic normality, that is
for some constant covariance matrix $\Sigma$. Coupled with a consistent estimator of $\Sigma$, asymptotic normality enables the construction of asymptotically valid confidence intervals for $\theta$ based on Gaussian or Student's t quantiles and asymptotically valid hypothesis tests, like the Wald test, based on chi-squared limits.
We would like our two-stage procedures to produce accurate estimates of $\theta_0$ even when the first stage nuisance estimates are relatively inaccurate. With this goal in mind, Chernozhukov2016 defined the notion of Neyman-orthogonal moments, inspired by the early work of Neyman1979. In our setting, the orthogonality condition of Chernozhukov2016 is implied by the following condition, which we will call first-order orthogonality:
Intuitively, first-order orthogonal moments are insensitive to small perturbations in the nuisance parameters and hence robust to small errors in estimates of these parameters. A main result of Chernozhukov2016 is that, if the moments $m$ are first-order orthogonal, then $o(n^{-1/4})$ error rates\footnote{In the sense of root mean squared error: $n^{1/4} \sqrt{\mathbb{E}[\|{h_0(X)-\hat{h}(X)}\|_2^2\mid \hat{h}]}\stackrel{p}{\to} 0$.} in the first stage estimation of $h_0$ are sufficient for $\sqrt{n}$-asymptotic normality of the estimates ${\hat{\theta}^{SS}}$ and ${\hat{\theta}^{CF}}$.
Our aim is to accommodate slower rates of convergence in the first stage of estimation by designing moments robust to larger nuisance estimation errors. To achieve this, we will introduce a generalized notion of orthogonality that requires higher-order nuisance derivatives of $m$ to be conditionally mean zero. We will make use of the following higher-order differential notation:
We are now equipped to define our notion of $S$-orthogonal moments:
We will often be interested in the special case of Definition {(ref)} in which $S$ is comprised of all vectors $\alpha \in \ensuremath{\mathbb N}^{\ell}$ with $\|\alpha\|_1\leq k$. This implies that all mixed nuisance derivatives of the moment of order $k$ or less are conditionally mean zero. We will refer to this special case as $k$-orthogonality or $k$-th order orthogonality.
The general notion of $S$-orthogonality allows for our moments to be more robust to errors in some nuisance functions and less robust to errors in others. This is particularly valuable when some nuisance functions are easier to estimate than others; we will encounter such an example in Section {(ref)}.
We will now show that $S$-orthogonality together with appropriate consistency rates for the first stage estimates of the nuisance functions imply $\sqrt{n}$-consistency and asymptotic normality of the two-stage estimates ${\hat{\theta}^{SS}}$ and ${\hat{\theta}^{CF}}$. Beyond orthogonality and consistency, our main Assumption (ref) demands identifiability, non-degeneracy, and regularity of the moments $m$, all of which are standard for establishing the asymptotic normality of $Z$-estimators.
We are now ready to state our main theorem on the implications of $S$-orthogonality for second stage $\sqrt{n}$-asymptotic normality. The proof can be found in Section {(ref)}.
A variety of standard sufficient conditions guarantee the consistency of ${\hat{\theta}^{SS}}$ and ${\hat{\theta}^{CF}}$. Our next result, proved in Section {(ref)}, establishes consistency under either of two commonly satisfied assumptions.
Our assumption on the first stage estimation rates, i.e., that $\forall \alpha \in \{ a \in \mathbb{N}^{\ell}: \|a\|_1 \leq k+1 \} \setminus S$ \balignst n^{1/2}\cdot \sqrt{\mathbb{E}\left[\prod_{i=1}^\ell |\hat{h}_i(X)-h_{0,i}(X)|^{2\alpha_i} \mid \hat{h} \right]} \stackrel{p}{\to} 0 \ealignst may seem complex, as it involves the interaction of the errors of multiple nuisance function estimates. In this section we give sufficient conditions that involve only the rates of individual nuisance function estimates and which imply our first stage rate assumptions. In particular, we are interested in formulating consistency rate conditions for each nuisance function $h_i$ with respect to an $\mathcal{L}^p$ norm,
We will make use of these sufficient conditions when applying our main theorem to the partially linear regression model in Section {(ref)}.
A simpler description of the sufficient conditions arises under $k$-orthogonality (Definition {(ref)}), since the set $\{ a \in \mathbb{N}^{\ell}: \|a\|_1 \leq k+1 \} \setminus S_k$ contains only vectors $\alpha$ with $\|\alpha\|=k+1$.
In the case of first-order orthogonality, Corollary {(ref)} requires that the first stage nuisance functions be estimated at a $o(n^{-1/4})$ rate with respect to the $\mathcal{L}^{4}$ norm. This is almost but not exactly the same as the condition presented in Chernozhukov2016, which require $o(n^{-1/4})$ consistency rates with respect to the $\mathcal{L}^{2}$ norm. Ignoring the expectation over $X$, the two conditions are equivalent.\footnote{We would recover the exact condition in Chernozhukov2016 if we replaced Assumption {(ref).(ref)} with the more stringent assumption that $\left|D^{\alpha}m(Z, \theta, h(X))\right| \leq \lambda_*$ a.s.} Moreover, in the case of $k$-orthogonality, Corollary {(ref)} requires $o(n^{-1/2(k+1)})$ rates with respect to the $\mathcal{L}^{2(k+1)}$ norm. More generally, $S$-orthogonality allows for some functions to be estimated slower than others as we will see in the case of the sparse linear model.
When second-order orthogonal moments satisfying Assumption (ref) are employed, Corollary {(ref)} implies that an $o(n^{-1/6})$ rate of nuisance parameter estimation is sufficient for $\sqrt{n}$-consistency of ${\hat{\theta}^{SS}}$ and ${\hat{\theta}^{CF}}$. This asymptotic improvement over first-order orthogonality holds the promise of accommodating more complex and higher-dimensional nuisance parameters. In this section, we detail both the limitations and the power of this approach in the partially linear regression (PLR) model setting popular in causal inference Chernozhukov2016.
Our first result shows that, under the PLR model, if the treatment noise, $\eta$, is conditionally Gaussian given $X$, then no second-order orthogonal moment can satisfy Assumption (ref), because every twice continuously differentiable $2$-orthogonal moment has $\mathbb{E}\left[\nabla_\theta m(Z, \theta_0, h_0(X))\right] = 0$ (a violation of Assumption (ref).3). The proof in Section {(ref)} relies on Stein's lemma.
In the following result, proved in Section {(ref)}, we establish that under mild conditions Assumption {(ref).(ref)} is necessary for the $\sqrt{n}$-consistency of ${\hat{\theta}^{SS}}$ in the PLR model.
We next show that, inversely, second-order orthogonal moments are available whenever the conditional distribution of treatment noise given $X$ is not a.s.\ Gaussian. Our proofs rely on a standard characterization of a Gaussian distribution, proved in Section {(ref)}:
We will focus on estimating the nuisance functions $q_0=f_0+\theta_0g_0$ and $g_0$ instead of the nuisance functions $f_0$ and $g_0$, since the former task is more practical in many applications. This is because estimating $q_0$ can be accomplished by carrying out an arbitrary non-parametric regression of $Y$ onto $X$. In contrast, estimating $f_0$ typically involves regressing $Y$ onto $(X,T)$, where $T$ is constrained to enter linearly. The latter might be cumbersome when using arbitrary ML regression procedures.
Our first result, established in Section {(ref)}, produces finite-variance 2-orthogonal moments when an appropriate moment of the treatment noise $\eta$ is known.
Our next result, proved in Section {(ref)}, addresses the more realistic setting in which we do not have exact knowledge of $\mathbb{E} \left[ \eta^{r} |X\right]$. We introduce an additional nuisance parameter and still satisfy an orthogonality condition with respect to these parameters.
We now consider deploying the PLR model in the high-dimensional linear regression setting, where $f_0(X) = \langle{X},{\beta_0}\rangle$ and $g_0(X) = \langle{X},{\gamma_0}\rangle$ for two $s$-sparse vectors $\beta_0, \gamma_0 \in \mathbb{R}^p$, $p$ tends to infinity as $n\to \infty$, and $(\eta, \epsilon, X)$ are mutually independent. Define $q_0 = \theta_0 \beta_0 + \gamma_0$. In this high-dimensional regression setting, Chernozhukov2016 showed that two-stage estimation with first-order orthogonal moments \balignt &m\left(Z, \theta, \langle{X},{q}\rangle,\langle{X},{\gamma}\rangle\right)=\\ \notag &\left(Y-\langle{X},{q}\rangle-\theta \left(T-\langle{X},{\gamma}\rangle\right) \right)\left(T-\langle{X},{\gamma}\rangle\right) \ealignt and Lasso estimates of the nuisance provides a $\sqrt{n}$-asymptotically normal estimator of $\theta_0$ when $s = o({n^{\frac{1}{2}}}{/\log p})$. Our next result, established in Appendix {(ref)}, shows that we can accommodate $s = o({n^{\frac{2}{3}}}{/\log p})$ with an explicit set of higher-order orthogonal moments.
We perform an experimental analysis of the second order orthogonal estimator of Theorem {(ref)} with $r=3$ for the case of estimating treatment effects in the PLR model with high-dimensional sparse linear nuisance functions. We compare our estimator with the double ML estimator (labeled `dml' in our figures) based on the first-order orthogonal moments (ref) of Chernozhukov2016. Our experiments are designed to simulate demand estimation from pricing and purchase data, where non-Gaussian treatment residuals are standard. Here, our covariates $X$ correspond to all collected variables that may affect a pricing policy. A typical randomized experiment in a pricing policy takes the form of random discounts from a baseline price as a company offers random discounts to customers periodically to gauge demand level. In this case, the treatment residual -- the unexplained fluctuation in price -- is decidedly non-Gaussian and specifically follows a discrete distribution over a small number of price points. Python code recreating all experiments is available at \url{https://github.com/IliasZadik/double_orthogonal_ml}.
\paragraph{Experiment Specification} We generated $n$ independent replicates of outcome $Y$, treatment $T$, and confounding covariates $X$. The confounders $X$ have dimension $p$ and have independent components from the $N(0,1)$ distribution. The treatment is a sparse linear function of $X$, $T= \ensuremath{\langle \gamma_0, X \rangle} + \eta$, where only $s$ of the $p$ coefficients of $\gamma_0$ are non-zero. The $x$-axis on each plot is the number of non-zero coefficients $s$. Moreover, $\eta$ is drawn from a discrete distribution with values $\{0.5, 0, -1.5, -3.5\}$ taken respectively with probabilities $(.65, .2, .1, .05)$. Here, the treatment represents the price of a product or service, and this data generating process simulates random discounts over a baseline price. Finally, the outcome is generated by a linear model, $Y = \theta_0 T + \ensuremath{\langle \beta_0, X \rangle} + \epsilon$, where $\theta_0=3$ is the treatment effect, $\beta_0$ is another sparse vector with only $s$ non-zero entries, and $\epsilon$ is drawn independently from a uniform $U(-\sigma_\epsilon, \sigma_\epsilon)$ distribution. Importantly, the coordinates of the $s$ non-zero entries of the coefficient $\beta_0$ are the same as the coordinates of the $s$ non-zero entries of $\gamma_0$. The latter ensures that variables $X$ create a true endogeneity problem, i.e., that $X$ affects both the treatment and the outcome directly. In such settings, controlling for $X$ is important for unbiased estimation.
To generate an instance of the problem, the common support of both $\gamma_0$ and $\beta_0$ was generated uniformly at random from the set of all coordinates, and each non-zero coefficient was generated independently from a uniform $U(0,5)$ distribution. The first stage nuisance functions were fitted for both methods by running the Lasso on a subsample of $n/2$ sample points. For the first-order method all remaining $n/2$ points were used for the second stage estimation of $\theta_0$. For the second-order method, the moments $\mathbb{E}[\eta^2]$ and $\mathbb{E}[\eta^3]$ were estimated using a subsample of $n/4$ points as described in Theorem {(ref)}, and the remaining $n/4$ sample points were used for the second stage estimation of $\theta_0$. For each method we performed cross-fitting across the first and second stages, and for the second-order method we performed nested cross-fitting between the $n/4$ subsample used for the $\mathbb{E}[\eta^2]$ and $\mathbb{E}[\eta^3]$ estimation and the $n/4$ subsample used for the second stage estimation. The regularization parameter $\lambda_n$ of each Lasso was chosen to be $\sqrt{\log(p)/n}$.
For each instance of the problem, i.e., each random realization of the coefficients, we generated $2000$ independent datasets to estimate the bias and standard deviation of each estimator. We repeated this process over $100$ randomly generated problem instances, each time with a different draw of the coefficients $\gamma_0$ and $\beta_0$, to evaluate variability across different realizations of the nuisance functions.
\paragraph{Distribution of Errors with Fixed Sparsity} In Figure (ref), we display the distribution of estimates based on orthogonal moments and second-order orthogonal moments for a particular sparsity level $s=100$ and for $n=5000$ and $p=1000$. We observe that both estimates are approximately normally distributed, but the orthogonal moment estimation exhibits significant bias, an order of magnitude larger than the variance.
\paragraph{Bias-Variance Tradeoff with Varying Sparsity} Figure (ref) portrays the median quantities (solid lines) and maximum and minimum of these quantities (error bars) across the $100$ different nuisance function draws as a function of the support size for $n=5000$, $p=1000$, and $\sigma_\epsilon=1$.
\paragraph{Varying $n$ and $p$} In Figure {(ref)}, we display how performance varies with $n$ and $p$. Due to computational considerations, for this parameter exploration, we only used a single problem instance for each $(n, p, s)$ triplet rather than $100$ instances as in the exploration above. We note that for $n=2000, p=5000$ the breaking point of our method is around $s=100$, while for $n=5000, p=2000$ it is around $s=550$. For $n=10000, p=1000$ our method performs exceptionally well even until $s=800$.
\paragraph{Varying $\sigma_\epsilon$} Finally Figure {(ref)} displays performance as the variance $\sigma_\epsilon$ of the noise $\epsilon$ grows.
Our aim in this work was to conduct accurate inference for fixed-dimensional target parameters in the presence of high-dimensional or nonparametric nuisance. To achieve this, we introduced a notion of $k$-th order orthogonal moments for two-stage $Z$-estimation, generalizing the first-order Neyman orthogonality studied in Chernozhukov2016. Given $k$-th order orthogonal moments, we established that estimating nuisance at an $o(n^{-1/(2k+2)})$ rate suffices for $\sqrt{n}$-consistent and asymptotically normal estimates of target parameters. We then studied the PLR model popular in causal inference and showed that a valid second-order orthogonal moment exists if and only if the treatment residual is not normally distributed. In the high-dimensional linear nuisance setting, these explicit second-order orthogonal moments tolerate significantly denser nuisance vectors than those accommodated by Chernozhukov2016. We complemented our results with synthetic demand estimation experiments showing the benefits of second-order orthogonal moments over standard Neyman-orthogonal moments.