Assumption-lean falsification tests of rate double-robustness of double-machine-learning estimators

Assumption-lean Falsification Tests of Rate Double-Robustness of Double-Machine-Learning Estimators

\newrefsection

{ Key words: Econometrics, Causal Inference, Machine Learning, Doubly-Robust Functionals, Higher-Order $U$-Statistics, Higher-Order Influence Functions}

\allowdisplaybreaks

Introduction

Suppose a data analyst has constructed and published a nominal $1 - \alpha$ large sample Wald confidence interval (CI) for a mean-square continuous linear functional $\psi$ of a conditional expectation $b (x) = \mathsf{E} [Y | X = x]$ with the Wald CI centered by a doubly-robust Double Machine Learning (DML) estimator $\widehat{\psi}_{1}$ of $\psi$ chernozhukov2018double (see Section (ref) for a formal definition of $\widehat{\psi}_{1}$). The estimator $\widehat{\psi}_{1}$ will depend on estimates $\widehat{b}(x)$ and $\widehat{p}(x)$ of two functions of $x$: $b(x)$ itself and the function $p(x)$ occurring in the Riesz representer of the linear functional. Were it achievable, our goal would be to construct an assumption-lean (i.e. essentially assumption-free) empirical test, with non-trivial power against certain alternatives, of the null hypothesis that the true asymptotic coverage for $\psi$ of the above nominal $1 - \alpha$ Wald CI is greater than or equal to $1 - \alpha$ under repeated sampling. By definition, assumption-lean tests make no complexity-reducing assumptions (such as smoothness or sparsity) on $b (x)$ or $p (x)$. If such a test rejects (with a very small p-value), we would have falsified (or more precisely, have strong evidence) that the true coverage of the published Wald CI is less than nominal. Henceforth, we will say a large sample Wald CI is valid if and only if the above null hypothesis is true. Unfortunately, following robins1997toward, such tests do not exist; for intuition, see Section (ref).

We therefore adopt the following less ambitious, but partially achievable goal: Our (new) goal is to construct an empirical assumption-lean test that can falsify an analyst's justification for the claim that their nominal Wald CI centered at a DML estimator is valid. If falsified, the analyst should then retract any claim of validity. We refer to our goal as partially achievable because our test can only falsify certain types of justification. To formally characterize which types, we first must review the properties of DML estimators.

A necessary and (essentially) sufficient condition for the validity of a Wald CI for $\psi$ centered at a DML estimator is that the (asymptotic) bias of the estimator is $o_{p}(n^{-1/2})$. chernozhukov2018double showed that a sufficient (but not necessary) condition for validity is that the (weighted) $L_{2}(\mathsf{P})$ rate of convergence $n^{-\kappa_{b}}$ of $\widehat{b}$ to $b$ multiplied by the rate of convergence $n^{-\kappa_{p}}$ of $\widehat{p}$ to $p$ is $o_{p}(n^{-1/2})$ or, equivalently, $\kappa_{b} + \kappa_{p} > 1 / 2$. Conditions in similar spirit also appeared in earlier works in the econometrics and statistics literature, such as chen2008semiparametric, chen2015sieve. This property was termed “rate double-robustness” in smucler2019unifying.

The main technical result of our paper is that we construct an assumption-lean empirical test, with power against certain alternatives, of the null hypothesis that $\kappa_{b} + \kappa_{p} > 1 / 2$, or in words, “that rate double-robustness is true”. Thus, our test can potentially falsify the justification of any analyst who uses, either explicitly or implicitly, “rate double-robustness” to justify the validity of her Wald CI. An example would be an analyst who (i) makes explicit, restrictive assumptions on both the complexities of the functions $b$ and $p$ (e.g. in terms of smoothness or sparsity) and on the algorithms used in their fitting followed by (ii) an appeal to theorems that guarantee $\kappa_{b} + \kappa_{p} > 1 / 2$ under these restrictions. For instance, when $b$ and $p$ are fit by minimizing a (weighted) penalized empirical (squared $L_{2}$-) loss (see (ref)) using deep neural networks farrell2021deep, chen2020causal, xu2022deepmed, $L_{2}(\mathsf{P})$-convergence rates satisfying $\kappa_{b} + \kappa_{p} > 1 / 2$ can be proved by assuming $b$ and $p$ live in sufficiently smooth H\"{o}lder spaces schmidt2020nonparametric\footnote{If our test rejects the hypothesis that $\kappa_{b} + \kappa_{p} > 1 / 2$, it also rejects the hypothesis that any assumed collection of restrictions on complexity of and fitting algorithms for $b$ and $p$ that imply $\kappa_{b} + \kappa_{p} > 1 / 2$ are all true.}. A second example, especially common in the applied literature, is an analyst who reports a nominal Wald CI centered on a DML estimator without any explicit discussion of its validity, other than citing chernozhukov2018double. We also regard such an analyst's implicit justification for her Wald CI as being by appeal to rate double-robustness.

On the other hand, there are “justifications” for Wald CI validity that cannot be falsified by our assumption-lean tests rejecting $\kappa_{b} + \kappa_{p} > 1 / 2$ with a very small p-value. Specifically, some recent papers have proposed novel DML estimators that, under very restrictive assumptions on both $b$ and $p$ and on the algorithms used in their estimation, have bias $o_{p} (n^{-1/2})$ and thus can center valid Wald CIs, even though $\kappa_{b} + \kappa_{p} < 1 / 2$ newey2018cross, kennedy2020towards. Therefore an analyst who justifies the validity of her Wald CIs by appeal to these restrictive assumptions is not at all surprised to learn that the null hypothesis $\kappa_{b} + \kappa_{p} > 1 / 2$ is false. These novel estimators are discussed both later in the Introduction and in Section (ref).

In the remainder of the paper, we restrict the functional $\psi$ to the class of Mixed-Bias or Doubly-Robust (DR) functionals rotnitzky2021characterization. This class strictly includes both (i) the class of mean-square continuous functionals that can be written as an expectation of an affine functional of a conditional expectation studied by chernozhukov2022automatic and (ii) the class of functionals studied by robins2008higher\footnote{The class of Mixed-Bias or DR functionals considered in this paper is itself contained in the class of functionals discussed in Section 5 of chernozhukov2022locally, that constitutes the most general class of functionals for which there exist first-order DR estimators.}. Many DR functionals are of substantive scientific and economic interest, see the examples after Definition (ref). Our unified treatment of the entire DR functional class requires that we use rather abstract notation. In order to prevent such abstract notation from making the reader miss the forest for the trees, we shall complete the introduction by using a familiar functional to motivate both our goals and our methodology. Furthermore, detailed regularity conditions will be suppressed in the Introduction to facilitate the exposition.

The motivation for our paper is best described by the following inferential quandary faced many (perhaps dozens) of times daily by data analysts employed by large tech companies. The quandary arises when the analyst needs to estimate a population average effect $\mathsf{E} [Y (a = 1)] - \mathsf{E} [Y (a = 0)]$ of a dichotomous treatment $A$ on a response $Y$ from observational data. Here $Y (a)$ is the counterfactual outcome under treatment level $a$. Often, data on a very high-dimensional vector $X$ (with dimension $d$) of pretreatment covariates are available and deemed sufficient for ignorability $Y (a) \mathop{\perp\!\!\!\!\perp} A | X$ to hold, thereby identifying\footnote{The minus sign plays no essential role and is added for notational convenience that will be made clear in Definition (ref).} $\psi^{a} \equiv -\mathsf{E} [Y (a)]$ as $- \mathsf{E} [b_{a} (X)]$ with $b_{a} (x) = \mathsf{E} [Y | X = x, A = a]$ when positivity holds. chernozhukov2018double argued persuasively that in an effort to obtain valid inference (i.e. confidence intervals) for $\psi^{a}$ [and thus for the average treatment effect (ATE) $- \psi^{a = 1} + \psi^{a = 0}$] with very high dimensional $X$, one should use (cross-fitted) DML estimators $\widehat{\psi}_{\mathsf{cf}, 1}^{a}$ to center nominal $1 - \alpha$ large sample Wald CI $\widehat{\psi}_{\mathsf{cf}, 1}^{a} \pm z_{\alpha / 2} \widehat{\mathsf{s.e.}} (\widehat{\psi}_{\mathsf{cf}, 1}^{a})$ where $z_{\alpha / 2}$ is the $\alpha / 2$ upper-quantile of a standard normal random variable and $\widehat{\mathsf{s.e.}}(\widehat{\psi}_{\mathsf{cf}, 1}^{a})$ is the estimator of the standard error $\mathsf{s.e.}(\widehat{\psi}_{\mathsf{cf}, 1}^{a})$ of $\widehat{\psi}_{\mathsf{cf}, 1}^{a}$ given in Proposition (ref) below. As a result, DML estimators rapidly became the standard in high tech.

For DR functionals, DML estimators combine the benefits of cross-fitting (cf), double robustness, and machine learning of nuisance parameters chernozhukov2018double, chernozhukov2022locally. Henceforth, for notational convenience, we remove the $a$ index by restricting to the case $a = 1$ so, for example, $\psi^{a=1}$ becomes $\psi$ and $b_{a}$ becomes $b$. Then, to compute a DML estimator, the data is randomly divided into two (or more) samples -- the estimation sample of size $n$ and the training or nuisance/training sample of size $n_{\mathsf{tr}} = N - n$ with $1 - c > n / N > c$ for some $c \in (0, 1)$. To simplify the exposition we take $c = 0.5$. Estimators $\widehat{b} (\cdot)$ and $\widehat{p} (\cdot)$ of $b (\cdot)$ and inverse propensity score $p (\cdot) \coloneqq 1 / \mathsf{E}[A|X = \cdot]$, where $p^{-1}$ is assumed to be strictly bounded between $(0, 1)$, are computed from the training sample data using modern black-box highly nonlinear machine learning algorithms, often deep neural networks. In semiparametric statistics literature, $b$ and $p$ are referred to as nuisance parameters/functions. Let $\mathsf{P}_{n}$ denote a sample average over the estimation sample. Then the doubly-robust one step estimator $\widehat{\psi}_{1} = \widehat{\psi} + \mathsf{P}_{n} [\widehat{\mathsf{IF}}_{1, \psi}] = \mathsf{P}_{n} [- \widehat{b} (X) - A \widehat{p} (X) (Y - \widehat{b} (X))]$ is constructed by adding to an initial plug-in estimator $\widehat{\psi} = \mathsf{P}_{n} [-\widehat{b} (X)]$ the sample average of an estimate of the first order influence function $\mathsf{IF}_{1, \psi} = - b (X) - A p (X) (Y - b (X)) - \psi$ of $\psi$\footnote{In the econometrics literature, $\widehat{\psi}$ is referred to as a first stage estimator, $\widehat{\psi}_{1}$ as the second stage estimator, and $\mathsf{P}_{n} [- A \widehat{p} (X) (Y - \widehat{b}(X))]$ is the debiasing term that makes $\widehat{\psi}_{1}$ a doubly-robust estimator satisfying Neyman orthogonality chernozhukov2018double, chernozhukov2022locally.}. The cross-fit DML estimator $\widehat{\psi}_{\mathsf{cf}, 1}$ is the arithmetic average of $\widehat{\psi}_{1}$ and $\bar{\widehat{\psi}}_{1}$, where $\bar{\widehat{\psi}}_{1}$ is computed like $\widehat{\psi}_{1}$ but with the roles of training and estimation samples switched. We define a standard DML estimator to be a DML estimator where $\widehat{b}$ and $\widehat{p}$ are separately estimated, each using data from the entire training sample; see Section (ref). This allows us to distinguish standard DML estimators from, for instance, “DCDR estimators” of newey2018cross that estimate $b$ and $p$ from separate non-overlapping subsamples of the training sample.

The analyst's inferential quandary is how to justify the claim that $\widehat{\psi}_{\mathsf{cf}, 1}\pm z_{\alpha /2}\widehat{\mathsf{s.e.}}(\widehat{\psi}_{\mathsf{cf}, 1})$ is a valid large sample $1-\alpha $ Wald CI. For it to be valid the following are necessary: (i) $\widehat{\psi}_{\mathsf{cf}, 1}$ is asymptotically normal, (ii) $\widehat{\mathsf{s.e.}} (\widehat{\psi}_{\mathsf{cf}, 1}) / \mathsf{s.e.}(\widehat{\psi}_{\mathsf{cf}, 1})$ converges to $1$ in probability and (iii) the (asymptotic) bias of $\widehat{\psi}_{\mathsf{cf}, 1}$ is of smaller order than $\mathsf{s.e.} (\widehat{\psi}_{\mathsf{cf}, 1})$. Since $\widehat{\psi}_{\mathsf{cf}, 1}$ is the average of $\widehat{\psi}_{1}$ and its “twin” $\bar{\widehat{\psi}}_{1}$, it must be the case that $\widehat{\psi}_{1} \pm z_{\alpha / 2} \widehat{\mathsf{s.e.}} (\widehat{\psi}_{1})$ is also a valid large sample $1 - \alpha$ Wald CI for $\psi$, therefore satisfying (i) - (iii) with $\widehat{\psi}_{1}$ substituted for $\widehat{\psi}_{\mathsf{cf}, 1}$. Since $\mathsf{s.e.} (\widehat{\psi}_{1})$ is order $n^{-1/2}$, it is necessary that its bias is $o_{p} (n^{-1/2})$ to satisfy (iii) above. As discussed in the literature newey2018cross, by far the most difficult of the three assumptions to satisfy and thus to justify is (iii). As mentioned above, for a standard DML estimator, if $\kappa_{b} + \kappa_{p} > 1 / 2$ holds, then (iii) holds. Specifically, that $\kappa_{b} + \kappa_{p} > 1 / 2$ implies the bias of $\widehat{\psi}_{1}$ is $o_{p} (n^{-1/2})$ is a consequence of applying Cauchy-Schwarz (CS) inequality to upper bound the bias, as we show next. The exact conditional bias of $\widehat{\psi}_{1}$ given the training sample, denoted as $\mathsf{Bias} (\widehat{\psi}_{1})$, is $\mathsf{E} [A (\widehat{b} (X) - b (X)) (\widehat{p} (X) - p (X))] = \int p^{-1}(x) (\widehat{b} (x) - b (x)) (\widehat{p} (x) - p (x)) \mathrm{d} F (x)$, where $F$ denotes the distribution of $X$. Here and henceforth, unless stated otherwise, all expectations will be understood to be conditional on the training sample, although that fact is suppressed in the notation for brevity. It then follows from the CS inequality that

We refer to the RHS of the above display as the (conditional on the training sample) Cauchy-Schwarz (CS) bias, denoted as $\mathsf{CSBias}(\widehat{\psi}_{1})$, of $\widehat{\psi}_{1}$. More generally, for any pair of positive functions $w = (w_{b}, w_{p})$ of $x$ strictly bounded from above and below, define

Then applying H\"{o}lder inequality, as in (ref) later in our paper, we have that $\mathsf{CSBias}^{w}(\widehat{\psi}_{1})$ and $\mathsf{CSBias} (\widehat{\psi}_{1})$ are equal up to a multiplicative positive constant, except that when $w_{b} = w_{p} = p^{-1}$, the equality is exact. Hence $\mathsf{CSBias}^{w} (\widehat{\psi}_{1}) = o_{p} (n^{-1/2})$, if and only if $\mathsf{CSBias}(\widehat{\psi}_{1})=o_{p}(n^{-1/2})$, if and only if rate double-robustness holds ($\kappa_{b} + \kappa_{p} > 1 / 2$, where we refer to $n^{-\kappa_{b}}$ and $n^{-\kappa_{p}}$ as the $p^{-1}$-weighted $L_{2} (\mathsf{P})$ convergence rates of $\widehat{b}$ and $\widehat{p}$ to $b$ and $p$). Thus $\mathsf{CSBias}^{w}(\widehat{\psi}_{1})=o_{p}(n^{-1/2})$ also implies $\mathsf{Bias}(\widehat{\psi}_{1})=o_{p}(n^{-1/2})$.

Main technical contributions

If we can empirically reject the null hypothesis $\mathsf{NH}_{0, \mathsf{CS}}: \mathsf{CSBias} (\widehat{\psi}_{1}) = o_{p}(n^{-1/2})$ encoding rate double-robustness, then $\kappa_{b} + \kappa_{p} \leqslant 1 / 2$, and we falsify the justification of the validity of a Wald CI centered at the standard DML estimator $\widehat{\psi}_{1}$. The main technical contribution of this paper is to construct a test of $\mathsf{NH}_{0, \mathsf{CS}}$.

The reader might rightfully complain at this point that finite sample tests of asymptotic hypotheses such as $\mathsf{NH}_{0, \mathsf{CS}}$ that concern rates of convergence cannot be constructed. To overcome this problem, we pair each asymptotic null hypothesis of interest with a natural non-asymptotic null hypothesis and then, by convention, formally declare the asymptotic hypothesis (not) rejected if the paired non-asymptotic null hypothesis is (not) rejected. Then we say that the asymptotic null hypothesis has been operationalized by its paired non-asymptotic null hypothesis. In particular, we pair the asymptotic null hypothesis $\mathsf{NH}_{0, \mathsf{CS}}$ with the following non-asymptotic null hypothesis

and construct an $\alpha^{\dag}$-level falsification test of $\mathsf{H}_{0, \mathsf{CS}} (\delta)$, where $\delta > 0$ is chosen by the analyst. With such an operationalized pairing

we will, by convention, declare $\mathsf{NH}_{0, \mathsf{CS}}$ (not) rejected if $\mathsf{H}_{0, \mathsf{CS}} (\delta)$ is (not) rejected. The larger $\delta$ that one chooses, the stronger evidence that rejection of $\mathsf{H}_{0, \mathsf{CS}} (\delta)$ has against the asymptotic null hypothesis $\mathsf{NH}_{0, \mathsf{CS}}$, but the less power one has to reject $\mathsf{H}_{0, \mathsf{CS}} (\delta)$. Suppose the analyst insists that her justification is $\mathsf{CSBias}^{w} (\widehat{\psi}_{1}) = o (n^{- 1 / 2})$ as defined in (ref), then, by the aforementioned equivalence between $\mathsf{CSBias} (\widehat{\psi}_{1})$ and $\mathsf{CSBias}^{w} (\widehat{\psi}_{1})$, we can use the same operationalized pairing (ref) to empirically falsify $\mathsf{CSBias}_{\theta}^{w} (\widehat{\psi}_{1}) = o (n^{- 1 / 2})$. We remark that one could, in principle, choose $\delta$ to be a diminishing sequence as a function of the sample size $n$.

Since we want our test to be assumption-lean, we make essentially no assumptions on the nuisance functions $b$ or $p$, their estimates $\widehat{b}$ or $\widehat{p}$, or the algorithms used to construct $\widehat{b}$ and $\widehat{p}$ from the training sample. To avoid relying on assumptions on $\widehat{b}$ and $\widehat{p}$, the falsification test and its properties are established by conditioning on the training sample, so the training sample is treated as fixed and statements such as $\mathsf{CSBias} (\widehat{\psi}_{1}) = o_{p} (n^{-1/2})$ become $\mathsf{CSBias} (\widehat{\psi}_{1}) = o (n^{-1/2})$. By arguing as in robins1997toward or ritov2014bayesian, in absence of complexity-reducing assumptions on $b$ or $p$, there is no uniformly consistent estimator of either $\mathsf{Bias} (\widehat{\psi}_{1})$ or $\mathsf{CSBias} (\widehat{\psi}_{1})$. However, we will exhibit a functional, denoted as $\mathsf{Bias}_{k} (\widehat{\psi}_{1})$, that is uniformly consistently estimable using a third-order $U$-statistic [derived using the theory of Higher-Order Influence Functions (HOIFs) robins2008higher, liu2017semiparametric] such that one can reject the “$k$-projected” null hypothesis

with nontrivial power. Here $\mathsf{Bias}_{k} (\widehat{\psi}_{1})$ is defined as

where for any $h \in L_{2} (\mathsf{P}_{F})$, $\Pi [h | p^{-1 / 2} \bar{\mathsf{z}}_{k}]$ denotes the population projection of $h$ onto the linear span of a (user-selected) $k$-dimensional dictionary (or basis functions) $p^{-1 / 2} \bar{\mathsf{z}}_{k} = p^{-1 / 2} (\mathsf{z}_{1}, \cdots, \mathsf{z}_{k})^{\top}$. We now show that $|\mathsf{Bias}_{k} (\widehat{\psi}_{1})|$ is a lower bound for $\mathsf{CSBias} (\widehat{\psi}_{1})$ as follows:

where the first inequality again follows from CS inequality and the second inequality follows from the fact that projection contracts norms\footnote{The reason why we do not focus on a different $k$-projected null hypothesis

will be explained in Remark (ref) of Section (ref). In principle, though, one can also consider testing $\mathsf{H}_{0, k, \mathsf{CS}} (\delta)$.}. Hence if $\mathsf{H}_{0, k} (\delta)$ is false, then $\mathsf{H}_{0, \mathsf{CS}} (\delta)$ is false and if $\mathsf{H}_{0, \mathsf{CS}} (\delta)$ is true, then $\mathsf{H}_{0, k} (\delta)$ is true; but the converses of the above two clauses do not necessarily hold. In fact, no test, ours included, can be a consistent test of $\mathsf{H}_{0, \mathsf{CS}} (\delta)$ (that is, no test can have power against all alternatives to $\mathsf{H}_{0, \mathsf{CS}} (\delta)$) unless one makes further possibly incorrect complexity-reducing assumptions on the nuisance functions of $b$ and $p$ and their estimates $\widehat{b}$ and $\widehat{p}$. This again follows from the argument in robins1997toward or ritov2014bayesian. But we also provide an intuitive explanation in Section (ref).

To further illustrate our approach, we can rewrite $\mathsf{Bias}_{k} (\widehat{\psi}_{1})$ as follows (see Appendix (ref)):

where $\Sigma_{k} \equiv \mathsf{E} [A \bar{\mathsf{z}}_{k} (X) \bar{\mathsf{z}}_{k} (X)^{\top}]$. When $\Sigma_{k}$ is known, $\mathsf{Bias}_{k} (\widehat{\psi}_{1})$ can be unbiasedly estimated by the following second-order $U$-statistic:

$\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$ is the second-order influence function of the functional $\mathsf{Bias}_{k} (\widehat{\psi}_{1})$; see Appendix (ref) for a more precise definition. Since $\Sigma_{k}$ is generally unknown, we replace $\Sigma_{k}$ by an estimate $\Sigma_{k}$ from the training sample. However, $\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$ is no longer unbiased for $\mathsf{Bias}_{k} (\widehat{\psi}_{1})$. As a consequence, to protect the level of our test, we replace $\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$ by a third-order $U$-statistic to correct for the additional bias; see Section (ref) for more details.

A preliminary version of the above idea has appeared in liu2020nearly, but the results therein are mainly for (1) the case when the distribution of the potentially high-dimensional covariates $X$ is known, or the so-called {\it semisupervised} setting, and (2) one special case of DR functionals, the expected conditional covariance. In this article, we advance the literature in the following regards.

Literature overview

To the best of our knowledge, the assumption-lean falsification test as constructed in this paper is new in the literature (except for its precursor liu2020nearly), and has different purposes from specification tests newey1985maximum in the econometrics literature. Thus we first mention a subset of the fast-growing literature in statistics and econometrics on estimating and drawing statistical inference for (certain members of) DR functionals using standard DML estimators; due to space limitation, see farrell2015robust, chernozhukov2018double, smucler2019unifying, bradic2019sparsity, chernozhukov2022locally and references therein.

Some recent works also consider further refinement of standard DML estimators newey2018cross, kline2020leave, bradic2019minimax, mcgrath2022undersmoothing, kennedy2020towards. The main distinction of these nonstandard DML estimators from the standard ones is to estimate $b$ and $p$ from separate subsamples of the training sample. Earlier in the introduction, we used the word “novel” rather than “nonstandard” in describing these DML estimators. Under very restrictive, specific complexity-reducing assumptions on $b$ and $p$ and on the algorithms used in their estimation, the estimators may have bias $o (n^{- 1 / 2})$ yet rate double-robustness fails to hold newey1990semiparametric, newey1994large. In the absence of such restrictive, specific complexity-reducing assumptions and fitting algorithms, it is unclear if these nonstandard DML estimator still outperform standard one either in theory or in practice. In fact, a lower bound established recently in balakrishnan2023fundamental shows that standard DML estimators are minimax optimal under an assumption-lean model that imposes no complexity-reducing assumptions on the nuisance functions $b$ or $p$. As a result, this article mainly focuses on the standard DML estimators but these nonstandard ones will also be considered briefly in Section (ref).

Last but not least, we remark that our work is also closely related to the literature on $\sqrt{n}$-consistent estimation and inference for low-dimensional parameters (implicitly) defined via (conditional) moment restrictions involving nonparametric nuisance functions; e.g., see ai2003efficient, ai2007estimation, ai2012semiparametric, chen2015sieve, to name a few. Such parameters encompass the DR functionals studied in this paper, and can be applied to endogeneity settings ai2003efficient, angrist1996identification, tchetgen2020introduction. Extending our framework (specifically the theory of higher-order influence functions) to such more complicated parameters is still an open problem.

Organization of the paper

The remainder of the paper is arranged as follows. In Section (ref), we describe the mathematical setup formally and review the definition and properties of DR functionals recently characterized in rotnitzky2021characterization. In Section (ref), we review the statistical properties of their standard DML estimators, based upon which we motivate and formally define our approach.

The main result of this paper, Section (ref), is to construct a valid $\alpha^{\dag}$-level falsification test of $\mathsf{H}_{0, \mathsf{CS}} (\delta)$ based on a third-order $U$-statistic, denoted as $\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$ (see (ref)), which is the estimated third-order influence function of $\mathsf{Bias}_{k} (\widehat{\psi}_{1})$ but with $\Sigma_{k}$ replaced by $\widehat{\Sigma}_{k}$, following the notation used in robins2008higher; also see Appendix (ref) for derivations.

In Section (ref), we study if the proposed falsification test could be also meaningful for nonstandard DML estimators, whose bias could be $o (n^{-1/2})$ even when $\kappa_{b} + \kappa_{p} \leqslant 1 / 2$ (rate double-robustness is violated). We first argue that the $k$-projected null hypothesis $\mathsf{H}_{0, k} (\delta)$ is a natural hypothesis to falsify. Then we construct a valid $\alpha^{\dag}$-level test of $\mathsf{H}_{0, k} (\delta)$ by modifying the falsification test in Section (ref) using higher order $U$-statistics. Since higher order $U$-statistics are computational costly, we also propose an early-stopping strategy that takes the analyst's computational budget into account. In Section (ref) we present results of simulation studies to evaluate the finite sample performance of our methods. Section (ref) concludes with a discussion of some open problems. Many of the technical details are deferred to the Appendix.

Formal setup and an illustration of our approach

The formal setup is as follows. We observe $N$ i.i.d. copies of the data vector $O = (W, X)$ drawn from some unknown probability distribution $\mathsf{P}_{\theta}$ belonging to a locally nonparametric model

parameterized by the parameter $\theta = (b, p, \theta \setminus \{b, p\})$ with $b, p, \theta \setminus \{b, p\}$ variation independent. To stress the dependence on $\theta$, from here on, we will attach $\theta$ to many symbols that have appeared, such as $\psi (\theta)$, $\mathsf{E}_{\theta}$, $\mathsf{Bias}_{\theta} (\widehat{\psi}_{1})$, $\mathsf{CSBias}_{\theta} (\widehat{\psi}_{1})$, $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$, etc. $X$ is a $d$-dimensional random vector with compact support (with distribution function $F$) whose density $f$ is bounded away from $0$ and $\infty$ on its support, and $d$ is allowed to increase with $N$. Here the maps $b: x \mapsto b(x) \in \mathcal{B}$ and $p: x \mapsto p(x) \in \mathcal{P}$ have range bounded and contained in $\mathbb{R}$. $\mathcal{M}$ is locally nonparametric in the sense that the tangent space for the model at each $\theta \in \Theta$ is equal to $L_{2} (\mathsf{P}_{\theta})$, e.g. when $b, p$ belong to H\"{o}lder balls with certain smoothness.

To avoid extraneous technical issues, we assume that observed data $O$ is bounded with probability 1 (see Remark (ref) for further discussion). We consider functionals (i.e. parameters) $\psi: \theta \mapsto \psi (\theta)$ that possess a (first order) influence function\footnote{The term “influence function” when used without further qualification is to be understood to be the first order influence function.} ichimura2022influence $\mathsf{IF}_{1, \psi} (\theta) = \mathsf{if}_{1, \psi} (O; \theta)$ (and thus a positive and finite semiparametric variance bound newey1990semiparametric) and are contained in the mixed bias or doubly-robust class of functionals of rotnitzky2021characterization defined as follows. Under the locally nonparametric models defined above, the influence function $\mathsf{IF}_{1, \psi} (\theta)$ of $\psi (\theta)$ with respect to the tangent space $L_{2} (\mathsf{P}_{\theta})$ is unique.

To understand the implication of DR functionals, let $\mathsf{P}_{n}$ be the empirical mean operator over the estimation sample and suppose $\theta'$ were an estimate of $\theta$ from the training sample (regarded as fixed, i.e. non-random). It then follows that the one step estimator $\psi (\theta') + \mathsf{P}_{n} \left[ \mathsf{IF}_{1, \psi} (\theta') \right]$ is doubly-robust scharfstein1999rejoinder, robins2001comments, bang2005doubly. That is, by (ref), it is unbiased for $\psi (\theta)$ under $\mathsf{P}_{\theta}$ if either $b = b'$ or $p = p'$. Because of this fact, we will use the term DR functional in this paper, as is done in much of the current literature, instead of the “mixed bias” terminology employed in rotnitzky2021characterization. To ease notation, we will restrict consideration to DR functionals for which $\mathsf{P}_{\theta} (S_{bp} \geqslant 0) = 1$. For a DR functional $\psi^{\dag} (\theta)$ of substantive interest for which $\mathsf{P}_{\theta} (S_{bp} \leqslant 0) = 1$, we will instead analyze $\psi (\theta) = - \psi^{\dag} (\theta)$.

Let $W = (Y, A)$. Below are some examples of DR functionals that are of substantive interest in economics and statistics.

To save space, we refer interested readers to rotnitzky2021characterization for further examples, including Average Treatment Effect on the Treated (ATT) (Example 8 therein), ATE under sensitivity analysis models (Example 3 therein) and etc.

Before proceeding further, we collect some frequently used notation. We use $\mathsf{E}_{\theta} [\cdot], \mathsf{var}_{\theta} [\cdot]$ and etc. to denote the expectation, variance, and etc. with respect to $\mathsf{P}_{\theta}$. The data is randomly divided into an estimation sample and a training (equivalently, nuisance) sample of size $n = N / 2$. To avoid notational clutter, all expectations, variances, and probabilities are conditional on the training sample unless otherwise stated. For a (random) vector $V$, $\Vert V \Vert_{\theta} \equiv \mathsf{E}_{\theta} [V^{\otimes 2}]^{1/2} = \mathsf{E}_{\theta} [V^{\top} V]^{1/2}$ denotes its $L_{2} (\mathsf{P}_{\theta})$ norm conditioning on the training sample, $\Vert V \Vert \equiv (V^{\otimes 2})^{1/2} = (V^{\top} V)^{1/2}$ its $\ell_{2}$ norm and $\Vert V \Vert_{\infty}$ its $L_{\infty} (\mathsf{P}_{\theta})$ norm. For any matrix $M$, $\Vert M \Vert$ will be reserved for its operator norm. Given an integer $k$, and a random vector $\bar{\mathsf{z}}_{k} = \bar{\mathsf{z}}_{k} (X)$, $\Pi_{\theta} [\cdot | \bar{\mathsf{z}}_{k}]$ denotes the population linear projection operator onto the linear space spanned by $\bar{\mathsf{z}}_{k}$ conditioning on the training sample, and $\Pi_{\theta}^{\perp} [\cdot | \bar{\mathsf{z}}_{k}] = \left( \mathsf{I} - \Pi_{\theta} \right) \left[ \cdot | \bar{\mathsf{z}}_{k} \right]$ is the projection onto the ortho-complement of $\bar{\mathsf{z}}_{k}$ in the Hilbert space $L_{2} (\mathsf{P}_{F})$ where $\mathsf{P}_{F}$ denotes the marginal law of $X$ under $\theta$. That is, for a random variable $V$,

The following common asymptotic notations are used throughout the paper: $x \lesssim y$ (equivalently $x = O (y)$ or $y = \Omega (x)$) denotes that there exists some constant $C > 0$ such that $x \leqslant C y$, $x \asymp y$ (equivalently $x = \Theta (y)$) means there exist some constants $c_{1} > c_{2} > 0$ such that $c_{2} |y| \leqslant |x| \leqslant c_{1} |y|$. $x = o (y)$ or $y = \omega (x)$ or $y \gg x$ or $x \ll y$ is equivalent to $\lim_{x, y \rightarrow \infty} \frac{x}{y} = 0$. For a random variable $X_{n}$ with law $\mathsf{P}$ possibly depending on the sample size $n$, $X_{n} = O_{\mathsf{P}} (a_{n})$ denotes that $X_{n} / a_{n}$ is bounded in $\mathsf{P}$-probability, and $X_{n} = o_{\mathsf{P}} (a_{n})$ means that $\lim_{n \rightarrow \infty} \mathsf{P} (|X_{n} / a_{n} | \geqslant \epsilon) = 0$ for every positive $\epsilon$.

Influence functions of DR functionals

\leavevmode

In this section, we generalize our discussion in the Introduction on $\psi (\theta) = - \mathsf{E}_{\theta} [Y (1)]$ under ignorability to the DR functionals and provide a more detailed exposition. We first review the influence functions of DR functionals (see Definition (ref)) established in rotnitzky2021characterization. Although the definition of DR functional is quite abstract, rotnitzky2021characterization derived the (nonparametric) influence functions of DR functionals by directly leveraging the form of the product bias appeared in Definition (ref). We summarize their results below for the sake of completeness.

We first give some examples of the influence functions $\mathsf{IF}_{1,\psi} (\theta)$ and Riesz representers of DR functionals. Similar moment conditions such as (ref) can be traced back to newey1994asymptotic, newey1997convergence, newey1998undersmoothing, newey2004twicing, chernozhukov2022automatic.

\allowdisplaybreaks We now briefly comment on the three parts of Proposition (ref) as they are all important for future development of the paper. Eq. (ref) exhibits the general formula of influence functions of DR functionals, which are the basic building blocks for standard DML estimators and most nonstandard DML estimators. Also, our framework heavily relies on higher order influence functions, which are derived from the (first-order) influence functions. Eq. (ref) and (ref) together suggest natural loss functions that could be used to fit the nuisance functions $b$ and $p$ from data. To see why, we first make the following important observation:

The proof of Lemma (ref) can be found in Appendix (ref). Based on Lemma (ref), one can often establish rates of convergence of $\widehat{b}, \widehat{p}$ to $b, p$ in $L_{2} (\mathsf{P}_{\theta})$ norm by solving the following minimization problem from the training sample:

where $\mathcal{F}$ is the set of functions computable by the machine learning algorithm. This is because the convergence properties of $\widehat{b}$ and $\widehat{p}$ are often established by excess risk bound that connects the empirical loss (ref) to the expected $\lambda$-weighted $L_{2} (\mathsf{P}_{\theta})$-loss (ref) under certain complexity-reducing assumptions. Importantly, given positive weight functions $w_{b}, w_{p}$ over $X$, that are strictly bounded from above and below, any $w_{b}$- (or $w_{p}$-) weighted $L_{2} (\mathsf{P}_{\theta})$-loss and $\lambda$-weighted $L_{2} (\mathsf{P}_{\theta})$-loss are equivalent up to constants under our assumption on $\lambda$ in Definition (ref):

and similarly for $\mathsf{E}_{\theta} [\lambda (X) \{\widehat{p} (X) - p (X)\}^{2}]$. In Appendix (ref), we also point out several other possible choices, including unweighted $L_{2} (\mathsf{P}_{\theta})$-loss, cross-entropy loss, and adversarial loss minimization strategies, that are also equivalent to the $\lambda$-weighted $L_{2} (\mathsf{P}_{\theta})$-loss up to multiplicative or additive constants.

For DR functionals, $\mathsf{CSBias}_{\theta} (\widehat{\psi}_{1})$, the Cauchy-Schwarz bias of $\widehat{\psi}_{1}$, is defined as

which is an upper bound of (the absolute value of) $\mathsf{Bias}_{\theta} (\widehat{\psi}_{1})$ to be defined in (ref), by CS inequality. Hence we have

that is, they are equivalent up to a multiplicative positive constant.

Standard DML estimators

\leavevmode

The following algorithm defines the standard DML estimators $\widehat{\psi}_{1}$ and $\widehat{\psi}_{\mathsf{cf}, 1}$ of a general DR functional $\psi (\theta)$ satisfying the regularity conditions given in Proposition (ref).

The next proposition provides asymptotic properties of standard DML estimators $\widehat{\psi}_{1}$ (and $\widehat{\psi}_{\mathsf{cf}, 1}$) of DR functionals. Its proof is straightforward and can be found in chernozhukov2018double, smucler2019unifying or farrell2021deep.

\allowdisplaybreaks

A prelude to our approach

\leavevmode

As briefly described in Section (ref), our main goal is to develop assumption-lean empirical methods to falsify $\mathsf{NH}_{0, \mathsf{CS}}: \mathsf{CSBias}_{\theta} (\widehat{\psi}_{1}) = o (n^{- 1 / 2})$ by falsifying its operationalized pair $\mathsf{H}_{0, \mathsf{CS}} (\delta): \mathsf{CSBias}_{\theta} (\widehat{\psi}_{1}) < \mathsf{s.e.}_{\theta} (\widehat{\psi}_{1}) \delta$.

However, as we argued in the Introduction, $\mathsf{CSBias}_{\theta} (\widehat{\psi}_{1})$ and $\mathsf{Bias}_{\theta} (\widehat{\psi}_{1})$ are generally not uniformly consistently estimable without complexity-reducing assumptions on $b$ and $p$ or their estimates $\widehat{b}$ and $\widehat{p}$. This is evident from the following decomposition of $\mathsf{Bias}_{\theta} (\widehat{\psi}_{1})$ because we do not have control over $\Pi_{\theta}^{\perp} [\lambda^{1 / 2} (\widehat{b} - b) | \lambda^{1 / 2} \bar{\mathsf{z}}_{k}]$ and $\Pi_{\theta}^{\perp} [\lambda^{1 / 2} (\widehat{p} - p) | \lambda^{1 / 2} \bar{\mathsf{z}}_{k}]$ without restrictive complexity-reducing assumptions on $b$ and $p$ or on $\widehat{b}$ and $\widehat{p}$:

where $\mathsf{TB}_{\theta, k} (\widehat{\psi}_{1}) \coloneqq \mathsf{E}_{\theta} [\Pi_{\theta}^{\perp} [\lambda^{1 / 2} (\widehat{b} - b) | \lambda^{1 / 2} \bar{\mathsf{z}}_{k}] (X) \Pi_{\theta}^{\perp} [\lambda^{1 / 2} (\widehat{p} - p) | \lambda^{1 / 2} \bar{\mathsf{z}}_{k}] (X)]$, which was referred to as the truncation bias in robins2008higher, liu2020nearly (also see Appendix (ref)), and

Here $\bar{\mathsf{z}}_{k}$ is a vector of $k$-dimensional vector of dictionary chosen by the analyst satisfying mild regularity conditions (see Condition (ref) in Section (ref)) and $\Sigma_{k} \coloneqq \mathsf{E}_{\theta} [\lambda (X) \bar{\mathsf{z}}_{k} (X) \bar{\mathsf{z}}_{k} (X)^{\top}] \equiv \mathsf{E}_{\theta} [S_{bp} \bar{\mathsf{z}}_{k} (X) \bar{\mathsf{z}}_{k} (X)^{\top}]$\footnote{To relate to the discussion in Section (ref), $S_{bp} = A$ for the functional $\psi (\theta) = - \mathsf{E}_{\theta} [Y (1)]$.} is the population Gram matrix of $\lambda^{1 / 2} \bar{\mathsf{z}}_{k}$.

Since we generally know neither the sign nor the magnitude of $\mathsf{TB}_{\theta, k} (\widehat{\psi}_{1})$, tests of $\mathsf{H}_{0} (\delta)$ based on $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$ fail to protect the nominal level uniformly. But fortunately, as can be shown in the same fashion as in (ref), the quantity $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$ is a lower bound of $\mathsf{CSBias}_{\theta} (\widehat{\psi}_{1})$, making it possible to construct nominal-level tests of $\mathsf{H}_{0, \mathsf{CS}} (\delta)$. The test statistic of $\mathsf{H}_{0, \mathsf{CS}} (\delta)$, introduced in Section (ref), is based on estimators of $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$, derived using the theory of higher-order influence functions (HOIFs) (which are higher-order $U$-statistics) developed in a series of papers by some of the authors robins2008higher, robins2017minimax, liu2017semiparametric. Due to space limitation, we refer the interested readers to robins2008higher or van2014higher for a more comprehensive review. {\it En route} to constructing these HOIF estimators and associated tests, we require access to only the study data and the functions $\widehat{b}$ and $\widehat{p}$ obtained by analysts. However, our tests are constructed without: i) refitting, modifying, or even having knowledge of the machine learning algorithms that have been employed to compute $\widehat{b}, \widehat{p}$ from the training sample, and ii) requiring any assumptions at all (aside from a few standard, quite weak assumptions given later) -- in particular, without making any assumptions about the smoothness or sparsity of the nuisance functions $b$ or $p$. The key to achieve i) is by conditioning on the training sample data. For a related discussion, see Remark (ref).

Formally, our proposed approach begins by noticing that $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$ can be unbiasedly estimated by the following infeasible second order $U$-statistic

with standard error of order $O \left( \frac{\sqrt{k}}{n} \vee \frac{1}{\sqrt{n}} \right)$ (see Theorem (ref) for more details). In fact, $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$ is the second-order influence function of $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$\footnote{This is why we adopt the notation $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$ following the convention in robins2008higher.}; see Appendix (ref) for more details. $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$ is infeasible because in general $\Sigma_{k}^{-1}$ is unknown. The unbiasedness follows directly from equation (ref) in Proposition (ref); the proof of the bound on the standard error of $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$ is deferred to Appendix (ref), and can also be found in liu2020nearly, which mostly focused on the infeasible case by assuming $\Sigma_{k}^{-1}$ to be known. Since $\Sigma_{k}^{-1}$ is generally unknown, we propose to estimate $\Sigma_{k}^{-1}$ by $\widehat{\Sigma}_{k}^{-1}$ from the training sample data, where $\widehat{\Sigma}_{k} \coloneqq \mathsf{P}_{n_{\mathsf{tr}}} [S_{b p} \bar{\mathsf{z}}_{k} (X) \bar{\mathsf{z}}_{k} (X)^{\top}]$ is simply the empirical Gram matrix estimator. In particular, throughout this paper, we choose $1 \ll k \ll n$. The reason for this choice is discussed in detail in the following remark.

Main results: A valid falsification test of null hypothesis $\mathsf{H}_{0, \mathsf{CS}} (\delta)$

Having introduced $\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$ as an estimator of $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$, we are now prepared to construct a valid nominal level-$\alpha^{\dag}$ falsification test of the non-asymptotic null hypothesis paired with $\mathsf{NH}_{0, \mathsf{CS}}$:

What constitutes a valid test of $\mathsf{H}_{0, \mathsf{CS}} (\delta)$?

\leavevmode

Before proceeding to the main result, we note that the infeasible falsification test of $\mathsf{H}_{0, \mathsf{CS}} (\delta)$ proposed in liu2020nearly can be generalized to the class of DR functionals as follows\footnote{We also provide relevant results in Appendix (ref).}:

where $\widehat{\mathsf{s.e.}} [\widehat{\psi}_{1}]$ (see (ref)) and $\widehat{\mathsf{s.e.}} [\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})]$ (see (ref) in Appendix (ref)) are estimators of the standard errors $\mathsf{s.e.}_{\theta} [\widehat{\psi}_{1}]$ and $\mathsf{s.e.}_{\theta} [\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})]$. Here one choose the cut-off $\varsigma_{k} = z_{\alpha^{\dag} / 2}$ by normal approximation. The asymptotic validity of $\widehat{\chi}_{2, k} (\Sigma_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta)$ relies on the unbiasedness of $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$, but the unbiasedness can be relaxed as in the proposition below.

We did not give a formal proof as the argument for Proposition (ref) is quite simple. Since $\mathsf{s.e.}_{\theta} [\widehat{\mathsf{Bias}}_{\theta, k} (\widehat{\psi}_{1})]$ can be as small as of order $\frac{\sqrt{k}}{n}$, one needs the bias $\mathsf{E}_{\theta} [\widehat{\mathsf{Bias}}_{\theta, k} (\widehat{\psi}_{1}) - \mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})]$ to be dominated by this order under the null $\mathsf{H}_{0, \mathsf{CS}} (\delta)$ or the null $\mathsf{H}_{0, k} (\delta)$ to protect the level of the test. The final clause of Proposition (ref) will be most relevant to Section (ref). We also refer readers to Appendix (ref) for more details.

Bias-reduced estimator of $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$

\leavevmode

The most natural estimator of $\widehat{\mathsf{Bias}}_{\theta, k} (\widehat{\psi}_{1})$ is $\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$ and the corresponding test is

However, we {\it cannot} prove that the above test is asymptotically valid for $\mathsf{H}_{0, \mathsf{CS}} (\delta)$, because Theorem (ref) below will show the bias of $\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$ for estimating $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$ is

under $\mathsf{H}_{0, \mathsf{CS}} (\delta)$. This obtained upper bound $O \left( \frac{\sqrt{k \mathsf{log} k}}{n} \right)$ of the bias due to estimating $\Sigma_{k}^{-1}$ exceeds that which is needed to protect the level of the test $\widehat{\chi}_{2, k} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta)$, as given in Condition (1) of Proposition (ref).

Fortunately, Theorem (ref) will also show that the following third-order $U$-statistic estimator of $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$

can reduce the upper bound of the bias due to estimating $\Sigma_{k}^{-1}$ from $O \left( \frac{\sqrt{k \mathsf{log} k}}{n} \right)$ in (ref) to

as long as we choose $k$ such that $k \mathsf{log} k = o (n)$.

To avoid clutter in the remainder of this paper, we introduce the following additional notation for various $L_{q}$-type norms of functions or their weighted $L_{2} (\mathsf{P}_{\theta})$-projections, for $q = 2, 4$:

which will also be useful for the later development of this paper. Note that we define $\mathsf{CSBias}_{\theta, k} (\widehat{\psi}_{1}) \equiv \mathbb{L}_{\theta, 2, \widehat{b}, k} \mathbb{L}_{\theta, 2, \widehat{p}, k}$, generalizing the definition of $\mathsf{CSBias}_{\theta, k} (\widehat{\psi}_{1})$ given in (ref) to the whole class of DR functionals.

We also need to further impose the following weak regularity conditions (Condition (ref)).

Note that Condition (ref) should also be compared to the following slightly stronger Condition (ref) in liu2020nearly, which is the same as Condition (ref) except (3) shall be replaced by the following (3'):

Finally, we summarize the statistical properties of $\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$ and $\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$ in Theorem (ref) below:

The proposed assumption-lean falsification test of $\mathsf{H}_{0, \mathsf{CS}} (\delta)$

\leavevmode

All the previous discussions in this section culminate in (1) the following test of $\mathsf{H}_{0, \mathsf{CS}} (\delta)$:

and (2) the following theorem showing its asymptotically validity:

The proof of Theorem (ref) is a direct consequence of Theorem (ref) and is deferred to Appendix (ref). Since $\widehat{\chi}_{3, k} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta)$ is not a consistent test, we say $\widehat{\chi}_{3, k} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta)$ can only falsify the null hypothesis $\mathsf{H}_{0, \mathsf{CS}} (\delta)$.

What if nonstandard DML estimators are used in practice?

Motivation for testing $\mathsf{H}_{0, k} (\delta)$

\leavevmode

Several recent works newey2018cross, bradic2019minimax, kline2020leave, kennedy2020towards, mcgrath2022undersmoothing have exhibited nonstandard DML estimators, denoted also as $\widehat{\psi}_{1}$ in this section (to avoid introducing new notation at this point), that improve upon standard DML estimators, with $\mathsf{Bias}_{\theta} (\widehat{\psi}_{1}) = o (n^{-1/2})$, even though $\mathsf{CSBias}_{\theta} (\widehat{\psi}_{1}) \gg n^{-1/2}$, i.e. $\kappa_{b} + \kappa_{p} \leqslant 1 / 2$ [recall that $n^{-\kappa_{b}}$ and $n^{-\kappa_{p}}$ are rates of convergence of $\widehat{b}$ and $\widehat{p}$ to $b$ and $p$ in (weighted) $L_{2} (\mathsf{P}_{\theta})$ norm]. However, to obtain these results, they all (1) use very special nuisance function estimators with $\widehat{b}$ and $\widehat{p}$ computed from separate non-overlapping subsamples of the training sample, and (2) assume very specific complexity-reducing assumptions such as H\"{o}lder smoothness newey2018cross, mcgrath2022undersmoothing or (approximate) sparsity bradic2019minimax. The problem for our methodology when applied to such estimators is that, if (1) and (2) above were true, our test $\widehat{\chi}_{3, k} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta)$ might still reject because $\kappa_{b} + \kappa_{p} \leqslant 1 / 2$, even though $\mathsf{Bias}_{\theta} (\widehat{\psi}_{1}) = o (n^{-1/2})$ holds.

In this section, we study one possible approach to extend our methodology to apply to nonstandard DML estimators. We will consider an approach in which we simply wish to test the null hypothesis $\mathsf{H}_{0, k} (\delta)$. This can be viewed as testing if the bias of $\widehat{\psi}_{1}$ is large along the {\it direction} of the basis/dictionary $\bar{\mathsf{z}}_{k}$ chosen by the analyst. If the analyst has a good grasp of the direction along which $\mathsf{Bias}_{\theta} (\widehat{\psi}_{1})$ may be large, say when $b, p$ belong to H\"{o}lder balls with certain smoothness, then she can test $\mathsf{H}_{0, k} (\delta)$ with $\bar{\mathsf{z}}_{k}$ chosen to be wavelet or B-spline series. If otherwise, we then suggest the analyst tests $\mathsf{H}_{0, k} (\delta)$ with a variety of choices of $\bar{\mathsf{z}}_{k}$, as a form of sensitivity analyses robins2003general. If $\mathsf{H}_{0, k} (\delta)$ is rejected with some $\bar{\mathsf{z}}_{k}$ (with a very small p-value), then there is strong evidence that $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1}) \gtrsim n^{-1/2}$ under the operationalized pairing between $\mathsf{H}_{0, k} (\delta)$ and $\mathsf{NH}_{0, k}: \mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1}) = o (n^{-1/2})$. Recall from Section (ref) that $\mathsf{Bias}_{\theta} (\widehat{\psi}_{1})$ can be decomposed as follows:

Furthermore, $\mathsf{TB}_{\theta, k} (\widehat{\psi}_{1})$ cannot be uniformly consistently estimated under the assumption-lean model being considered here. Even so, $\mathsf{Bias}_{\theta} (\widehat{\psi}_{1})$ will not be $o (n^{-1/2})$, and therefore its associated Wald CI centered at $\widehat{\psi}_{1}$ will not be valid, unless $\mathsf{TB}_{\theta, k} (\widehat{\psi}_{1})$ and $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$ happen to have leading order terms of the same magnitude but opposite signs. As it seems quite fortuitous for such a cancellation to occur, an analyst might agree to retract her claim that the Wald CI is valid. Hence in this section we focus on $\mathsf{H}_{0, k} (\delta)$ and investigate how to test this particular null hypothesis.

Issues with $\widehat{\chi}_{3, k} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta)$ as a test of $\mathsf{H}_{0, k} (\delta)$ and a rescue by HOIFs

\leavevmode

In this section, we begin by showing that $\widehat{\chi}_{3, k} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta)$ is not necessarily an $\alpha^{\dag}$-level test of $\mathsf{H}_{0, k} (\delta)$ for nonstandard DML estimators (neither for standard DML estimators, which does not contradict the claim about the validity of $\widehat{\chi}_{3, k} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta)$ in Theorem (ref) because it is shown to be valid under $\mathsf{H}_{0, \mathsf{CS}} (\delta)$, a different null hypothesis from $\mathsf{H}_{0, k} (\delta)$). Recall from Theorem (ref) that the bias $\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$ for $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$ was upper bounded by:

Under $\mathsf{H}_{0, k} (\theta)$, $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1}) \lesssim n^{- 1 / 2}$, but because $\mathbb{L}_{\theta, 2, \widehat{b}, k} \mathbb{L}_{\theta, 2, \widehat{p}, k} \equiv \mathsf{CSBias}_{\theta, k} (\widehat{\psi}_{1}) \geqslant \mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$ by CS inequality, we cannot ensure $\mathbb{L}_{\theta, 2, \widehat{b}, k} \mathbb{L}_{\theta, 2, \widehat{p}, k} \lesssim n^{-1/2}$. Hence one cannot guarantee $\mathsf{EB}_{3, \theta, k} = o \left( \frac{\sqrt{k}}{n} \right)$, as required in Proposition (ref) to ensure the validity of $\widehat{\chi}_{3, k} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta)$ as a test of $\mathsf{H}_{0, k} (\delta)$, without making further unverifiable assumptions on $b, p, \widehat{b}, \widehat{p}$ or $\bar{\mathsf{z}}_{k}$. A natural solution would be to further reduce the bias of $\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$ due to estimating $\Sigma_{k}^{-1}$ by $\widehat{\Sigma}_{k}^{-1}$, which motivates the following $m$-th order $U$-statistic estimator of $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$:

$\widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1})$ can further reduce the bias [and $\widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\Sigma_{k}^{-1})$ is the $m$-th order influence function of $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$ (see HOIF-related theory in Appendix (ref))]. In particular, we have the following theorem on the statistical properties of $\widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1})$, similar to Theorem (ref). Again, the bias and variance bounds under Condition (ref) are proved in liu2017semiparametric. Under Condition (ref), see the proof of Lemma (ref) in Appendix (ref). Recall from Section (ref) that Condition (ref) is weaker than Condition (ref).

The proposed assumption-lean test of $\mathsf{H}_{0, k} (\delta)$ and early-stopping

\leavevmode

Based on the above discussion, we propose the following nominal $\alpha^{\dag}$-level test of $\mathsf{H}_{0, k} (\delta)$: for any fixed $m \geqslant 3$,

Theorem (ref) immediately implies the following:

Formally, for a fixed integer $M > 0$ that is determined by the analyst's computational budget, define

If $\widehat{\chi}_{M, k}^{es} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta)$ fails to reject $\mathsf{H}_{0, k} (\delta)$ at some $m \leqslant M$, we stop the test at $m$ and claim that we fail to reject $\mathsf{H}_{0, k} (\delta)$. This “early stopping” procedure is again an asymptotically valid $\alpha^{\dag}$-level test, under slightly stronger assumptions than those in Theorem (ref):

The power of the above “early-stopping” procedure is more challenging to characterize, which we leave as future work.

Monte Carlo experiments

In the Monte Carlo (MC) experiments, we focus on the parameter (up to a minus sign) that was extensively discussed in the Introduction, $\psi (\theta) \equiv \mathsf{E}_{\theta} [Y (a = 1)] \equiv \mathsf{E}_{\theta} [b (X)]$ under ignorability. Recall that $p(X) = 1 / \mathsf{E}_{\theta} [A | X]$ is the inverse propensity score and $b (X) = \mathsf{E}_{\theta} [Y | A = 1, X]$ is the conditional mean of the outcome in the treatment group. For this parameter, $\Sigma_{k} = \mathsf{E}_{\theta} [A \bar{\mathsf{z}}_{k}(X) \bar{\mathsf{z}}_{k}(X)^{\top}]$ and $\widehat{\Sigma}_{k} = n^{-1} \sum_{i \in \mathsf{tr}} A_{i} \bar{\mathsf{z}}_{k} (X_{i}) \bar{\mathsf{z}}_{k} (X_{i})^{\top}$. We choose $\psi (\theta) \equiv 0$.

We consider two simulation setups. In simulation setup I, we draw $N=100,000$ i.i.d. $X_{j}$ for $j = 1, \ldots, 4$ (so $d = 4$). The marginal density $f_{j}$ of each $X_{j}$ is supported on $[0, 1]$ with $f_{j} \in \text{H\"{o}lder} (0.1 + c)$ for some small $c > 0$, as defined in Appendix (ref). The correlation between each pair of $X_{j}$ and $X_{k}$, with $j \neq k$, is introduced based on the algorithm described in Appendix (ref). We then simulate $Y$ and $A$ according to the following data generating mechanism:

and

where $h_{b} (\cdot; 0.25)$ and $h_{p} (\cdot; 0.25)$ have the forms as defined in Appendix (ref) and hence both belong to $\text{H\"{o}lder} (0.25 + c)$ for some very small $c > 0$. The numerical values for $\left( \tau_{b, j}, \tau_{p, j} \right)_{j = 1}^{4}$ are provided in Table (ref). We fix half of the $N = 100,000$ samples as the training sample so $n_{\mathsf{tr}} = n = 50,000$ and only consider the randomness from the estimation sample in the simulation. In simulation setup II, we consider the same data generating mechanism as in setup I except that we choose $b (X) = \sum_{j = 1}^{4} \tau_{b, j} h_{b} (X_{j}; 0.6)$ and $1 / p (X) \equiv \mathsf{expit} \left\{ \sum_{j = 1}^{4} \tau_{p, j} h_{p} (X_{j}; 0.6) \right\}$ where $h_{b} (\cdot; 0.6)$ and $h_{p} (\cdot; 0.6)$ have the forms as defined in Appendix (ref) and hence both belong to $\text{H\"{o}lder} (0.6 + c)$ for some very small $c > 0$.

We choose D12 (or equivalently db6) Daubechies wavelets at resolutions $\ell \in (6, 7, 8)$ to form the dictionary

with the corresponding $k' \in \{2^{6} = 64, 2^{7} = 128, 2^{8} = 256\}$ and $k \in \{64 \cdot 4 = 256, 128 \cdot 4 = 512, 256 \cdot 4 = 1024\}$. To compute the oracle statistics and tests, we evaluate $\Sigma _{k}$ through MC integration by simulating $L = 10^{7}$ independent $(A, X)$ from the true data generating law. To investigate the finite sample performance of the statistical procedures developed in this article, all the summary statistics of the MC experiments are calculated based on 100 replicates. We estimate the nuisance functions $1 / p (x)$ and $b (x)$ using generalized additive models (GAMs). In particular, the smoothing parameters were selected by generalized cross validation, the default setup in $\mathsf{gam}$ function from $\mathsf{R}$ package $\mathsf{mgcv}$. We choose db6 father wavelets to construct HOIF estimators when analyzing the simulated data.

Finite sample performance of tests for $\mathsf{H}_{0, k} (\delta)$

In this section, we consider testing the null hypothesis $\mathsf{H}_{0, k} (\delta)$. Henceforth we investigate the finite sample performance of the oracle statistics and tests $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$, $\widehat{\psi}_{2, k} (\Sigma_{k}^{-1})$, and $\widehat{\chi}_{2, k} (\Sigma_{k}^{-1})$ where $\widehat{\psi}_{2, k} (\Sigma_{k}^{-1}) = \widehat{\psi}_{1} - \widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$, together with the statistics and tests relying on $\widehat{\Sigma}_{k}^{-1}$: $\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$, $\widehat{\mathbb{IF}}_{33, k} (\widehat{\Sigma}_{k}^{-1})$, $\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$, $\widehat{\chi}_{3, k} (\widehat{\Sigma}_{k}^{-1})$, $\widehat{\psi}_{2, k} (\widehat{\Sigma}_{k}^{-1})$, $\widehat{\psi}_{3, k} (\widehat{\Sigma}_{k}^{-1})$, and $\widehat{\chi}_{2, k} (\widehat{\Sigma}_{k}^{-1})$, where $\widehat{\psi}_{2, k} (\widehat{\Sigma}_{k}^{-1}) = \widehat{\psi}_{1} - \widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$ and $\widehat{\psi}_{2, k} (\widehat{\Sigma}_{k}^{-1}) = \widehat{\psi}_{1} - \widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$

First, we check the asymptotic normalities of $\frac{\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})}{\widehat{\mathsf{s.e.}} [\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})]}$, $\frac{\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})}{\widehat{\mathsf{s.e.}} [\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})]}$, $\frac{\widehat{\mathbb{IF}}_{33, k} (\widehat{\Sigma}_{k}^{-1})}{\widehat{\mathsf{s.e.}} [\widehat{\mathbb{IF}}_{33, k} (\widehat{\Sigma}_{k}^{-1})]}$, and $\frac{\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})}{\widehat{\mathsf{s.e.}} [\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})]}$ through normal qq-plots displayed in Figures (ref) (simulation setup I) and (ref) (simulation setup II). We observe that the distributions of most of these statistics are close to normal at different $k$'s ($k = 256$: left panels; $k = 512$: middle panels; $k = 1024$: right panels).

In the simulation, we use nonparametric bootstrap to estimate the standard errors of $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$, $\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$, $\widehat{\mathbb{IF}}_{33, k} (\widehat{\Sigma}_{k}^{-1})$ and $\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$, as described in Appendix (ref). We study if the estimated standard errors of $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1}) $, $\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$, $\widehat{\mathbb{IF}}_{33, k} (\widehat{\Sigma}_{k}^{-1})$, and $\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$ by nonparametric bootstrap are close to their true standard errors (calibrated by the MC standard deviations from 100 replicates in the simulation). We use $B = 100$ bootstrap samples to compute the bootstrapped standard errors as the estimated standard errors for all four statistics. In Tables (ref) (simulation setup I) and (ref) (simulation setup II), we display the MC standard deviations (the upper numerical values in each cell), accompanied with the MC averages of the estimated standard errors (the lower numerical values outside the parenthesis in each cell) and MC standard deviations of the estimated standard errors (the lower numerical values inside the parenthesis in each cell) of all three statistics at $k = 256$ (left panel), $k = 512$ (middle panel) and $k = 1024$ (right panel). From Table (ref) and Table (ref), we observe that the estimated standard errors only slightly differ from the MC standard deviations.

Then we investigate (1) the finite sample performance of $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$, $\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$, and $\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$ and evaluate how close they are to $\mathsf{Bias}_{\theta} (\widehat{\psi}_{1})$, which is evaluated based on the MC bias of 100 replicates, and (2) the rejection rate of the tests $\widehat{\chi}_{2, k} (\Sigma_{k}^{-1}; z_{0.10 / 2}, \delta)$, $\widehat{\chi}_{2, k} (\widehat{\Sigma}_{k}^{-1}; z_{0.10 / 2}, \delta)$ and $\widehat{\chi}_{3, k} (\Sigma_{k}^{-1}; z_{0.10 / 2}, \delta)$ for the null hypothesis $\mathsf{H}_{0, k} (\delta)$. The numerical results are shown in Tables (ref) (simulation setup I) and (ref) (simulation setup II).

Concluding remarks

In this paper, we developed a valid assumption-lean test, based on third-order $U$-statistics, that can empirically falsify the justification for a Wald CI centered at a standard DML estimator $\widehat{\psi}_{1}$ covering the underlying DR functional $\psi (\theta)$ at the nominal rate (see Section (ref)). When nonstandard DML estimators are used, we also develop a test that is based on higher-order $U$-statistics, which are in fact HOIFs of $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$. We mention a few interesting future directions to end our manuscript. First, it is important to extend the proposed approach to allow endogeneity angrist1996identification, ai2003efficient, chen2005measurement, newey2003instrumental, ai2007estimation, chen2013optimal, breunig2019simple, chen2018optimal, chen2016methods, when some auxiliary variables are available for point-identifying the causal effects. Machine learning, including deep learning, has also been applied to these problems in recent years chen2023efficient, kompa2022deep. In Appendix (ref), we document the HOIFs for the running example $\psi (\theta) = - \mathsf{E}_{\theta} [Y (a = 1)]$ in Section (ref) under the so-called “proximal causal learning” framework tchetgen2020introduction. This framework has been shown to be closely related to other approaches dealing with endogeneity, e.g. synthetic controls abadie2010synthetic, shi2021theory or quadratic functionals of nonparametric instrumental variable regression breunig2019simple. Based on the form of these HOIFs, it is straightforward to generalize our method to scenarios under endogeneity. Second, as mentioned earlier in Literature Overview, parameters implicitly defined via (conditional) moment restrictions involving nonparametric nuisance functions include, as special cases, DR functionals and certain parameters related to instrumental variables and proximal causal inference ai2003efficient, ai2007estimation, ai2012semiparametric. Developing the theory of higher-order influence functions for these parameters will be a natural next step. Third, extending our framework to heterogeneous treatment effect chernozhukov2017generic, kennedy2022minimax could also be of practical interest, given its significance in personalized decision making. Fourth, another interesting direction is to explore if it is possible to use other bias correction strategies to construct the assumption-lean falsification test, such as the (iterative) bootstrap approach investigated in cattaneo2018kernel, cattaneo2019two, and koltchinskii2020estimation. Finally, we are also considering to directly learn a data-driven representation $\widetilde{\bar{\mathsf{z}}}_{k}$ via the penultimate layer ansuini2019intrinsic, damian2022neural or distillation ha2021adaptive of a deep neural network trained to predict the residuals of the nuisance function estimators, hoping to increase the chance of rejection when the bias or the CS bias indeed exceeds $n^{-1/2}$ in real-world settings.

\printbibliography

\newgeometry{margin = 0.5in, paperwidth=11in, paperheight=12in} \eject \pdfpagewidth=11in \pdfpageheight=12in \oddsidemargin +0.2in \evensidemargin +0.0in \topmargin 5pt \linespread{1.5}\parskip .05in

\allowdisplaybreaks