Automatic Debiased Machine Learning via Riesz Regression

Automatic Debiased Machine Learning Via Riesz Regression

Introduction

Many parameters of interest depend on regressions. Examples include treatment effects, regression decompositions, and policy effects. Often, a regression may be high dimensional, depending on many variables. For example there may be many covariates for treatment effects. Machine learning methods such as neural nets, random forests, and Lasso can be used to estimate parameters of interest that depend on high dimensional regressions.

A general problem with estimating parameters of interest using machine learning is that machine learners are biased by regularization and/or model selection. This bias may pass through when the learner is plugged into a formula for a parameter of interest and make the parameter estimator highly biased. This problem can be avoided by using Neyman orthogonal estimating equations where machine learners have zero first-order effect. Cross-fitting, a form of sample splitting, can also help.

The orthogonal estimating equations for regressions depend on a Riesz representer $\alpha_0$ that must be estimated. The primary innovation of this paper is to provide an automatic estimator of $\alpha_0$ that uses only the definition of the parameter of interest and the regression but does not require knowing a formula for $\alpha_0$. We give an objective function with expectation that is minimized at $\alpha_0$ that depends only the parameter of interest. We refer to minimization of this objective function as a Riesz regression, being equivalent to minimizing the expected squared deviation from $\alpha_0$. Neural nets, random forests, and other methods can be used for this Riesz regression. Using the Riesz regression estimator in the bias correction completes an algorithm that 1) specifies the parameter of interest; 2) specifies a learner of the unknown regression; and 3) uses the Riesz regression estimator of $\alpha_0$ determined by steps 1) and 2).

A second innovation of this paper is to construct and derive properties of estimators that depend on generalized regressions, which minimize an expected loss over some linear set of functions. These generalized regressions include conditional means, least squares projections, functions that minimize quasi-likelihoods, and quantile regressions. Debiasing for generalized regressions depends on a weighted version of the Riesz representer. We give a weighted Riesz regression that only uses the parameter of interest and the generalized regression for bias correction.

A third contribution of this paper is finite sample mean square error bounds for Reisz regressions. These bounds are obtained using the critical radius of functions of $\alpha$ on which the objective function depends and approximation error bounds for the unknown $\alpha_0$. A fourth contribution is convergence rates for neural net Riesz regressions. These are based on known results on critical radius and approximation error for neural nets and the finite sample bounds given here.

In work that followed up on the first version of this paper (chernozhukov2022riesznet) we found that using Riesz regressions to debias neural net and random forest estimators of the average treatment effect was much more accurate than state of the art methods based on inverse propensity score weighting, in a Monte Carlo study. Both the automatic neural net and random forest debiasing also led to accurate confidence intervals in those experiments.

Automatic debiasing for Lasso and reproducing kernel Hilbert space regressions was previously given by chernozhukov2022automatic and singh2022kernel respectively. The estimator of $\alpha_0$ given here goes beyond these to provide automatically debiasing for generalized regressions based on neural nets, random forests, and other machine learners. These innovations allow researchers to use any of a wide variety of automatically debiased machine learners learners to estimate parameters of interest that depend on generalized regressions. For example, automatic debiased machine learning with neural nets could be especially useful for parameters that depend on high dimensional, nonlinear generalized regressions.

This paper builds on recent work on Neyman orthogonal scores and debiased machine learning. We use model free orthogonal estimating equations like those of chernozhukov2022locally that are the sum of an identifying moment function and a bias adjustment (influence function) term for generalized regression from ichimura2022influence. Those papers did not give the Riesz regression. Finite sample mean square error bounds for a general learner of $\alpha_0$ are obtained by applying results of foster2019orthogonal that characterize error bounds in terms of critical radius and approximation. The rate of convergence for neural net Riesz regression use critical radius and approximation rate results given in farrell2021deep. Additional neural net rate conditions could be obtained using yarotsky2018optimal. The learner of $\alpha_0$ differs from those of farrell2021deep, farrell2021deep2 in using the Riesz regression rather than a known form for $\alpha_0$.

We also build upon ideas in classical semi- and nonparametric learning theory with low dimensional regressions using traditional smoothing methods (van1991differentiable, bickel1993efficient, newey1994asymptotic, robins1995semiparametric, van2000asymptotic), that do not apply to machine learners. The orthogonal estimating equations given in chernozhukov2022locally and used here build on previous work on nonparametric orthogonal moment functions by levit1975efficiency,hasminskii1978nonparametric,bickel1988estimating,newey2004twicing. Targeted maximum likelihood (van2006targeted) based on machine learners has been considered by van2011targetedbook and large sample theory given by luedtke2016optimal.

In section (ref) we give the Reisz regression and an automatic debiased machine learning algorithm for parameters that depend linearly on a nonparametric regression, including examples. Section (ref) also gives finite sample mean square error bounds for the general and neural net Riesz regression and asymptotic inference results for parameters that are linear functionals of a nonparametric regression. Section (ref) extends the estimation methods and theory to nonlinear functionals of generalized regressions. In section (ref) we illustrate the usefulness of our methods with an empirical application. Section (ref) presents the results of our simulation exercises.

Average Linear Effects for a Conditional Mean

To highlight the innovation provided by the Riesz regression, we first consider average linear effects that depend on a conditional mean. In section (ref) we consider the general setting of nonlinear functions of generalized regressions.

Parameters of Interest

We consider data that consists of i.i.d. observations $W_1, \ldots, \allowbreak W_n$, each having CDF $F_0$. A data observation $W$ includes an outcome variable $Y$ and regressors $X$. In this section, we focus on parameters that depend on the conditional mean of $Y$ given $X$. We will denote a possible such regression function by $\gamma$, with $\gamma_0(x) = \operatorname{E}[Y \mid X = x]$ being the true regression function.

The parameter of interest $\theta_0$ has the form

where $m(w, \gamma)$ is a functional that depends on a data observation $w$ and a possible regression function $\gamma$. For now, we assume that the expected functional $\gamma \mapsto \operatorname{E}[m(W, \gamma)]$ is linear and continuous in $\gamma$, meaning that there is a constant $C>0$ with $|\operatorname{E}[m(W, \gamma)]|^{2}\leq C\operatorname{E}[\gamma(X)^2]$ for all $\gamma$ with $\operatorname{E}[\gamma(X)^2] < \infty$. Under this assumption, there exists a function $v_m$ with $\operatorname{E}[v_m(X)^2] < \infty$ such that

The existence of this $v_m$ follows from the Riesz representation theorem, and it is equivalent to the semiparametric variance bound for $\theta_0$ being finite (see newey1994asymptotic, hirshberg2021augmented, chernozhukov2019global). For these reasons $v_m$ is often referred to as the Riesz representer. In this Section, where the parameter of interest depends on a nonparametric regression, the Riesz representer $v_m = \alpha_0$ needs to be estimated for debiased machine learning. In section (ref), where $\gamma$ may be a generalized regression, $\alpha_0$ will be a weighted version of the Riesz regression.

There are a variety of important, empirically relevant parameters of interest that have this form. We illustrate with some familiar examples to help highlight and motivate the Riesz regression:

Estimation

We will base estimation of $\theta_0$ on a Neyman orthogonal estimating equation, i.e. score, where first step estimation has zero first order effects, that is also doubly robust in having expectation zero if either $\gamma=\gamma_0$ or $\alpha=\alpha_0$. This score is

as in chernozhukov2022locally, where taking expectations gives, for any $\alpha, \gamma$,

Here we see that at the true parameter value $\theta_0$, the expectation of the score $\psi(W,\gamma,\alpha,\theta_0)$ differs from zero only to second order and equals zero if either $\gamma=\gamma_0$ or $\alpha=\alpha_0$. Thus the score is Neyman orthogonal and doubly robust in that it has zero expectation if either $\gamma=\gamma_0$ or $\alpha=\alpha_0$.

Estimation of $\alpha_0$ is essential to construction of a debiased machine learner of the parameter of interest. The primary innovation of this paper is to give an extremum characterization of $\alpha_0$ and use this to estimate $\alpha_0$. This extremum characterization is given by

where the third equality holds because $\operatorname{E}[\alpha_0(X)^2]$ does not depend on $\alpha$ and $\alpha_0=v_m$, and the fourth equality follows from equation (ref) with $\gamma=\alpha$. This characterization can be used to estimate $\alpha_0$ by replacing the expectation with the sample average and minimizing over some set of possible $\alpha$ functions.

We call a resulting estimator of $\alpha_0$ a Riesz regression, motivated by minimization of the least squares objective function in equation (ref). This estimator is automatic in dependng only on the function $m(w,\gamma)$ that determines the parameter of interest and in not requiring the form of $\alpha_{0}$. In particular, this method does not depend on plugging in non-parametric estimates of components of $\alpha_{0}$. This feature is useful when $\alpha_0$ does not have a simple form. For causal parameters such as those of Examples 1-3, the Riesz regression avoids inverting a learner of a conditional probability or a pdf. Instead, the Riesz regression learns $\alpha_{0}$ directly.

Our estimation strategy for the parameter of interest is to combine the Reisz regression estimator of $\alpha_0$ and estimation of $\gamma_0$ in the Neyman orthogonal score with the use of cross-fitting to reduce overfitting bias.\footnote{See newey2018cross for more on the advantages of cross-fitting.} The outline of our estimation strategy is as follows:

Our estimation strategy is very general, allowing for any choice of learner $\hat{\gamma}_{\ell}$ and any Riesz regression $\hat{\alpha}_{\ell}$ encoded in the class of functions $\mathcal{A}_n$. Special kinds of Riesz regressions have been given in previous literature. These include linear combinations of a dictionary of functions $(b_1(x), \ldots, b_p(x))'$, and $p$ large, with an $L_1$ penalty in the loss function (chernozhukov2022automatic), or functions embedded in a reproducing kernel Hilbert space (singh2021debiased). chernozhukov2020adversarial allowed for the estimation of $\alpha_0$ in arbitrary function spaces, but proposed a computationally harder minimax loss formulation. A primary innovation of this paper is to provide the Riesz regression for automatic estimation of $\alpha_0$ and corresponding asymptotic theory.

As an example, below we will give primitive conditions for a neural net Riesz regression. A general neural net takes the form \[x\overset{f_{1}}{\longmapsto}H^{(1)}\overset{f_{2}}{\longmapsto}\cdots\overset{f_{m}}{\longmapsto}H^{(m)}\] where $H^{(l)}=\{H_{k}^{(l)}\}_{k=1}^{K_{l}}$ are called neurons, $x$ is the original finite-dimensional input, and the function $f_{l}$ maps one layer of neurons to the next as in \[f_{l}:v\overset{f_{l}}{\longmapsto}\{H_{k}^{(l)}(v)\}_{k=1}^{K_{l} }:=(1,\{\sigma(v^{\prime}\beta_{k,l})\}_{k=2}^{K_{l}}),\] where each $\beta_{k,l}$ is a $K_{l-1}$ vector of parameters and $\sigma(u)$ is a nonlinear activation function. We will focus on the case where $\sigma(u)$ is the RELU function $\sigma(u)=\max\{0,u\}$. An important special case is a multilayer perceptron (MLP) network where the number of neurons $K_{l}=K$ is the same for each layer, for which results were recently given by farrell2021deep. Sparse versions of this specification, where many of the elements of the coefficient vectors $\beta_{k,l}$ may be zero, have also been considered recently by schmidt2020nonparametric. yarotsky2018optimal gave other neural net specifications with good approximation properties.

In the setting of Example (ref), a neural net Riesz regression would be constructed as $\hat\alpha_\ell(d, z) = \alpha(d, z; \hat\beta_\ell)$, where: \[\hat\beta_\ell = \arg\min_{\beta} \bigg[ \sum_{i \notin I_\ell} \left\lbrace - 2[\alpha(1, Z_i; \beta) - \alpha(0, Z_i; \beta) ] + \alpha(D_i, Z_i; \beta)^2 \right\rbrace + \Lambda_r(\beta) \bigg] \] for some penalty function $\Lambda_r(\beta)$ (e.g., L1, L2, or the elastic net). Because $D$ is binary, a convenient neural net architecture in this case could be a bi-headed MLP, $\alpha(d, z; \beta, \delta_0, \delta_1) = d g(z; \beta)'\delta_1 + (1 - d) g(z; \beta)'\delta_0$, where $g(z; \beta)$ is an MLP. An even more flexible specification would be to have $\alpha(d, z; \beta_0, \beta_1) = d g(z, \beta_0) + (1 - d) g(z, \beta_1)$, i.e. an MLP for the case $d = 1$ and another MLP for the case $d = 0$.

For Example (ref), a neural net Riesz regression is $\hat\alpha_\ell(d, z) = \alpha(d, z; \hat\beta_\ell)$, where: \[\hat\beta_\ell = \arg\min_{\beta} \bigg[ \sum_{i \notin I_\ell} \left\lbrace - 2[\partial_d \alpha(D_i, Z_i; \beta)] + \alpha(D_i, Z_i; \beta)^2 \right\rbrace + \Lambda_r(\beta) \bigg]. \] In particular, notice that the loss function involves taking a derivative of the neural net with respect to one of the inputs. A convenient parametrization of the neural net in this case is a locally linear function $\alpha(d, z; \phi, \beta) = \phi(d, z)'g(z; \beta)$, where $\phi(d, z)$ is a dictionary of known, differentiable basis functions, and $g(z; \beta)$ is a neural net (e.g. an MLP). In that case, $\partial_d \alpha(d, z; \beta) = [\partial_d \phi(d, x)]'g(z; \beta)$. This approach was used in chernozhukov2022riesznet to construct a random forest estimator of $\alpha_0$, exhibiting good performance in Monte Carlo simulations.

Large Sample Inference for Linear Effects of Regression

In this Section, we give mean square convergence rates for learners $\hat{\alpha}_{\ell}$ and $\sqrt{n}$-consistency and asymptotic normality results for the learner $\hat{\theta}$ of the object of interest and its asymptotic variance estimator $\hat{V}$. We first derive convergence rates for $\hat{\alpha}_{\ell}.$

Convergence Rates for $\hat{\alpha}_{\ell}$

In this subsection we suppress the $\ell$ subscript for notational convenience. We consider the problem of estimating \[ \alpha_{0}=\arg\min_{\alpha}\operatorname{E}[-2m(W,\alpha)+\alpha(X)^{2}], \] where we have used the extremum characterization of $\alpha_0$ in equation (ref). For any random variable $a(W)$ let $\left\Vert a\right\Vert =\sqrt{\operatorname{E}[a(W)^{2}]}$ and $\left\Vert a\right\Vert _{\infty}=\sup _{w\in\mathcal{W}}\left\vert a(w)\right\vert.$ \vsedit{For simplicity of exposition we will only consider the case where the estimate is defined over a growing sieve space $\mathcal{A}_n$ and no regularization is used, i.e. $\Lambda_r(\alpha)=0$:

Our estimation rate can easily be extended to regularized estimation with appropriate regularization weight.} We assume that $m(W,\alpha)$ is mean square continuous in the following sense:

Define:

Define: \[ \alpha^{\ast}=\arg\min_{\alpha\in\mathcal{A}_n}\operatorname{E}[-2m(W,\alpha)+\alpha(X)^{2}] \] to be the best approximation of $\alpha_{0}$ by an element of $\mathcal{A}_n$.

See e.g. foster2019orthogonal for the definition of the critical radius used in the statement this result. To use Theorem (ref) to obtain a mean square convergence rate for $\hat{\alpha}$ it is important to know the critical radius and the rate at which $\Vert\alpha_{\ast}-\alpha_{0}\Vert$ shrinks as the approximating set $\mathcal{A}_n$ becomes richer. For example, farrell2021deep have recently obtained such results for deep, ReLU neural nets. We can apply their results to obtain a mean square rate for such a learner of $\alpha_{0}$ when $x$ is a $d$ dimensional input for the multilayer perceptron (MLP) network with $m$ layers and width $K.$

The convergence rate depends on the smoothness of the function $\alpha_{0}(x),$ as specified in the following result. Specifically we assume that the support of $X$ is contained in a Cartesian product $\mathcal{X}$ of compact intervals and $\alpha_{0}(X)$ can be extended to a function that is continuously differentiable on $\mathcal{X}$ and has $\beta$ continuous derivatives.

When $\alpha_{0}$ is smooth enough, in that $\beta$ is large enough, the upper bound on $\left\Vert \hat{\alpha} -\alpha_{0}\right\Vert $ in Corollary (ref) gives a mean square convergence rate that can be close to, but less than $n^{-1/2}$. Such rate can be obtained by choosing the width $K$ and depth $m$ to approximately balance the two terms in Corollary (ref), $K \asymp n^{\frac{d}{2(\beta + d)}}\ln^2(n)$, $m \asymp \ln(n)$, as in farrell2021deep, in which case \[ \left\Vert \hat{\alpha}-\alpha_{0}\right\Vert^{2} = O_p\left(n^{-\frac{\beta}{\beta + d} }\ln^8(n)\right). \] Faster rates could be obtained using the neural nets of yarotsky2018optimal or the sparse neural nets of schmidt2020nonparametric. We focus on Corollary (ref) for the MLP neural net because it is a widely used architecture in practice, and because the rates obtained are fast enough for the estimators of the parameter of interest to be asymptotically normal.

Large Sample Inference for $\theta_{0}$

We use additional regularity conditions to show asymptotic normality of $\hat{\theta}$ and consistency of the asymptotic variance estimator $\hat{V}$. We will first give a general result for $\hat{\theta}$ that applies to any $\hat{\alpha}_{\ell}$ and does not rely on Theorem (ref) for a convergence rate for $\hat{\alpha}_{\ell}$. Similarly, any regression learner $\hat{\gamma}_\ell$ can be used here as long as its mean-square convergence rate is fast enough, as formalized below. Such convergence rate results are available for shallow (chen1999improved) and deep (yarotsky2018optimal, schmidt2020nonparametric, farrell2021deep) neural nets, random forests (syrgkanis2020estimation), LASSO (bickel2009simultaneous), boosting (luo2022highdimensional) and other high-dimensional methods.

The following assumption imposes a few additional regularity conditions. Let $\sigma_0^2(X) = \operatorname{E}[(Y - \gamma_0(X))^2 \mid X]$ denote the conditional variance of $Y$ given $X$.

Next, we require mean square consistency of $\hat{\gamma}_{\ell}$ and $\hat{\alpha}_{\ell},$ that the product of their mean-square convergence rates is smaller than $n^{-1/2},$ and a boundedness condition for $\hat{\alpha}_{\ell}$.

Part (i) implies that both $\hat{\gamma}_\ell$ and $\hat{\alpha}_\ell$ are consistent in mean square. Part (ii) captures an important tradeoff between the rates of convergence for $\hat{\gamma}_\ell$ and $\hat{\alpha}_\ell$. In settings where the regression can be estimated at a relatively fast rate of convergence, the learner for the debiasing function can converge more slowly, and vice versa, as long as the product of their mean-square convergence rates vanishes faster than $n^{-1/2}$. The results we have obtained for the neural net learner $\hat{\alpha}_{\ell}$ can be used to verify these conditions and we do so in Corollary (ref) to follow. The mean square convergence of $\hat{\gamma}_{\ell}$ is a primitive condition for this paper and allows use of a wide variety of $\hat{\gamma}_{\ell}$ in the construction of the estimator.

We have the following large sample inference result under these conditions.

Next we use Theorem (ref) and Corollary (ref) to formulate regularity conditions when $\hat{\alpha}_{\ell}$ is the neural net learner of $\alpha_{0}$ in section (ref). Let \[ \epsilon_{\alpha n}^2=K^{2}m^{2}\ln(K^{2}m)\ln(n)/n+[Km\sqrt{\ln(K^{2} m)}]^{-2(\beta/d)+\varepsilon}. \] This $\epsilon_{\alpha n}^2$ is taken from the upper bound for $\left\Vert\hat{\alpha}_{\ell}-\alpha_{0}\right\Vert ^{2}$ in Corollary (ref) and so characterizes the mean square convergence rate of the automatic neural net learner $\hat{\alpha}_{\ell}$.

Average Effects for Generalized Regressions

Linear Effects

In this section we extend the results to parameters that depend on functions $\gamma_{0}$ other than the conditional mean, that we refer to as generalized regressions. Suppose that $\gamma_0$ is defined as the solution to a general $M$-estimation problem:

where $\Gamma$ is a closed (in mean square) linear subspace of $L^2(X)$. For example, when $\ell(W, \gamma) = \frac{1}{2} (Y - \gamma(X))^2$ is the square loss and $\Gamma = L^2(X)$, then $\gamma_0(X) = \operatorname{E}[Y \mid X]$ and we recover the case of regression.

By the first order condition of the minimization problem (ref), $\gamma_0$ satisfies

for some functional $\rho(W, \gamma)$, typically a generalized notion of the (negative) derivative of the loss function $\ell(W, \gamma)$. In the case of regression, we can take $\rho(W, \gamma) = Y - \gamma(X)$ to be the non-parametric residual. For other statistical problems, we will refer to the function $\rho(W, \gamma)$ as a generalized residual. The results of this section will apply to any $\gamma_0$ that is identified by an orthogonality condition as in (ref), even beyond $M$-estimation problems.

This setting covers many interesting features of the conditional distribution of $Y$ given $X$. First, suppose that $\Gamma = L^2(X)$, so that the functional form of $\gamma_0$ is unrestricted. For example, when $\rho(W ,\gamma) = \tau - 1(Y < \gamma(X))$ for $0 < \tau < 1$, then $\gamma_0(x)$ is the $\tau$-th conditional quantile of $Y$ given $X = x$. When $\rho(W,\gamma)=\lambda(\gamma(X))[Y-\mu(\gamma(X))]$ for a link function $\mu(a)$ and another function $\lambda(a)$, this corresponds to the first order conditions of a generalized linear model nelder1972generalized. For binary $Y \in \{0, 1\}$, $\mu(a)$ the standard logistic CDF, and $\lambda(a) \equiv 1$, for instance, this set up corresponds to a (non-parametric) logistic regression, where $\gamma_0(X) = \mu^{-1}(\Pr(Y = 1 \mid X)) = \ln (\Pr(Y = 1 \mid X)/ \Pr(Y = 0 \mid X))$ corresponds to the log-odds.

The set $\Gamma$ could be used to encode parametric or semi-parametric restrictions on $\gamma_0$. One example is $X=(X_{1},X_{2},...)$ and $\Gamma$ the mean square closure of finite linear combinations of $X$. This corresponds to a high (infinite) dimensional, approximately sparse $\Gamma$, where the orthogonality condition is equivalent to $\operatorname{E}[X_{j} \rho(W,\gamma_{0})]=0$ for all $j$. In a case considered also by hirshberg2021augmented and farrell2021deep, $\Gamma$ is the mean square closure of $\{a(X_{1})+X_{2}^{\prime}b(X_{1})\}$ where $a(X_{1})$ is a scalar function and $b(X_{1})$ a vector of functions, each having unrestricted functional form. We could also take $\Gamma$ to be the mean-square closure of additive functions $a_{1}(X_{1})+a_{2}(X_{2})$, where $X_{1}\,\ $and $X_{2}$ are distinct components of $X$. In some cases, the resulting $\gamma_0$ will have a projection interpretation: for instance, when $\ell(W, \gamma)$ is the square loss, then $\gamma_0 = \arg\min_{\gamma \in \Gamma} \operatorname{E}[(\operatorname{E}[Y \mid X] - \gamma(X))^2]$ is the best approximation to $\operatorname{E}[Y \mid X]$ in $\Gamma$ in the mean-square sense.

For now, we also continue to assume that the parameter of interest has the form $\theta_{0} = \operatorname{E}[m(W,\gamma_{0})]$, where $\gamma\mapsto m(W,\gamma)$ is linear and $\operatorname{E}[m(W,\gamma)]$ is mean square continuous on $\Gamma$. We will relax the linearity assumption in section (ref). We will extend the results of the previous section by modifying (ref). Define, for any $\gamma, \alpha \in \Gamma$, a score:

where we have replaced $y - \gamma(x)$ with the generalized residual $\rho(W,\gamma)$.

This score satisfies $\operatorname{E}[\psi(W, \gamma_0, \alpha, \theta_0)] = \operatorname{E}[\alpha(X) \rho(W, \gamma_0)] = 0$ for any $\alpha \in \Gamma$ by (ref), and hence it is Neyman-orthogonal with respect to $\alpha$. To get Neyman orthogonality with respect to $\gamma$ we need to find a function $\alpha_0$ that satisfies

where $v_m(X)$ is the Riesz representer in (ref) and, for a scalar $a$,

that we assume exists. We further assume that we can normalize the sign of $\rho(W,\gamma)$ so that $v_{\rho}(X)\leq0$, as will hold when $\operatorname{E}[\rho(W, \gamma_{0} + a) \mid X]$ is monotonically decreasing in $a$. For example, when $\rho(W,\gamma)=Y-\gamma(X)$ as in Section 2 we have $v_{\rho}(X)=-1$. Also, when $\rho(W,\gamma)=p-1(Y<\gamma(X))$ then $v_{\rho}(X)=-f_{Y\mid X} (\gamma_{0}(X)\mid X),$ the negative of the conditional pdf of $Y$ given $X$ evaluated at $y=\gamma_{0}(X)$.\footnote{Note that since $\gamma_0$ corresponds to the $p$ quantile of $Y\mid X$, if we denote with $\gamma_0(p, X)$ the $p$-th conditional quantile and with $\partial_p\gamma_0(p, X)$ its derivative with respect to $p$, then we have $v_\rho(X) = - (\partial_p\gamma_0(p, X))^{-1}$.} The Neyman-orthogonality condition above includes $v_\rho(X)$, which was previously equal to $-1$. Here $v_\rho(X)$ is needed to account for the effect of $\gamma$ on the residual $\rho(W,\gamma)$.

A key innovation of our work is to note that Equation (ref) can be viewed as the first order condition to the following optimization problem:

where the second equality follows by linearity of expectations, and the third equality follows by equation (ref). Thus, $\alpha_{0}$ minimizes the expectation of an objective function that depends on $\alpha$ only through the functional of interest $m(W,\alpha)$ and $\alpha(X)$. As with equation (ref), minimizing this objective function does not require any knowledge of the form of $\alpha_0$.

When $v_{\rho}(X)\neq-1$ the $\alpha_0$ will not be the Riesz representer $v_m(X)$. Instead, $\alpha_0$ can be interpreted as minimizing weighted least squares criterion that depends on the Riesz representer. As shown in ichimura2022influence,

Thus $\alpha_0(X)$ minimizes a weighted least square criterion with weight $-v_{\rho}(X)$ and the variable being predicted given by $-v_m(X)/v_{\rho}(X)$. For this reason we refer to $\alpha_0(X)$ as a weighted Riesz regression.

Though the objective functions of equations (ref) and (ref) differ only by a constant, only equation (ref) possesses the desirable properties that we set out to accomplish of depending solely on known functions of $\alpha$. The objective in equation (ref) was not given in ichimura2022influence.

\vsedit{In some cases there will be a function $\bar{v}_{\rho}(W)$ such that $\operatorname{E}[\bar{v}_{\rho }(W)\mid X]=v_{\rho}(X).$ By iterated expectations, the objective function is not affected by replacing $v_{\rho}(X)$ with $\bar{v}_{\rho}(W),$ because \[ \operatorname{E}[-2m(W,\alpha)-v_{\rho}(X)\alpha(X)^{2}]=\operatorname{E}[-2m(W,\alpha)-\bar{v}_{\rho} (W)\alpha(X)^{2}]. \] In practice, it may be easier to minimize the objective function that depends on $\bar{v}_{\rho}(W)$ to avoid having to estimate $v_{\rho}(X)=\operatorname{E}[\bar{v}_{\rho}(W)\mid X]$. For this reason, we focus on a sample objective function that depends on an estimator $\hat{v}_\rho(W)$ of $\bar{v}_\rho(W)$ that is allowed to take $W$ as input, instead of just $X$.}

To obtain an estimate of $\alpha_{0}$, we replace the sample criterion in step (2b) of the algorithm in subsection (ref) with:

for $\mathcal{A}_{n}\subset\Gamma$, where $\hat{v}_\rho(W)$ is an estimator \vsedit{of $v_\rho(X)$, in the sense that:

} We refer to this $\hat\alpha_\ell$ as a weighted Riesz regression estimator. When $v_\rho(x)$ is known, we can use $\hat{v}_\rho(W)=v_\rho(X)$, and the expectation of this objective function is (ref), plus a penalty. Steps (3) and (4) are modified accordingly to:

Below we give regularity conditions and a theorem to extend the results of subsection (ref) to the case where $v_\rho(x)$ is unknown and needs to be estimated. \vsedit{For simplicity of exposition we will only consider the case where the estimator is defined over a growing sieve space $\mathcal{A}_n$ and no regularization is used, i.e. $\Lambda_r(\alpha)=0$:

Our estimation rate can easily be extended to regularized estimation with appropriate regularization weight.} Let $\mathcal{V}_n$ denote the function space in which the estimator $\hat{v}_\rho$ is restricted to lie in. Let $\alpha_{\ast}$ be any function in $\mathcal{A}_n$ (e.g. we will typically consider $\alpha_{\ast}=\inf_{\alpha\in \mathcal{A}_n} \|\alpha_{\ast}-\alpha_0\|$, but $\alpha_{\ast}$ can in fact be any function that is not chosen based on the samples). Define: \vsedit{

}

We remark that the uniform upper bound of $1$ can be replaced by any constant upper bound $b$ and the rate that we achieve will be identical, up to an extra multiplicative factor $b$, via a standard re-scaling argument (i.e. applying our result to re-scaled version of the original problem and then scaling back the guarantee).

Moreover, we note that if a separate sample was used to estimate $\hat{v}_{\rho}$ and not the same as the one that was used for $\hat{\alpha}$, then we can weaken the theorem to only require $\delta_n$ to upper bound the critical radius of $\mathrm{star}(\mathcal{A}_n-\alpha_{*})$ and not $\mathrm{star}(\sqrt{\mathcal{V}_n}\cdot (\mathcal{A}_n-\alpha_{*}))$. Note that a sufficient condition for Equation (ref) is that $|v_\rho(X)|\geq \lambda$, almost surely. However, for most function spaces $\mathcal{A}_n$, this condition can be satisfied by more benign assumptions. For this it is crucial that we only invoked the property at the difference of two functions that both lie in $\mathcal{A}_n$ and not for instance for $\alpha - \alpha_0$ (since $\alpha_0$ can potentially lie outside of the space). For instance, if the functions in $\mathcal{A}_n$ are are high-dimensional linear functions $\phi(X)'\beta$, then Equation (ref) is satisfied, with $\lambda = \mu/C$, if:

since then, if we let $\alpha=\phi(\cdot)'\beta$ and $\alpha_{\ast}=\phi(\cdot)'\beta_{\ast}$ and $\nu=\beta - \beta_{\ast}$, then:

The following assumption provides regularity conditions on the residual $\rho(w, \gamma)$ and the functional of interest $m(w, \gamma)$ for the generalized regression case.

The next condition allows for $\rho(W,\gamma)$ to be nonlinear in $\gamma.$

This assumption imposes the usual faster than $n^{-1/4}$ convergence rate for $\hat{\gamma}_{\ell}$ when $\rho(w,\gamma)$ is nonlinear in $\gamma$ but does not require that rate when $\rho(W,\gamma)$ is linear in $\gamma$.

We have the following large sample inference result under these conditions.

Nonlinear Effects of Multiple Regressions

Some important objects of interest are expectations of nonlinear functionals of multiple regressions. In this Section we give Auto-DML for such effects. Such effects have the form $\theta_{0}=\operatorname{E}[m(W,\gamma_{0})]$ where $m(w,\gamma)$ is nonlinear in a possible value $\gamma$ of multiple generalized regressions $(\gamma_{1}(X_{1}),...,\gamma_{J}(X_{J}))^{\prime}$ with regressors $X_{j},$ residual $\rho_{j}(W,\gamma_{j})$, and $\Gamma_{j}$ specific to each regression $\gamma_{j}(X_{j})$. The corresponding orthogonal score like is like that of subsection (ref) except that the bias correction is a sum of $J$ terms with the $j^{th}$ term being the bias correction for the learner of $\gamma_{j}$. Similarly to newey1994asymptotic, pg. 1357, the orthogonal score is

\vscomment{We should explain here using the same "first-order" correction intuition, why we are using the directional derivative as the moment. Currently we don't offer any intuition}

Each of the terms in the bias correction can be estimated by the product of a learner $\hat{\alpha}_{j\ell}(X_{j})$ and the residual $\rho_{j}(W,\hat{\gamma}_{j\ell})$, but now the learner $\hat{\alpha}_{j\ell}(X_{j})$ differs from the one given in section (ref) in the way needed to account correctly for nonlinearity of $m(W,\gamma)$ in $\gamma$. The difference is that in the objective function for $\hat{\alpha}_{j\ell}(X_{j})$ the functional of interest $m(w,\alpha)$ is replaced by an estimated Gateaux derivative with respect to the $j^{th}$ component of $\gamma.$ Let \[ \hat{D}_{j}(W,\alpha_{j})=\left. \frac{d}{d\tau}m(W,\hat{\gamma }_{\ell}+\tau e_{j}\alpha_{j})\right\vert _{\tau=0} \] be such a Gateaux derivative estimator, where $e_{j}$ denotes the $j$-th column of the identity matrix. This derivative will often be straightforward to calculate as an analytic derivative with respect to the scalar $\tau.$ When $m(w,\gamma)$ is linear in a single $\gamma$ this derivative just evaluates $m(W_{i}, \gamma)$ at $\gamma=\alpha$ giving the $m(W,\alpha)$ of subsection (ref).

To obtain $\hat{\alpha}_{j\ell}(X_{j})$ we also make use of an estimated derivative $\hat{v}_{\rho j}(W_{i})$ of $\rho_{j}(W,\gamma_{j})$ with respect to $\gamma_{j}$ at $\hat{\gamma}_{j\ell}.$ Then $\hat{\alpha}_{j\ell}$ is given by

where $\mathcal{A}^{j}_n$ is the set of approximating functions for $\alpha_{j}.$ As with linear $m(w, \gamma)$ this $\hat{\alpha}_{j\ell}$ depends just on $m(w,\gamma)$ and the first step. Thus $\hat{\alpha}_{j\ell}$ is automatic, in the same way as in section (ref), in only requiring $m(w,\gamma)$ and the regression residual $\rho_{j}(W_{i},\gamma_{j})$ for its construction.

Below we give two examples of this setting:

It is straightforward to obtain a convergence rate for $\hat{\alpha}_{j\ell}$ analogous to Theorem (ref). The following result does so while accounting for the presence of $\hat{\gamma}$ in $\hat{D}_{j}(W_{i},\alpha_{j}).$ For notational convenience we suppress the $j$ subscripts.

The construction of $\hat{\theta}$ is analogous to that in subsection (ref) with the bias correction term being the sum of terms for each $\gamma_{j}$ in $\gamma.$ That is,

It is straightforward to specify conditions for asymptotic normality of $\hat{\theta}$ and consistency of $\hat{V}$ by combining the conditions of section (ref) with the convergence rate result of Corollary (ref). For relative simplicity we give a result only for neural net learners. We also assume for simplicity that each $X_{j}$ has the same dimension $d.$

This condition is analogous to Assumption (ref).

This condition is analogous to Assumption (ref).

Empirical Application

To illustrate our methods, we study whether applicant race is a significant predictor of banks' mortgage denial decisions. Following munnell1996mortgage, we use the publicly available Boston Home Mortgage Disclosure Act (HDMA) dataset. The dataset contains information on 2,925 mortgage applications made in 1990 in the Greater Boston metropolitan area. We restrict attention to black and white applicants in single-family households (excluding other racial minorities and multi-family residences), which reduces our sample size to 2,380 observations.

Our outcome of interest is an indicator $Y = 1$ if the mortgage application was denied. Our regressor of interest is an indicator $D = 1$ if the applicant is black. We also have access to a vector of covariates, which we denote by $Z$, containing financial and other characteristics of the applicant that banks may factor into their mortgage denial decisions. These include monthly debt to income (DTI) ratio; monthly housing expenses to income (HTI) ratio; loan to assessed property value (LTV) ratio; a categorical variable for “bad” consumer credit score with 6 categories (1 if no slow payments or delinquencies, 2 if one or two slow payments or delinquencies, 3 if more than two slow payments or delinquencies, 4 if insufficient credit history for determination, 5 if delinquent credit history with payments 60 days overdue, and 6 if delinquent credit history with payments 90 days overdue); a categorical variable for “bad” mortgage credit score with 4 categories (1 if no late mortgage payments, 2 if no mortgage payment history, 3 if one or two late mortgage payments, and 4 if more than two late mortgage payments); an indicator for public record of credit problems including bankruptcy, charge-offs, and collective actions; an indicator for denial of application for mortgage insurance; three indicators for self-employed, single, and high school graduate, the 1989 Massachusetts unemployment rate in the applicant's industry, and an indicator for whether the unit is a condominium.

Table (ref) reports the sample means and standard deviations of the variables used in the analysis. The probability of being denied a mortgage is 19 percentage points higher for black applicants than for white applicants. However, black applicants are also more likely to have financial and socio-economic characteristics linked to mortgage denial, as Table (ref) shows. For example, black applicants have higher (worse) consumer and mortgage credit indices on average, and are more likely to have a public record of credit problems and to be single. We would like to test whether the racial differences in probability of mortgage denial persist once we control for these covariates.

To showcase the versatility of our method, we present results for three estimands:

We estimate these parameters using AutoDML, where both $\hat{\gamma}$ and $\hat{\alpha}$ are neural net learners. For the difference in probability, we have $v_\rho(X) = -1$, since $\gamma_0$ is a conditional mean. For the average difference in log-odds and the average difference in odds, $v_\rho(X) = -\lambda(\gamma_0(X))$, where $\lambda(t) := e^{-t}/(1 - e^{-t})^2$ is the standard logistic PDF, which we estimate by replacing $\gamma_0$ with a preliminary estimate $\hat\gamma_{\text{prel}}$, also based on a neural net learner. We describe the architecture and training hyperparameter choice in detail in Appendix (ref).

Table (ref) presents the results of our main analysis. Once we control for covariates, the difference in probability of mortgage denial decreases from 19 to 8 percentage points. If we look at the average log-odds or odds instead, we observe differences of 0.829 or 0.157, respectively. These differences are all estimated to be statistically different from 0 at the 1% significance level.

A slight modification of our method allows us to estimate average differences for subgroups of applicants with certain characteristics (analogous to conditional average treatment effects or CATEs). Suppose we want to estimate an average effect for applicants with $Z_j = z$ for a particular covariate $Z_j$. To obtain these, we weight the Neyman orthogonal estimating equation (ref) as follows: \[\hat{\theta}(z) = \frac{1}{\sum_{i}^n \omega_i(z)}\sum_{\ell=1}^{L}\sum_{i\in I_{\ell}} \omega_i(z) \{m(W_ {i},\hat{\gamma}_{\ell})+\hat{\alpha}_{\ell}(X_{i} )\rho(W_{i},\hat{\gamma}_{\ell})\}.\] When $Z_j$ is a categorical variable, we take $\omega_i(z) = 1\{Z_j = z\}$. When $Z_j$ is continuous, we take $\omega_i(z) = K((Z_i - z)/h)$ for a kernel function $K$ and a small but fixed bandwith $h$.\footnote{chernozhukov2019global analyze a localized version of this parameter, that is, the limit as $h \to 0$, which is beyond the scope of this paper.} Figure (ref) presents the racial differences in probability, average log-odds and average odds by values of the consumer credit index and the monthly DTI ratio. Remarkably, we estimate the racial differences in all three estimands to be higher for applicants with a delinquent credit history or with insufficient credit history (although the latter is quite imprecisely estimated). The racial differences appear to be constant for most of the range of the monthly DTI ratio variable, except values below 0.075 for which it is also imprecisely estimated.

Monte Carlo Simulations

Based on the Empirical Application

First, we analyze the performance of our method in the setting of our empirical application. We redraw the covariates $Z$ based on a generative adversarial network (GAN) trained on the real mortgage data. We use an elastic-net Logit, with penalties chosen by cross-validation, to estimate the outcome regression $\Pr(Y = 1 \mid D, Z)$ and the propensity score $\Pr(D = 1 \mid Z)$ in the real mortgage data, which we take as ground truth in our simulations. We present results for the three estimands on interest in Section (ref): the difference in probability, the average difference in log-odds and the average difference in odds.

Table (ref) presents simulation results over 1,000 draws for $n = 2,000$ and $n = 10,000$. The first column shows the non-parametric $R^2$ for $\gamma$, defined as $R^2(\gamma) = 1 - \operatorname{E}[(\hat{\gamma}(X) - \gamma_0(X))^2]/\mathrm{Var}(\gamma_0(X))$, where the expectation is evaluated over a test set not used to estimate $\hat{\gamma}$. The second column shows the same non-parametric $R^2$ metric for $\alpha$. We also give the mean absolute error (MAE), bias, standard deviation (sd), the average standard error to standard deviation ratio (se/sd) and the coverage of a 95% confidence interval (covg.).

For the estimands we study, the Riesz regression $\alpha_0$ can be characterized explicitly based on the Riesz representer for the Average Treatment Effect (ATE), \[v_m(X) = \frac{D}{\Pr(D = 1 \mid Z)} - \frac{1-D}{1 - \Pr(D = 1 \mid Z)}.\] For the difference in probability estimand, which is an average linear effect for a conditional mean (Section (ref)), $\alpha_0(X) = v_m(X)$. For the average difference in log-odds, which is an average linear effect for a generalized regression (Section (ref)), the weighted Riesz regression is $\alpha_0(X) = v_m(X)/(-v_\rho(X))$, where $v_\rho(X) = -\lambda(\gamma_0(X))$. Finally, for the average difference in odds, which is an average non-linear effect for a generalized regression (Section (ref)), we have $\alpha_0(X) = e^{\gamma_0(X)}v_m(X)/(-v_\rho(X))$.

Our automatic debasing method does not make use of this explicit characterization of $\alpha_0$. To benchmark our results, we compare its performance to an estimator that uses the explicit characterization of $\alpha_0$. In the ATE setting, this is known as an Augmented Inverse Propensity Weighting (AIPW) estimator, so we will refer to these as AIPW-like. To build the AIPW-like estimator of $\theta_0$ we plug learners of the outcome and treatment propensities $\Pr(Y = 1 \mid D, Z)$ and $\Pr(D = 1 \mid Z)$ into the formula for $\alpha_0$; we try both a non-parametric version based on neural nets (the same architecture and hyperparameters as our main specification) and a “well-specified” version, where we use the same elastic-net Logit that we used to build the ground truth.

A comparison between auto-DML and the AIPW-like benchmark sheds light on the advantages of our automatic approach. For large sample sizes, $n = 10,000$, both methods perform comparably well. There is no loss in efficiency between our main specification (which uses a non-parametric, NN-based method) and the correctly-specified AIPW estimator (which uses correctly specified parametric learners for $\gamma$ and $\alpha$). Our automatic debiasing method achieves close to nominal coverage, whereas the AIPW-like method gets worse coverage when we use the non-parametric, neural net specification. A reason for that could be that, in the explicit characterization of $\alpha_0$, we are plugging numbers that are close to zero into a denominator (such as the propensity score $\Pr(D = 1 \mid Z)$ or the logit pdf $v_\rho(X)$), so that estimation error amplifies. This is reflected into the large negative non-parametric $R^2$ for $\alpha$ when $n = 2,000$.

Consistent with these results, in work that followed up on the first version of this paper (chernozhukov2022riesznet), we found that our automatic debiasing method using neural net and random forest Riesz regressions performed much better than state of the art methods based on inverse propensity score weighting in Monte Carlo experiments. singh2023double also found that automatic debiased estimators performed better than plugin-based estimators in the setting of local average treatment effects.

Additional Simulations

We present an additional set of simulations when $\gamma_0$ is a quantile of the conditional distribution of $Y \mid X$ for a continuous outcome $Y$. For simplicity, we will focus on the conditional median. As discussed in Section (ref), this corresponds to the generalized residual $\rho(W, \gamma) = 0.5 - 1(Y < \gamma(X)),$ which is a sub-derivative of the “check” loss function $\ell(W, \gamma) = 0.5|Y - \gamma(X)|$.

We will consider four simulation settings. In the first setting, the object of interest will be the average difference in conditional median for a binary treatment $D$, that is, $\theta_0 = \operatorname{E}[\gamma_0(1, Z) - \gamma_0(0, Z)]$. We assume that $\Pr(D = 1 \mid Z) = \mathrm{logit}(-0.1 + 0.5 Z_1 - 0.2 Z_2)$, and that $Y = \mu_Y(X) + \epsilon$ with $\mu_Y(X) = 0.5D - 0.2 D \times Z_3 + 0.3 Z_2$ and $\epsilon \sim N(0, 1)$. The second setting considers the same data generating process, but we use a twice-differentiable quantile loss based on Epanechnikov kernel smoothing, due to he2023smoothed, rather than the non-differentiable check loss function. In the third setting, our object of interest is the average derivative of the conditional median with respect to a continuous treatment $D$, that is, $\theta_0 = \operatorname{E}[\partial_d \gamma_0(D, Z)]$, where we draw $D = \mu_D(Z) + \eta$ for $\mu_D(Z) = -0.1 + 0.5 Z_1 - 0.2 Z_2$ and $\eta \sim N(0, 1)$; the distribution of $Y \mid X$ is as before. The last setting uses the same data generating process, but it focuses on a non-linear parameter, the average derivative squared $\theta_0 = \operatorname{E}[(\partial_d \gamma_0(D, Z))^2]$. This parameter is useful, for example, in testing whether $\partial_d \gamma_0(D, Z) = 0$ with probability 1. In all settings, to work with a realistic data generating process for the covariates, we will draw $Z$ from the same GAN trained on the real mortgage data that we used in the previous subsection.

The results over 1,000 simulation draws are presented in Table (ref). In all the settings we consider, our estimator performs well, with low bias and coverage confidence intervals close to the nominal 95% level.