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

The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.

143,443 characters

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


\newrefsection

\title[Assumption-lean higher-order testing]{\normalfont A\lowercase{ssumption-lean Falsification Tests of Rate Double-Robustness of Double-Machine-Learning Estimators}}
\author{$\text{Lin Liu}^{1}$}
\thanks{1. Institute of Natural Sciences, MOE-LSC, School of Mathematical Sciences, CMA-Shanghai, SJTU-Yale Joint Center for Biostatistics and Data Science, Shanghai Jiao Tong University; Shanghai Artificial Intelligence Laboratory \href{[email removed]}{[email removed]}}
\author{$\text{Rajarshi Mukherjee}^{2}$}
\thanks{2. Department of Biostatistics, Harvard T. H. Chan School of Public Health, \href{[email removed]}{[email removed]}}
\author{$\text{James M. Robins}^{3}$}
\thanks{3. Department of Epidemiology and Department of Biostatistics, Harvard T. H. Chan School of Public Health, \href{[email removed]}{[email removed]}}
\thanks{This article is based on a talk given by the last author (JMR) in the ``\href{https://www.nber.org/sites/default/files/2020-08/CEME_Micro_2019.pdf}{2019 Conference for Celebrating Whitney Newey's Contributions to Econometrics}'' held at MIT, and is heavily inspired by works of Whitney Newey. The authors would like to thank three anonymous referees and \href{https://economics.yale.edu/people/xiaohong-chen}{Xiaohong Chen} for constructive comments and \href{https://gaofn.xyz/}{Fengnan Gao} (School of Mathematics and Statistics at the University College Dublin) for discussions that significantly improve the paper. LL was supported by NSFC Grants No.12101397 and 12090024, Shanghai Municipal Science and Technology Major Project No.2021SHZDZX0102, Shanghai Municipal Science and Technology Grants No.21ZR1431000 and 21JC1402900. RM was partially supported by NSF Grant EAGER-1941419. JMR was supported by the U.S. Office of Naval Research grant N000141912446, and National Institutes of Health (NIH) awards R01 AG057869 and R01 AI127271.}
\maketitle

\begin{abstract}
The class of doubly-robust (DR) functionals studied by Rotnitzky et al. (2021) is of central importance in economics and biostatistics. It 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 Chernozhukov et al. (2022b) and (ii) the class of functionals studied by Robins et al. (2008). The present state-of-the-art estimators for DR functionals $\psi$ are double-machine-learning (DML) estimators (Chernozhukov et al., 2018). A DML estimator $\widehat{\psi}_{1}$ of $\psi$ depends on estimates $\widehat{p} (x)$ and $\widehat{b} (x)$ of a pair of nuisance functions $p(x)$ and $b(x)$, and is said to satisfy ``rate double-robustness'' if the Cauchy--Schwarz upper bound of its bias is $o (n^{- 1/2})$. Rate double-robustness implies that the bias is $o (n^{-1/2})$, but the converse is false. Were it achievable, our scientific goal would have been to construct valid, assumption-lean (i.e. no complexity-reducing assumptions on $b$ or $p$) tests of the validity of a nominal $(1 - \alpha)$ Wald confidence interval (CI) centered at $\widehat{\psi}_{1}$. But this would require a test of the bias to be $o (n^{-1/2})$, which can be shown not to exist. We therefore adopt the less ambitious goal of falsifying, when possible, an analyst's justification for her claim that the reported $(1 - \alpha)$ Wald CI is valid. In many instances, an analyst justifies her claim by imposing complexity-reducing assumptions on $b$ and $p$ to ensure ``rate double-robustness''. Here we exhibit valid, assumption-lean tests of $\mathsf{H}_{0}$: ``rate double-robustness holds'', with non-trivial power against certain alternatives. If $\mathsf{H}_{0}$ is rejected, we will have falsified her justification. However, no assumption-lean test of $\mathsf{H}_{0}$, including ours, can be a consistent test. Thus, the failure of our test to reject is not meaningful evidence in favor of $\mathsf{H}_{0}$.
\end{abstract}
{\footnotesize \textbf{Key words:} Econometrics, Causal Inference, Machine Learning, Doubly-Robust Functionals, Higher-Order $U$-Statistics, Higher-Order Influence Functions}

\allowdisplaybreaks

\section{Introduction}
\label{sec: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$ \citep{chernozhukov2018double} (see Section \ref{sec:dml} 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 \textit{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 \citet{robins1997toward}, such tests do not exist; for intuition, see Section \ref{sec:prelude}.

We therefore adopt the following less ambitious, but \textit{partially} achievable goal: Our (new) goal is to construct an empirical assumption-lean test that can falsify an analyst's \textit{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 \textit{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})$. \citet{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 \citet{chen2008semiparametric, chen2015sieve}. This property was termed ``rate double-robustness'' in \citet{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 \textit{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 \eqref{loss_data}) using deep neural networks \citep{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 \text{H\"{o}lder}{} spaces \citep{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 \citet{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$ \citep{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{sec:nonstandard}.

In the remainder of the paper, we restrict the functional $\psi$ to the class of Mixed-Bias or Doubly-Robust (DR) functionals \citep{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 \citet{chernozhukov2022automatic} and (ii) the class of functionals studied by \citet{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 \citet{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{def:dr}. 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{def:dr}.} $\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. \citet{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{thm:drml} 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 \citep{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 \citep{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{sec:dml}. This allows us to distinguish standard DML estimators from, for instance, ``DCDR estimators'' of \citet{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 \citep{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
\begin{equation}
|\mathsf{Bias} (\widehat{\psi}_{1})| \leqslant \{\mathsf{E} [p^{-1} (X) (\widehat{b} (X) - b (X))^{2}]\}^{1 / 2} \{\mathsf{E} [p^{-1} (X) (\widehat{p} (X) - p (X))^{2}]\}^{1 / 2}.
\label{cs-heuristic}
\end{equation}
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
\begin{equation*}
\mathsf{CSBias}^{w} (\widehat{\psi}_{1})\equiv \{\mathsf{E} [w_{b} (X) (\widehat{b} (X) - b(X))^{2}]\}^{1/2} \{\mathsf{E} [w_{p} (X) (\widehat{p} (X) - p (X))^{2}]\}^{1/2}.
\end{equation*}
Then applying \text{H\"{o}lder}{} inequality, as in \eqref{holder} 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})$.

\subsection*{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 \textit{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
\begin{equation}
\label{h0-cs}
\mathsf{H}_{0, \mathsf{CS}} (\delta): \mathsf{CSBias} (\widehat{\psi}_{1}) \leqslant \delta \mathsf{s.e.} (\widehat{\psi}_{1})
\end{equation}
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
\begin{equation}
\label{pairing}
\begin{split}
\mathsf{NH}_{0, \mathsf{CS}}: \mathsf{CSBias} (\widehat{\psi}_{1}) = o_{p} (n^{- 1 / 2}) \sim \mathsf{H}_{0, \mathsf{CS}} (\delta): \mathsf{CSBias} (\widehat{\psi}_{1}) < \mathsf{s.e.} (\widehat{\psi}_{1}) \delta,
\end{split}
\end{equation}
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 \eqref{equiv}, 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 \eqref{pairing} 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 \citet{robins1997toward} or \citet{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) \citep{robins2008higher, liu2017semiparametric}] such that one can reject the ``$k$-projected'' null hypothesis
\begin{equation}
\label{h0k}
\mathsf{H}_{0, k} (\delta): |\mathsf{Bias}_{k} (\widehat{\psi}_{1})| \leqslant \delta \mathsf{s.e.} (\widehat{\psi}_{1})
\end{equation}
with nontrivial power. Here $\mathsf{Bias}_{k} (\widehat{\psi}_{1})$ is defined as
\begin{align*}
\mathsf{Bias}_{k} (\widehat{\psi}_{1}) \coloneqq \mathsf{E} [\Pi [p^{-1 / 2} (b - \widehat{b}) | p^{-1 / 2} \bar{\mathsf{z}}_{k}] (X) \Pi [p^{-1 / 2} (p - \widehat{p}) | p^{-1 / 2} \bar{\mathsf{z}}_{k}] (X)],
\end{align*}
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:
\begin{equation}
\label{csbias}
\begin{split}
& \vert \mathsf{Bias}_{k} (\widehat{\psi}_{1}) \vert \leqslant \underbrace{\{\mathsf{E} [\Pi [p^{-1 / 2} (b - \widehat{b}) | p^{-1 / 2} \bar{\mathsf{z}}_{k}] (X)^{2}]\}^{1 / 2} \{\mathsf{E} [\Pi [p^{-1 / 2} (p - \widehat{p}) | p^{-1 / 2} \bar{\mathsf{z}}_{k}] (X)^{2}]\}^{1 / 2}}_{\eqqcolon \mathsf{CSBias}_{k} (\widehat{\psi}_{1})} \leqslant \mathsf{CSBias} (\widehat{\psi}_{1}),
\end{split}
\end{equation}
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
\begin{equation}
\label{h0k-cs}
\mathsf{H}_{0, k, \mathsf{CS}} (\delta): \mathsf{CSBias}_{k} (\widehat{\psi}_{1}) \leqslant \delta \mathsf{s.e.} (\widehat{\psi}_{1}).
\end{equation}
will be explained in Remark \ref{rem:temper} of Section \ref{sec:main}. In principle, though, one can also consider testing $\mathsf{H}_{0, k, \mathsf{CS}} (\delta)$.}\label{ft:temper}.
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 \citet{robins1997toward} or \citet{ritov2014bayesian}. But we also provide an intuitive explanation in Section \ref{sec:prelude}.

To further illustrate our approach, we can rewrite $\mathsf{Bias}_{k} (\widehat{\psi}_{1})$ as follows (see Appendix \ref{app:oracle}):
\begin{align*}
\mathsf{Bias}_{k} (\widehat{\psi}_{1}) = \mathsf{E} [A (Y - \widehat{b} (X)) \bar{\mathsf{z}}_{k} (X)^{\top}] \Sigma_{k}^{-1} \mathsf{E} [\bar{\mathsf{z}}_{k} (X) (A \widehat{p} (X) - 1)]
\end{align*}
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:
\begin{align*}
\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1}) \equiv \frac{1}{n (n - 1)} \sum_{1 \leqslant i_{1} \neq i_{2} \leqslant n} A_{i_{1}} (Y_{i_{1}} - \widehat{b} (X_{i_{1}})) \bar{\mathsf{z}}_{k} (X_{i_{1}})^{\top} \Sigma_{k}^{-1} \bar{\mathsf{z}}_{k} (X_{i_{2}}) (A_{i_{2}} \widehat{p} (X_{i_{2}}) - 1).
\end{align*}
$\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{app:hoif} 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{sec:main} for more details.

A preliminary version of the above idea has appeared in \citet{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.
\begin{itemize}
\item First, we extend the results in \citet{liu2020nearly} to the more realistic case where $\Sigma_{k}$ is unknown. In particular, we construct a test of the null hypothesis of rate double-robustness based on a third-order $U$-statistic (which is the estimated third-order influence function of $\mathsf{Bias}_{k} (\widehat{\psi}_{1})$).
\item Second, the results in this paper require much weaker assumptions\footnote{Note that both the assumptions in this paper and in \citet{liu2020nearly} are easily checkable by the analysts; see Conditions \ref{cond:w} and \ref{cond:sw}.} on the dictionary $\bar{\mathsf{z}}_{k}$ than those in \citet{liu2020nearly}, allowing the proposed methodology to be used more broadly in substantive studies.
\item Third, we extend the results of \citet{liu2020nearly} for expected conditional covariance to the entire DR functional class, which requires that we derive the HOIFs of DR functionals. This could be of independent interest, considering the role of HOIFs in constructing rate-optimal estimators for smooth/differentiable functionals or their generalizations \citep{tchetgen2008minimax, kennedy2022minimax}. In Appendix \ref{app:prox_hoif}, we construct the HOIFs for the average treatment effect in the presence of unmeasured confounding (or under endogeneity) in the setting of proximal causal inference \citep{miao2018identifying, miao2018confounding}. This part is built on \citet{cui2023semiparametric}'s derivation of the first-order influence function under the assumptions underlying proximal causal inference\footnote{Based on personal communication, \href{https://statistics.wharton.upenn.edu/profile/ett/}{Eric Tchetgen Tchetgen}, \href{https://sites.google.com/view/yifancui}{Yifan Cui}, and colleagues have also independently derived HOIFs of this parameter under the proximal causal learning framework in an unpublished manuscript.}. A more thorough discussion is deferred to Section \ref{sec:conclusion} and Appendix \ref{app:prox_hoif}.
\end{itemize}


\subsection*{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 \citet{liu2020nearly}), and has different purposes from specification tests \citep{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 \citet{farrell2015robust, chernozhukov2018double, smucler2019unifying, bradic2019sparsity, chernozhukov2022locally} and references therein.

Some recent works also consider further refinement of standard DML estimators \citep{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 \citep{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 \citet{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{sec:nonstandard}.

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 \citet{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 \citep{ai2003efficient, angrist1996identification, tchetgen2020introduction}. Extending our framework (specifically the theory of higher-order influence functions) to such more complicated parameters is still an open problem.

\subsection*{Organization of the paper}
The remainder of the paper is arranged as follows. In Section \ref{sec:review}, we describe the mathematical setup formally and review the definition and properties of DR functionals recently characterized in \citet{rotnitzky2021characterization}. In Section \ref{sec:dr}, 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{sec:main}, 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 \eqref{if2233}), 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 \citet{robins2008higher}; also see Appendix \ref{app:hoif} for derivations.

In Section \ref{sec:nonstandard}, 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{sec:main} 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{sec:simulation} we present results of simulation studies to evaluate the finite sample performance of our methods. Section \ref{sec:conclusion} concludes with a discussion of some open problems. Many of the technical details are deferred to the Appendix.



\section{Formal setup and an illustration of our approach}
\label{sec:review}
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
\begin{equation*}
\mathcal{M} = \left\{ \mathsf{P}_{\theta}; \theta = (b, p, \theta \setminus \{b, p\}) \in \Theta = \mathcal{B} \times \mathcal{P} \times \Theta \setminus \{\mathcal{B}, \mathcal{P}\} \right\}
\end{equation*}
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 \text{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{rem:w} 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.} \citep{ichimura2022influence} $\mathsf{IF}_{1, \psi} (\theta) = \mathsf{if}_{1, \psi} (O; \theta)$ (and thus a positive and finite semiparametric variance bound \citep{newey1990semiparametric}) and are contained in the mixed bias or doubly-robust class of functionals of \citet{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.

\begin{definition}[Definition of a mixed bias or doubly-robust functional (DR functional) (Definition 1 of \citet{rotnitzky2021characterization})]
\label{def:dr}
$\psi (\theta)$ is a doubly-robust functional if, for each $\theta \in \Theta$ there exists $b: x \mapsto b (x) \in \mathcal{B}$ and $p: x \mapsto p (x) \in \mathcal{P}$ such that (i) $\theta = (b, p, \theta \setminus \{b, p\})$ and $\Theta = \mathcal{B} \times \mathcal{P} \times \Theta \setminus \{\mathcal{B}, \mathcal{P}\}$ and (ii) for any $\theta, \theta'$
\begin{equation}
\psi (\theta) - \psi (\theta^{\prime}) + \mathsf{E}_{\theta} \left[ \mathsf{IF}_{1, \psi} (\theta^{\prime}) \right] = \mathsf{E}_{\theta} \left[ S_{bp} (b (X) - b^{\prime} (X)) (p (X) - p^{\prime} (X)) \right] \label{eq:drbias}
\end{equation}
where $S_{bp} \equiv s_{bp} (O)$ with $o \mapsto s_{bp} (o)$ a known function that does not depend on $\theta$ or $\theta'$ satisfying either $\mathsf{P}_{\theta} (S_{bp} \geqslant 0) = 1$ or $\mathsf{P}_{\theta} (S_{bp} \leqslant 0) = 1$. $b$ and $p$ are called nuisance functions in the semiparametric statistics literature, a terminology we also adopt in this paper. Finally, we let $$\lambda (x) \coloneqq \mathsf{E}_{\theta} [S_{bp} |X = x]$$ and assume $\lambda$ to be strictly bounded from above and below.
\end{definition}

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 \citep{scharfstein1999rejoinder, robins2001comments, bang2005doubly}. That is, by \eqref{eq:drbias}, 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 \citet{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.

\begin{enumerate}
\item The counterfactual mean of $Y$ when $\{0, 1\}$-valued $A$ is set to $1$ (under strong ignorability), $\mathsf{E}_{\theta} [Y (a = 1)]$ is identified from the distribution of $W = (A, A Y, X)$ by the DR functional $\psi^{\dag} (\theta) = \mathsf{E}_{\theta} [b (X)]$ with $b (x) = \mathsf{E}_{\theta} [Y | X = x, A = 1]$, $p (x) = 1 / \mathsf{E}_{\theta} [A | X = x]$, $S_{bp} = - A$, and $\lambda (x) = - p (x)^{-1}$. Here $p(x)$ is the inverse of the propensity score $\pi (x) = \mathsf{E}_{\theta} [A | X = x]$. Because $S_{bp} = - A$, we instead analyze the DR functional
\begin{equation*}
\psi (\theta) = - \psi^{\dag} (\theta) = - \mathsf{E}_{\theta} [b (X)] = - \mathsf{E}_{\theta} [Y (a = 1)]
\end{equation*}
for which $S_{bp} = A$. We will use $\psi (\theta) = - \mathsf{E}_{\theta} [b (X)] = - \mathsf{E}_{\theta} [Y (a = 1)]$ as a running example below\footnote{\citet{chernozhukov2022automatic} used an alternative decomposition under which they included our $A, X$ in their $X$. See Example 1 of \citet{rotnitzky2021characterization}.}. Average treatment effect, one of the most intensively studied causal parameters, is thus a difference of two DR functionals.


\item The expected conditional covariance
\begin{equation*}
\psi (\theta) = \mathsf{E}_{\theta} [(Y - b (X)) (A - p (X))]
\end{equation*}
with $b (x) = \mathsf{E}_{\theta} [Y | X = x], p (x) = \mathsf{E}_{\theta} [A | X = x]$ is a DR functional (equivalently, a MB functional) with $S_{bp} = \lambda (X) \equiv 1$. When $Y$ and $b (X)$ are replaced by $A$ and $p (X)$, $\psi (\theta)$ reduces to the expected conditional variance $\phi (\theta) = \mathsf{E}_{\theta} [(A - p (X))^{2}]$. The expected conditional covariance has been popularized recently by \citet{shah2020hardness} for the problem of conditional independence testing, under a different name ``generalized covariance measure''. The expected conditional covariance also relates to the regression coefficient $\tau$ in the following semiparametric regression: $Y = \tau A + b (X) + \varepsilon$ where $\varepsilon$ is mean zero conditional on $A$ and $X$ \citep{li2011higher}.
\end{enumerate}

To save space, we refer interested readers to \citet{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$,
\begin{equation}
\Pi _{\theta} [V | \bar{\mathsf{z}}_{k}] (x) = \bar{\mathsf{z}}_{k}^{\top} (x) \{\mathsf{E}_{\theta} [\bar{\mathsf{z}}_{k} (X) \bar{\mathsf{z}}_{k} (X)^{\top}]\}^{-1} \mathsf{E}_{\theta} [\bar{\mathsf{z}}_{k} (X) V], \Pi _{\theta}^{\perp} [V | \bar{\mathsf{z}}_{k}] (x) = V - \Pi_{\theta} [V | \bar{\mathsf{z}}_{k}] (x). \label{eq:projection}
\end{equation}

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$.

\subsection{Influence functions of DR functionals}\leavevmode
\label{sec:dr}

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{def:dr}) established in \citet{rotnitzky2021characterization}. Although the definition of DR functional is quite abstract, \citet{rotnitzky2021characterization} derived the (nonparametric) influence functions of DR functionals by directly leveraging the form of the product bias appeared in Definition \ref{def:dr}. We summarize their results below for the sake of completeness.

\begin{proposition}[A summary of Theorems 1 and 2 of \citet{rotnitzky2021characterization}]
\label{thm:if1}
Suppose (1) $\psi (\theta)$ for $\theta = \left( b, p, \theta \setminus \{b, p\} \right) \in \Theta \equiv \mathcal{B} \times \mathcal{P} \times \Theta \setminus \{\mathcal{B}, \mathcal{P}\}$ is a DR functional (equivalently MB functional) according to Definition \ref{def:dr}; (2) the regularity Condition \ref{cond:r} in Appendix \ref{app:regular} (Condition 1 of \citet{rotnitzky2021characterization}) holds; (3) $b (X)$, $p (X)$, $\lambda (X) b (X)$ and $\lambda (X) p (X)$ belong to $L_{2} (\mathsf{P}_{F})$ for all $\theta \in \Theta$.

Then there exists a statistic $S_{0}$ and linear maps $h \mapsto m_{1} (O, h)$ for $h \in \mathcal{B}$ and $h \mapsto m_{2} (O, h)$ for $h \in \mathcal{P}$ independent of $\theta$ satisfying the following:

\begin{itemize}
\item The influence function of $\psi (\theta)$ is given by:
\begin{equation}  \label{eq:if1}
\mathsf{IF}_{1, \psi} (\theta) = \mathcal{H} (b, p) - \psi (\theta)
\end{equation}
where $\mathcal{H} (b, p) \coloneqq S_{bp} b (X) p (X) + m_{1} (O, b) + m_{2} (O, p) + S_{0}$;
\item Furthermore, the following first-order doubly-robust moment conditions hold:
\begin{align}
\label{eq:mean_zero}
\left\{ \begin{array}{c}
\mathsf{E}_{\theta} [S_{bp} h (X) p (X) + m_{1} (O, h)] = 0 \text{ \ for all } h \in \mathcal{B}, \\
\mathsf{E}_{\theta} [S_{bp} b (X) h (X) + m_{2} (O, h)] = 0 \text{ \ for all } h \in \mathcal{P}.
\end{array} \right.
\end{align}
And $\psi (\theta) = \mathsf{E}_{\theta} [m_{1} (O, b) + S_{0}] = \mathsf{E}_{\theta} [m_{2} (O, p) + S_{0}] = \mathsf{E}_{\theta} [- S_{bp} b (X) p (X) + S_{0}]$.
\item Further suppose the maps $h \mapsto \mathsf{E}_{\theta} [m_{1} (O, h)]$ and $h \mapsto \mathsf{E}_{\theta} [m_{2} (O, h)]$ for $h \in L_{2} (\mathsf{P}_{F})$ are continuous and linear with Riesz representers $\mathcal{R}_{1}(X) \equiv \mathcal{R}_{1}(X; \theta)$ and $\mathcal{R}_{2}(X) \equiv \mathcal{R}_{2}(X; \theta)$\footnote{\label{fn:riesz}The Riesz representer $\mathcal{R} (X)$ of a continuous linear functional $h \mapsto \mathsf{E} [m (O, h)]$ for $h \in L_{2} (\mathsf{P}_{F})$ is, by definition, the function of $X$ satisfying $\mathsf{E} [m (O, h)] = \mathsf{E} \left[ \mathcal{R} (X) h (X) \right]$; see \citet{chen2007large, chen2015sieve, chernozhukov2022debiased, chernozhukov2022automatic, rotnitzky2021characterization}.} respectively. Then $b (X) = - \mathcal{R}_{2} (X) / \lambda (X)$ and $p (X) = - \mathcal{R}_{1} (X) / \lambda (X)$; moreover, $b$ and $p$ are the (global) minimizers of the following minimization problems
\begin{equation}
\label{loss}
\begin{split}
b (\cdot) & = \mathsf{arg \; min}_{h \in L_{2} (\mathsf{P}_{F})} \mathsf{E}_{\theta} \left[ S_{bp} \frac{h (X)^{2}}{2} + m_{2} (O, h) \right], \\
p (\cdot) & = \mathsf{arg \; min}_{h \in L_{2} (\mathsf{P}_{F})} \mathsf{E}_{\theta} \left[ S_{bp} \frac{h (X)^{2}}{2} + m_{1} (O, h) \right].
\end{split}
\end{equation}
\end{itemize}
\end{proposition}

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 \eqref{eq:mean_zero} can be traced back to \citet{newey1994asymptotic, newey1997convergence, newey1998undersmoothing, newey2004twicing, chernozhukov2022automatic}.

\begin{example}[Some examples of Riesz representers for DR functionals]
\label{eg:riesz}\leavevmode

\begin{enumerate}
\item $\psi (\theta) = - \mathsf{E}_{\theta} [b (X)] = - \mathsf{E}_{\theta} [Y (a = 1)]$. Then $\mathsf{IF}_{1, \psi} (\theta) = \mathcal{H} (b, p) - \psi (\theta)$ with
\begin{equation*}
\mathcal{H} (b, p) = - \left( - A b (X) p (X) + b (X) + A Y p (X) \right).
\end{equation*}
Here $S_{bp} = A$ and $\lambda (X) = p (X)^{-1}$, $m_{1} (O, b) = - b(X)$, $m_{2} (O, p) = - A Y p (X)$ and $S_{0} = 0$. Then $\mathcal{R}_{1} (X) = -1$ and $\mathcal{R}_{2} (X) = - b (X) / p (X)$ since
\begin{equation*}
\left\{ \begin{array}{l}
\mathsf{E}_{\theta} \left[ m_{1} (O, h) \right] = \mathsf{E}_{\theta} [- h (X)], \\
\mathsf{E}_{\theta} \left[ m_{2} (O, h) \right] = \mathsf{E}_{\theta} [- A Y h (X)] = \mathsf{E}_{\theta} \left[ - \left\{ p (X) \right\}^{-1} b (X) h (X) \right].
\end{array} \right.
\end{equation*}

\item $\psi (\theta) = \mathsf{E}_{\theta} \left[ (Y - b (X)) (A - p (X)) \right]$ -- expected conditional covariance between $A$ and $Y$ given $X$ with $b (X) = \mathsf{E}_{\theta} [Y | X], p (X) = \mathsf{E}_{\theta} [A | X]$: $\mathsf{IF}_{1, \psi} (\theta) = \mathcal{H} (b, p) - \psi (\theta)$ with $\mathcal{H} (b, p) = b (X) p (X) - A b (X) - Y p (X) + A Y$, $S_{bp} = 1$ and $\lambda (X) = 1$, $m_{1} (O, b) = - A b (X)$, $m_{2} (O, p) = - Y p (X)$, and $S_{0} = A Y$. Then $\mathcal{R}_{1} (X) = - p (X)$ and $\mathcal{R}_{2} (X) = - b(X)$ since
\begin{equation*}
\left\{ \begin{array}{l}
\mathsf{E}_{\theta} \left[ m_{1} (O, h) \right] = \mathsf{E}_{\theta} [- A h(X)] = \mathsf{E}_{\theta} [- p(X) h(X)], \\
\mathsf{E}_{\theta} \left[ m_{2} (O, h) \right] = \mathsf{E}_{\theta} [- Y h(X)] = \mathsf{E}_{\theta} [- b(X) h(X)].
\end{array} \right.
\end{equation*}
\end{enumerate}
\end{example}

\allowdisplaybreaks
We now briefly comment on the three parts of Proposition \ref{thm:if1} as they are all important for future development of the paper. Eq. \eqref{eq:if1} 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. \eqref{eq:mean_zero} and \eqref{loss} 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:
\begin{lemma}
\label{lem:loss}
The minimization problems given in \eqref{loss}, with the function class $L_{2} (\mathsf{P}_{F})$ replaced by some $\mathcal{F} \subset L_{2} (\mathsf{P}_{F})$
\begin{equation}
\label{loss_f}
\begin{split}
\widetilde{b} (\cdot) = \mathsf{arg \; min}_{h \in \mathcal{F}} \mathsf{E}_{\theta} \left[ S_{bp} \frac{h (X)^{2}}{2} + m_{2} (O, h) \right], \widetilde{p} (\cdot) = \mathsf{arg \; min}_{h \in \mathcal{F}} \mathsf{E}_{\theta} \left[ S_{bp} \frac{h (X)^{2}}{2} + m_{1} (O, h) \right]
\end{split}
\end{equation}
are equivalent to:
\begin{equation}
\label{equi_loss_f}
\begin{split}
\widetilde{b} (\cdot) = \mathsf{arg \; min}_{h \in \mathcal{F}} \mathsf{E}_{\theta} \left[ \lambda (X) (b (X) - h (X))^{2} \right], \widetilde{p} (\cdot) = \mathsf{arg \; min}_{h \in \mathcal{F}} \mathsf{E}_{\theta} \left[ \lambda (X) (p (X) - h (X))^{2} \right].
\end{split}
\end{equation}
\end{lemma}
The proof of Lemma \ref{lem:loss} can be found in Appendix \ref{app:loss}. Based on Lemma \ref{lem:loss}, 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:
\begin{equation}
\label{loss_data}
\begin{split}
\widehat{b} = \underset{h \in \mathcal{F}}{\mathsf{arg \; min}} \mathsf{P}_{n_{\mathsf{tr}}} \left[ S_{bp} \frac{h^{2}}{2} + m_{2} (O, h) \right], \widehat{p} = \underset{h \in \mathcal{F}}{\mathsf{arg \; min}} \mathsf{P}_{n_{\mathsf{tr}}} \left[ S_{bp} \frac{h^{2}}{2} + m_{1} (O, h) \right].
\end{split}
\end{equation}
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 \eqref{loss_data} to the expected $\lambda$-weighted $L_{2} (\mathsf{P}_{\theta})$-loss \eqref{loss_f} 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{def:dr}:
\begin{equation}
\label{holder}
\inf_{x} \frac{\lambda (x)}{w_{b} (x)} \mathsf{E}_{\theta} [w_{b} (X) \{\widehat{b} (X) - b(X)\}^{2}] \leqslant \mathsf{E}_{\theta} [\lambda (X) \{\widehat{b} (X) - b(X)\}^{2}] \leqslant \sup_{x} \frac{\lambda (x)}{w_{b} (x)} \mathsf{E}_{\theta} [w_{b} (X) \{\widehat{b} (X) - b(X)\}^{2}]
\end{equation}
and similarly for $\mathsf{E}_{\theta} [\lambda (X) \{\widehat{p} (X) - p (X)\}^{2}]$. In Appendix \ref{app:loss}, 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
\begin{equation}
\label{cs-general}
\mathsf{CSBias}_{\theta} (\widehat{\psi}_{1}) \equiv \left\{ \mathsf{E}_{\theta} [\lambda (X) (\widehat{b} (X) - b (X))^{2}] \mathsf{E}_{\theta} [\lambda (X) (\widehat{p} (X) - p (X))^{2}] \right\}^{1 / 2}
\end{equation}
which is an upper bound of (the absolute value of) $\mathsf{Bias}_{\theta} (\widehat{\psi}_{1})$ to be defined in \eqref{bias}, by CS inequality. Hence we have
\begin{equation}
\label{equiv}
\mathsf{CSBias}_{\theta} (\widehat{\psi}_{1}) \asymp \mathsf{CSBias}_{\theta}^{w} (\widehat{\psi}_{1}) \equiv \{\mathsf{E}_{\theta} [w_{b} (X) \{\widehat{b} (X) - b(X)\}^{2}] \mathsf{E}_{\theta} [w_{p} (X) \{\widehat{p} (X) - p (X)\}^{2}]\}^{1 / 2},
\end{equation}
that is, they are equivalent up to a multiplicative positive constant.

\subsection{Standard DML estimators}\leavevmode
\label{sec:dml}

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{thm:if1}.

\begin{itemize}
\item[(i)] The $N$ study subjects are randomly split into two parts: an estimation sample of size $n$ and a training (nuisance) sample of size $n_{tr} = N - n$ with $n / N \approx 1 / 2$. Without loss of generality we shall assume that $i = 1, \ldots, n$ corresponds to the estimation sample.

\item[(ii)] Nuisance estimators $\widehat{b}$ and $\widehat{p}$ are separately constructed from the entire training sample data using machine learning, such as deep neural networks \citep{farrell2021deep, chen2020causal} and define $\widehat{\theta} \coloneqq (\widehat{b}, \widehat{p}, \theta \setminus \{b, p\})$.

\item[(iii)] Denote $\mathbb{IF}_{1} \equiv \mathbb{IF}_{1} (\theta) \coloneqq \frac{1}{n} \sum_{i = 1}^{n} \mathsf{IF}_{1, \psi} (\theta)$ and $\widehat{\mathbb{IF}}_{1} \equiv \mathbb{IF}_{1} (\widehat{\theta}) \coloneqq \frac{1}{n} \sum_{i = 1}^{n} \mathsf{IF}_{1, \psi} (\widehat{\theta})$. Then
\begin{equation*}
\begin{split}
\widehat{\psi}_{1} & = \widehat{\mathbb{IF}}_{1} + \psi (\widehat{\theta}) = \frac{1}{n} \sum_{i = 1}^{n} \mathcal{H}_{i} (\widehat{b}, \widehat{p}) = \frac{1}{n} \sum_{i = 1}^{n} \left( S_{bp, i} \widehat{b} (X_{i}) \widehat{p} (X_{i}) + m_{1} (O_{i}, \widehat{b}) + m_{2} (O_{i}, \widehat{p})+ S_{0, i} \right)
\end{split}
\end{equation*}
from $n$ subjects in the estimation sample and $\widehat{\psi}_{\mathsf{cf}, 1} = \dfrac{1}{2} \left( \widehat{\psi}_{1} + \overline{\widehat{\psi}}_{1} \right)$ where $\overline{\widehat{\psi}}_{1}$ is $\widehat{\psi}_{1}$ but with the training and estimation samples reversed.
\end{itemize}
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 \citet{chernozhukov2018double}, \citet{smucler2019unifying} or \citet[Theorem 3]{farrell2021deep}.
\begin{proposition}\label{thm:drml}
Under the conditions of Proposition \ref{thm:if1}, conditional on the training sample, if a) the bias of $\widehat{\psi}_{1}$
\begin{equation}
\label{bias}
\mathsf{Bias}_{\theta} (\widehat{\psi}_{1}) \coloneqq \mathsf{E}_{\theta} [S_{bp} (b (X) - \widehat{b} (X)) (p (X) - \widehat{p} (X))] \equiv \mathsf{E}_{\theta} [\lambda (X) (b (X) - \widehat{b} (X)) (p (X) - \widehat{p} (X))]
\end{equation}
is $o (n^{- 1 / 2})$ and b) $\widehat{b} (x)$ and $\widehat{p} (x)$ converge to $b (x)$ and $p (x)$ in $L_2 (\mathsf{P}_{\theta})$, then:
\begin{enumerate}
\item $\widehat{\psi}_{1} - \psi (\theta) = n^{-1} \sum_{i = 1}^n \mathsf{IF}_{1, \psi, i} (\theta) + o (n^{-1/2})$ and $\widehat{\psi}_{\mathsf{cf}, 1} - \psi (\theta) = N^{-1} \sum_{i = 1}^N \mathsf{IF}_{1, \psi, i} (\theta) + o_{p} (N^{-1/2})$. Further $n^{1/2} (\widehat{\psi}_{1} - \psi(\theta))$ converges conditionally and unconditionally to a normal distribution with mean zero; $\widehat{\psi}_{\mathsf{cf}, 1}$ is a regular, asymptotically linear estimator; i.e., $N^{1/2} (\widehat{\psi}_{\mathsf{cf}, 1} - \psi (\theta))$ converges unconditionally to a normal distribution with mean zero and variance equal to the semiparametric variance bound $\mathsf{var}_{\theta} \left[ \mathsf{IF}_{1, \psi} (\theta) \right]$.

\item The nominal $(1 - \alpha)$ Wald CIs
\begin{equation}
\mathsf{CI}_{\alpha} (\widehat{\psi}_{1}) \coloneqq \widehat{\psi}_{1} \pm z_{\alpha / 2} \widehat{\mathsf{s.e.}} [\widehat{\psi}_{1}], \ \mathsf{CI}_{\alpha} (\widehat{\psi}_{\mathsf{cf}, 1}) \coloneqq \widehat{\psi}_{\mathsf{cf}, 1} \pm z_{\alpha / 2} \widehat{\mathsf{s.e.}} [\widehat{\psi}_{\mathsf{cf}, 1}] \label{eq:ci}
\end{equation}
are asymptotically valid nominal $(1 - \alpha)$ CI for $\psi (\theta)$. Here $\widehat{\mathsf{s.e.}} [\widehat{\psi}_{1}] = \{\widehat{\mathsf{var}} [\widehat{\psi}_{1}]\}^{1 / 2}$ with
\begin{equation}
\label{var:dml}
\widehat{\mathsf{var}} [\widehat{\psi}_{1}] = \frac{1}{n^2} \sum_{i = 1}^{n} \left( \mathsf{IF}_{1, \psi, i} (\widehat{\theta}) - \frac{1}{n} \sum_{i = 1}^{n} \mathsf{IF}_{1, \psi, i} (\widehat{\theta}) \right)^{2}
\end{equation}
and $\widehat{\mathsf{s.e.}} [\widehat{\psi}_{\mathsf{cf}, 1}] = \frac{1}{2} \{\widehat{\mathsf{var}} [\widehat{\psi}_{1}] + \widehat{\mathsf{var}} [\overline{\widehat{\psi}}_{1}]\}^{1 / 2}$. The interval $\mathsf{CI}_{\alpha} (\widehat{\psi}_{1})$, unlike $\mathsf{CI}_{\alpha} (\widehat{\psi}_{\mathsf{cf}, 1})$, is also an asymptotically valid $(1 - alpha)$ CI conditional on the training sample.
\end{enumerate}
\end{proposition}

\allowdisplaybreaks
\subsection{A prelude to our approach}\leavevmode
\label{sec:prelude}

As briefly described in Section \ref{sec:introduction}, 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}$:
\begin{equation}
\label{decomposition}
\mathsf{Bias}_{\theta} (\widehat{\psi}_{1}) = \mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1}) + \mathsf{TB}_{\theta, k} (\widehat{\psi}_{1})
\end{equation}
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 \citet{robins2008higher, liu2020nearly} (also see Appendix \ref{app:hoif}), and
\begin{equation}\label{bias_k}
\begin{split}
\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1}) & \coloneqq \mathsf{E}_{\theta} [\Pi_{\theta} [\lambda^{1 / 2} (\widehat{b} - b) | \lambda^{1 / 2} \bar{\mathsf{z}}_{k}] (X) \Pi_{\theta} [\lambda^{1 / 2} (\widehat{p} - p) | \lambda^{1 / 2} \bar{\mathsf{z}}_{k}] (X)] \\
& \equiv \mathsf{E}_{\theta} [\lambda (X) (\widehat{b} (X) - b (X)) \bar{\mathsf{z}}_{k} (X)]^{\top} \Sigma_{k}^{-1} \mathsf{E}_{\theta} [\bar{\mathsf{z}}_{k} (X) \lambda (X) (\widehat{p} (X) - p (X))].
\end{split}
\end{equation}
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{cond:sw} in Section \ref{sec:main}) 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{sec:introduction}, $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 \eqref{cs-heuristic}, 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{sec:main}, 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 \citep{robins2008higher, robins2017minimax, liu2017semiparametric}. Due to space limitation, we refer the interested readers to \citet{robins2008higher} or \citet{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{rem:finite}.


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
\begin{equation}\label{if22}
\begin{split}
\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1}) \coloneqq \frac{1}{n (n - 1)} \sum_{1 \leqslant i_{1} \neq i_{2} \leqslant n} \left[ \mathcal{E}_{\widehat{b}, m_{2}} (\bar{\mathsf{z}}_{k}) (O) \right]_{i_{1}}^{\top} \Sigma_{k}^{-1} \left[ \mathcal{E}_{\widehat{p}, m_{1}} (\bar{\mathsf{z}}_{k}) (O) \right]_{i_{2}}, \text{ where} \\
\mathcal{E}_{\widehat{b}, m_{2}} (\bar{\mathsf{z}}_{k}) (O) \coloneqq S_{bp} \widehat{b} (X) \bar{\mathsf{z}}_{k} (X) + m_{2} (O, \bar{\mathsf{z}}_{k}), \mathcal{E}_{\widehat{p}, m_{1}} (\bar{\mathsf{z}}_{k}) (O) \coloneqq S_{bp} \widehat{p} (X) \bar{\mathsf{z}}_{k} (X) + m_{1} (O, \bar{\mathsf{z}}_{k})
\end{split}
\end{equation}
with standard error of order $O \left( \frac{\sqrt{k}}{n} \vee \frac{1}{\sqrt{n}} \right)$ (see Theorem \ref{thm:hoif_stats} 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 \citet{robins2008higher}.}; see Appendix \ref{app:hoif} 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 \eqref{eq:mean_zero} in Proposition \ref{thm:if1}; the proof of the bound on the standard error of $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$ is deferred to Appendix \ref{app:oracle}, and can also be found in \citet{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.

\begin{remark}[On the choice of $k$]\leavevmode
\label{rem:k}
We shall always choose $k$ to be much less than the sample size $n = N / 2$ of the estimation sample, for the standard error of $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$ to be smaller than or equal to the order $n^{-1/2}$ of the standard error of $\widehat{\psi}_{1}$, thereby creating the possibility of detecting, for any given $\delta > 0$, that $|\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})| / \mathsf{s.e.}_{\theta} (\widehat{\psi}_{1}) \geqslant \delta$, when the sample size $n$ is sufficiently large.

Moreover, we also need $k \ll n_{\mathsf{tr}} = n$, the sample size of the training sample, to ensure $\widehat{\Sigma}_{k}^{-1}$ to be a operator-norm consistent estimator of $\Sigma_{k}^{-1}$, without imposing, possibly incorrect, additional smoothness assumption on the distribution of $(S_{bp}, X)$ or sparsity assumptions on $\Sigma_{k}$ or $\Sigma_{k}^{-1}$ required for accurate estimation when $k \geqslant n$.

Finally, we take $k \gg 1$ to ensure the asymptotic normality of the $U$-statistic $\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$ (after appropriate scaling) \citep{bhattacharya1992class} to determine the rejection region of our falsification test to be studied in the next section.

One might view the choice of $k \ll n$ a disadvantage. On the surface, it seems to exclude high-dimensional statistical models. However, this is not the case: (1) The nuisance functions are allowed to be estimated by any procedure, be it high-dimensional or not, from the training sample; (2) $k \ll n$ ensures that our testing procedure can be a valid test without any complexity-reducing assumptions on the nuisance functions, including the distribution of $X$, that is required to evaluate $\Sigma_{k}^{-1}$. We cannot let $k \asymp n$ because, as we will see, our test needs the bias due to estimating $\Sigma_{k}$ by $\widehat{\Sigma}_{k}$ to be sufficiently small and recent developments along this direction \citep{cattaneo2018inference, kline2020leave, yadlowsky2022explaining, jiang2022new, jochmans2022heteroscedasticity} cannot be used because of the more stringent assumptions on $\bar{\mathsf{z}}_{k}$, such as sub-Gaussian tail conditions. The jackknife or bootstrap bias correction methods developed by \citet{cattaneo2019two, cattaneo2018kernel} could potentially be applied to our problem but since they only allow $k = O (\sqrt{n})$, the chance of rejecting $\mathsf{H}_{0, \mathsf{CS}} (\delta)$ could be diminished. It will be an interesting problem to study if their approach can be extended to allow $\sqrt{n} \ll k \ll n$, so as to be applied to our problem. We list this as an open problem in Section \ref{sec:conclusion}.
\end{remark}






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

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}}$:
\begin{align*}
\mathsf{H}_{0, \mathsf{CS}} (\delta): \frac{\mathsf{CSBias}_{\theta} (\widehat{\psi}_{1})}{\mathsf{s.e.}_{\theta} (\widehat{\psi}_{1})} \leqslant \delta.
\end{align*}

\subsection{What constitutes a valid test of $\mathsf{H}_{0, \mathsf{CS}} (\delta)$?}\leavevmode
\label{sec:what}

Before proceeding to the main result, we note that the infeasible falsification test of $\mathsf{H}_{0, \mathsf{CS}} (\delta)$ proposed in \citet{liu2020nearly} can be generalized to the class of DR functionals as follows\footnote{We also provide relevant results in Appendix \ref{app:oracle_test}.}:
\begin{equation}
\widehat{\chi}_{2, k} (\Sigma_{k}^{-1}; \varsigma_{k}, \delta) = \mathbbm{1} \left\{ \frac{|\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})|}{\widehat{\mathsf{s.e.}} [\widehat{\psi}_{1}]} - \varsigma_{k} \frac{\widehat{\mathsf{s.e.}} [\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})]}{\widehat{\mathsf{s.e.}} [\widehat{\psi}_{1}]} \geqslant \delta \right\}
\end{equation}
where $\widehat{\mathsf{s.e.}} [\widehat{\psi}_{1}]$ (see \eqref{var:dml}) and $\widehat{\mathsf{s.e.}} [\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})]$ (see \eqref{eq:if22_bvar} in Appendix \ref{app:bootstrap}) 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.

\begin{proposition}
\label{prop:master}
Given any estimator $\widehat{\mathsf{Bias}}_{\theta, k} (\widehat{\psi}_{1})$ of $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$ that, conditioning on the training sample data, satisfies: \\
(1) Under $\mathsf{H}_{0, \mathsf{CS}} (\delta)$: $\mathsf{E}_{\theta} [\widehat{\mathsf{Bias}}_{\theta, k} (\widehat{\psi}_{1}) - \mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})] = o \left( \frac{\sqrt{k}}{n} \right)$; \\
(2) $\mathsf{s.e.}_{\theta} [\widehat{\mathsf{Bias}}_{\theta, k} (\widehat{\psi}_{1})] \asymp \mathsf{s.e.}_{\theta} [\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})]$; \\
(3) Under $\mathsf{H}_{0, \mathsf{CS}} (\delta)$: $\frac{\widehat{\mathsf{Bias}}_{\theta, k} (\widehat{\psi}_{1}) - \mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})}{\mathsf{s.e.}_{\theta} [\widehat{\mathsf{Bias}}_{\theta, k} (\widehat{\psi}_{1})]} \overset{d}{\rightarrow} N (0, 1)$,
then (again conditioning on the training sample data) the test
\begin{equation}
\widehat{\chi}_{k} (\varsigma_{k}, \delta) = \mathbbm{1} \left\{ \frac{|\widehat{\mathsf{Bias}}_{\theta, k} (\widehat{\psi}_{1})|}{\widehat{\mathsf{s.e.}} [\widehat{\psi}_{1}]} - \varsigma_{k} \frac{\widehat{\mathsf{s.e.}} [\widehat{\mathsf{Bias}}_{\theta, k} (\widehat{\psi}_{1})]}{\widehat{\mathsf{s.e.}} [\widehat{\psi}_{1}]} \geqslant \delta \right\},
\end{equation}
where $\widehat{\mathsf{s.e.}} [\widehat{\mathsf{Bias}}_{\theta, k} (\widehat{\psi}_{1})]$ is a consistent estimator of $\mathsf{s.e.}_{\theta} [\widehat{\mathsf{Bias}}_{\theta, k} (\widehat{\psi}_{1})]$, rejects $\mathsf{H}_{0, \mathsf{CS}} (\delta)$ with probability less than $\alpha^{\dag}$ as $n \rightarrow \infty$ with $\varsigma_{k} = z_{\alpha^{\dag} / 2}$ under $\mathsf{H}_{0, \mathsf{CS}} (\delta)$. That is, we say that $\widehat{\chi}_{k} (z_{\alpha^{\dag} / 2}, \delta)$ is an asymptotically valid test of $\mathsf{H}_{0, \mathsf{CS}} (\delta)$.

The above statement also holds if one replaces $\mathsf{H}_{0, \mathsf{CS}} (\delta)$ by the $k$-projected null hypothesis $\mathsf{H}_{0, k} (\delta)$.
\end{proposition}

We did not give a formal proof as the argument for Proposition \ref{prop:master} 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{prop:master} will be most relevant to Section \ref{sec:nonstandard}. We also refer readers to Appendix \ref{app:sigma_test} for more details.

\subsection{Bias-reduced estimator of $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$}\leavevmode
\label{sec:soif}

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
\begin{equation}
\widehat{\chi}_{2, k} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta) = \mathbbm{1} \left\{ \frac{|\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})|}{\widehat{\mathsf{s.e.}} [\widehat{\psi}_{1}]} - z_{\alpha^{\dag} / 2} \frac{\widehat{\mathsf{s.e.}} [\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})]}{\widehat{\mathsf{s.e.}} [\widehat{\psi}_{1}]} \geqslant \delta \right\}.
\end{equation}
However, we {\it cannot} prove that the above test is asymptotically valid for $\mathsf{H}_{0, \mathsf{CS}} (\delta)$, because Theorem \ref{thm:soif_toif} 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
\begin{equation}
\label{eb2:upper_null}
\mathsf{EB}_{\theta, k, 2} \coloneqq \mathsf{E}_{\theta} [\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1}) - \mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})] = O \left( \frac{\sqrt{k}}{n} \sqrt{\mathsf{log} k} \right)
\end{equation}
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{prop:master}.

Fortunately, Theorem \ref{thm:soif_toif} will also show that the following third-order $U$-statistic estimator of $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$
\begin{equation}
\label{if2233}
\begin{split}
& \ \widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1}) \coloneqq \widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1}) + \widehat{\mathbb{IF}}_{33, k} (\widehat{\Sigma}_{k}^{-1}), \text{ where} \\
\widehat{\mathbb{IF}}_{33, k} (\widehat{\Sigma}_{k}^{-1}) \coloneqq & \ \frac{(n - 2)!}{n!} \sum_{1 \leqslant i_{1} \neq i_{2} \neq i_{3} \leqslant n} \left[ \mathcal{E}_{\widehat{b}, m_{2}} (\bar{\mathsf{z}}_{k}) (O) \right]_{i_{1}}^{\top} \widehat{\Sigma}_{k}^{-1} \left[ S_{bp} \bar{\mathsf{z}}_{k} (X) \bar{\mathsf{z}}_{k} (X)^{\top} - \widehat{\Sigma}_{k} \right]_{i_{3}} \widehat{\Sigma}_{k}^{-1} \left[ \mathcal{E}_{\widehat{p}, m_{1}} (\bar{\mathsf{z}}_{k}) (O) \right]_{i_{2}}
\end{split}
\end{equation}
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 \eqref{eb2:upper_null} to
\begin{equation}
\label{eb3:upper_null}
\mathsf{EB}_{\theta, k, 3} \coloneqq \mathsf{E}_{\theta} [\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1}) - \mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})] = O \left( \frac{\sqrt{k}}{n} \sqrt{\frac{k \mathsf{log} k}{n}} \right) = o \left( \frac{\sqrt{k}}{n} \right)
\end{equation}
as long as we choose $k$ such that $k \mathsf{log} k = o (n)$.

\begin{remark}
$\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$ can be viewed as a de-biased version of $\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$, and the bias is partially corrected by $\widehat{\mathbb{IF}}_{33, k} (\widehat{\Sigma}_{k}^{-1})$. As will be shown in Appendix \ref{app:hoif}, $\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\Sigma_{k}^{-1})$ is the third-order influence function of $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$ and $\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$ is its estimated version.
\end{remark}

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$:
\begin{align*}
\mathbb{L}_{\theta, 2, \widehat{b}} & \coloneqq \left\{ \mathsf{E}_{\theta} [\lambda (X) (\widehat{b} (X) - b (X))^{2}] \right\}^{1 / 2}, \mathbb{L}_{\theta, q, \widehat{b}, k} \coloneqq \left\{ \mathsf{E}_{\theta} [\{\Pi_{\theta} [\lambda^{1 / 2} (\widehat{b} - b) \vert \lambda^{1 / 2} \bar{\mathsf{z}}_{k}] (X)\}^{q}] \right\}^{1 / q}, \\
\mathbb{L}_{\theta, 2, \widehat{p}} & \coloneqq \left\{ \mathsf{E}_{\theta} [\lambda (X) (\widehat{p} (X) - p (X))^{2}] \right\}^{1 / 2}, \mathbb{L}_{\theta, q, \widehat{p}, k} \coloneqq \left\{ \mathsf{E}_{\theta} [\{\Pi_{\theta} [\lambda^{1 / 2} (\widehat{p} - p) \vert \lambda^{1 / 2} \bar{\mathsf{z}}_{k}] (X)\}^{q}] \right\}^{1 / q}, \\
\mathbb{L}_{\theta, \infty, \widehat{b}, k} & \coloneqq \Vert \Pi_{\theta} [\lambda^{1 / 2} (\widehat{b} - b) \vert \lambda^{1 / 2} \bar{\mathsf{z}}_{k}] \Vert_{\infty}, \mathbb{L}_{\theta, \infty, \widehat{p}, k} \coloneqq \Vert \Pi_{\theta} [\lambda^{1 / 2} (\widehat{p} - p) \vert \lambda^{1 / 2} \bar{\mathsf{z}}_{k}] \Vert_{\infty}, \mathbb{L}_{\theta, 2, \widehat{\Sigma}, k} \coloneqq \Vert \widehat{\Sigma}_{k} - \Sigma_{k} \Vert,
\end{align*}
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 \eqref{csbias} to the whole class of DR functionals.

We also need to further impose the following weak regularity conditions (Condition \ref{cond:sw}).

\begin{customthm}{W}
\leavevmode\label{cond:sw}

\begin{enumerate}
\item All the eigenvalues of $\Sigma_{k}$ are bounded away from 0 and $\infty$.

\item The observed data $O = (W, X)$, the true nuisance functions $b(X)$ and $p(X)$, the estimated nuisance functions $\widehat{b}(X)$ and $\widehat{p}(X)$ are bounded with $\mathsf{P}_{\theta}$-probability 1, and $\lambda (X)$ are bounded away from 0 and $\infty$ with $\mathsf{P}_{\theta}$-probability 1.

\item $\Vert \bar{\mathsf{z}}_{k}^{\top} \bar{\mathsf{z}}_{k} \Vert_{\infty} \leqslant B k$ for some constant $B > 0$.
\end{enumerate}
\end{customthm}

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

\begin{customthm}{S}
\leavevmode\label{cond:w}

\begin{enumerate} [label=(\arabic*')] \setcounter{enumi}{2}

\item $\Vert \bar{\mathsf{z}}_{k}^{\top} \bar{\mathsf{z}}_{k} \Vert_{\infty} \leqslant B k$ for some constant $B > 0$; in addition, $\mathbb{L}_{\theta, \infty, \widehat{b}, k} < \infty$ and $\mathbb{L}_{\theta, \infty, \widehat{p}, k} < \infty$.
\end{enumerate}
\end{customthm}

\begin{remark}
\label{rem:w}
Condition \ref{cond:w}(3') is stronger than Condition \ref{cond:sw}(3) and it holds for Cohen-Vial-Daubechies wavelets, local polynomial partition and B-spline series \citep{newey1997convergence, belloni2015some}. But it does not hold in general for Fourier or Legendre polynomial series of $X$ when $X$ is compactly supported \citep{belloni2015some}. We also refer interested readers to \citet[Section 6]{kennedy2020discussion} and \citet[Section 2.1]{liu2020rejoinder} for further motivation on why we decide to relax Condition \ref{cond:w} by Condition \ref{cond:sw}. However, as we will see, when $\Sigma_{k}$ is unknown and need to be estimated from the training sample, violation of Condition \ref{cond:w}(3') but not Condition \ref{cond:sw}(3), will incur a loss in the power of the test (the magnitude depending on the particular series used; see Appendix \ref{app:var}), but the validity of the test (i.e. level) is still guaranteed. Nevertheless, such loss in power happens (or equivalently, Condition \ref{cond:sw} holds but Condition \ref{cond:w} fails to hold), only if the $L_{\infty}$-norms of the projections $\Pi_{\theta} [\lambda^{1 / 2} (\widehat{b} - b) \vert \lambda^{1 / 2} \bar{\mathsf{z}}_{k}]$ and $\Pi_{\theta} [\lambda^{1 / 2} (\widehat{p} - p) \vert \lambda^{1 / 2} \bar{\mathsf{z}}_{k}]$ are not bounded even though those of $\lambda^{1 / 2} (\widehat{b} - b)$ and $\lambda^{1 / 2} (\widehat{p} - p)$ are bounded.
\end{remark}

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{thm:soif_toif} below:

\begin{theorem}\label{thm:soif_toif}
Under the conditions of Proposition \ref{thm:if1} and Condition \ref{cond:sw}, with $k, n \rightarrow \infty$ but $k \mathsf{log} k = o (n)$, conditioning on the training sample, the following hold on the event that $\widehat{\Sigma}_{k}$ is invertible:

\begin{enumerate}
\item The bias and variance of $\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$ as an estimator of $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$ are of the following order:
\begin{align*}
\left\vert \mathsf{EB}_{\theta, 2, k} \right\vert & \lesssim \mathbb{L}_{\theta, 2, \widehat{b}, k} \mathbb{L}_{\theta, 2, \widehat{p}, k} \mathbb{L}_{\theta, 2, \widehat{\Sigma}, k} \lesssim \mathbb{L}_{\theta, 2, \widehat{b}, k} \mathbb{L}_{\theta, 2, \widehat{p}, k} \sqrt{\frac{k \mathsf{log} (k)}{n}}, \\
\mathsf{var}_{\theta} \left[\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1}) \right] & \lesssim \frac{1}{n} \left\{ \frac{k}{n} + \mathbb{L}_{\theta, 2, \widehat{b}, k}^{2} + \mathbb{L}_{\theta, 2, \widehat{p}, k}^{2} \right\}
\end{align*}
where $\mathsf{EB}_{\theta, 2, k}$ is defined in equation \eqref{eb2:upper_null}.

\item The bias and variance of $\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1}) \coloneqq \widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1}) + \widehat{\mathbb{IF}}_{33, k} (\widehat{\Sigma}_{k}^{-1})$ are of the following order:
\begin{align*}
\left\vert \mathsf{EB}_{\theta, 3, k} \right\vert & \lesssim \mathbb{L}_{\theta, 2, \widehat{b}, k} \mathbb{L}_{\theta, 2, \widehat{p}, k} \mathbb{L}_{\theta, 2, \widehat{\Sigma}, k}^{2} \lesssim \mathbb{L}_{\theta, 2, \widehat{b}, k} \mathbb{L}_{\theta, 2, \widehat{p}, k} \frac{k \mathsf{log} (k)}{n}, \\
\mathsf{var}_{\theta} \left[ \widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1}) \right] & \lesssim \frac{1}{n} \left\{ \frac{k}{n} + \mathbb{L}_{\theta, 2, \widehat{b}, k}^{2} + \mathbb{L}_{\theta, 2, \widehat{p}, k}^{2} + k \mathbb{L}_{\theta, 2, \widehat{b}, k}^{2} \mathbb{L}_{\theta, 2, \widehat{p}, k}^{2} \right\}
\end{align*}
where $\mathsf{EB}_{\theta, 3, k}$ is defined in equation \eqref{eb2:upper_null}.

If, however, Condition \ref{cond:sw} is strengthened to Condition \ref{cond:w}
\begin{equation*}
\mathsf{var}_{\theta} \left[ \widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1}) \right] \lesssim \frac{1}{n} \left\{ \frac{k}{n} + \mathbb{L}_{\theta, 2, \widehat{b}, k}^{2} + \mathbb{L}_{\theta, 2, \widehat{p}, k}^{2} \right\}.
\end{equation*}

\item Under $\mathsf{H}_{0, \mathsf{CS}} (\delta)$:
\begin{equation*}
\frac{\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1}) - \mathsf{E}_{\theta} \left[ \widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1}) \right]}{\mathsf{s.e.}_{\theta} [\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})]} \text{ and } \frac{\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1}) - \mathsf{E}_{\theta}\left[\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})\right]}{\mathsf{s.e.}_{\theta} [\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})]} \overset{d}{\rightarrow} N (0, 1).
\end{equation*}
When we estimate $\mathsf{s.e.}_{\theta} [\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$ and $\mathsf{s.e.}_{\theta} [\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$ respectively by their consistent bootstrap estimators $\widehat{\mathsf{s.e.}} [\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$ and $\widehat{\mathsf{s.e.}} [\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$ defined in \eqref{eq:if22_bvar} and \eqref{bvar:if2233}, with $\Sigma_{k}$ replaced by $\widehat{\Sigma}_{k}$, in Appendix \ref{app:bootstrap}, we also have
\begin{equation*}
\frac{\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1}) - \mathsf{E}_{\theta} \left[ \widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1}) \right]}{\widehat{\mathsf{s.e.}} [\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})]} \text{ and } \frac{\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1}) - \mathsf{E}_{\theta}\left[\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})\right]}{\widehat{\mathsf{s.e.}} [\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})]} \overset{d}{\rightarrow} N (0, 1).
\end{equation*}
\end{enumerate}
\end{theorem}

\begin{remark}\leavevmode
\begin{itemize}
\item The bias and variance bounds under Condition \ref{cond:w} have been proved in \citet{liu2017semiparametric}. The change in the proof under Condition \ref{cond:sw} is given in Lemma \ref{lem:hoif_var} in Appendix \ref{app:var}. The normal approximation under $\mathsf{H}_{0, \mathsf{CS}} (\delta)$ follows from Slutsky theorem and asymptotic normality of $U$-statistics with diverging (with $n$) kernels given in Theorem 1 of \citet{bhattacharya1992class}.
\item The generic bias bounds in Theorem \ref{thm:soif_toif} imply those stated in \eqref{eb2:upper_null} and \eqref{eb3:upper_null} under $\mathsf{H}_{0, \mathsf{CS}} (\delta)$, as $$\mathbb{L}_{\theta, 2, \widehat{b}, k} \mathbb{L}_{\theta, 2, \widehat{p}, k} \leqslant \mathbb{L}_{\theta, 2, \widehat{b}} \mathbb{L}_{\theta, 2, \widehat{p}} = \mathsf{CSBias}_{\theta} (\widehat{\psi}_{1}) \lesssim n^{- 1 / 2}$$ when $\mathsf{H}_{0, \mathsf{CS}} (\delta)$ is true.
\end{itemize}
\end{remark}

\subsection{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)$:
\begin{equation}
\label{the_test}
\widehat{\chi}_{3, k} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta) \coloneqq \mathbbm{1} \left\{ \frac{|\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})|}{\widehat{\mathsf{s.e.}} [\widehat{\psi}_{1}]} - z_{\alpha^{\dag} / 2} \frac{\widehat{\mathsf{s.e.}} [\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})]}{\widehat{\mathsf{s.e.}} [\widehat{\psi}_{1}]} \geqslant \delta \right\},
\end{equation}
and (2) the following theorem showing its asymptotically validity:
\begin{theorem}
\label{thm:cs}
Assume all the conditions in Theorem \ref{thm:soif_toif}. Under $\mathsf{H}_{0, \mathsf{CS}} (\delta): \frac{\mathsf{CSBias}_{\theta} (\widehat{\psi}_{1})}{\mathsf{s.e.}_{\theta} (\widehat{\psi}_{1})} \leqslant \delta$,
\begin{equation}
\lim_{n \rightarrow \infty} \mathbb{P}_{\theta} \left( \widehat{\chi}_{3, k} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta) = 1 \right) \leqslant \alpha^{\dag}.
\end{equation}
That is, $\widehat{\chi}_{3, k} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta)$ is an asymptotically level-$\alpha^{\dag}$ test of $\mathsf{H}_{0, \mathsf{CS}} (\delta)$.
\end{theorem}

The proof of Theorem \ref{thm:cs} is a direct consequence of Theorem \ref{thm:soif_toif} and is deferred to Appendix \ref{app:sigma_test}. 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)$.

\begin{remark}
\label{rem:finite}
Because our results are conditioning on the training sample, the test $\widehat{\chi}_{3, k} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta)$ only relies on `asymptopia' to ensure that its rejection probability can be closely approximated by its Gaussian limiting distribution. If Berry-Esseen bound or tail inequalities with constants estimable from the data were available for $\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$, we could in principle eliminate our dependence on asymptopia. At the sample size used in our simulation (see Section \ref{sec:simulation}), the normal qqplots of $\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$ in Figure \ref{fig:qq_h0} suggest that normal approximation is reasonably close to the true distribution between $(-2, 2)$ under $\mathsf{H}_{0, \mathsf{CS}} (\delta)$.
\end{remark}

\begin{remark}
\label{rem:temper}
Now we return to the issue raised in footnote 5 on page 7 in the Introduction. Another natural choice for the $k$-projected null hypothesis is
\begin{align*}
\mathsf{H}_{0, k, \mathsf{CS}} (\delta): \mathsf{CSBias}_{\theta, k} (\widehat{\psi}_{1})^{2} \equiv \mathsf{E}_{\theta} \left[ \Pi_{\theta} [\lambda^{1 / 2} (\widehat{b} - b) | \lambda^{1 / 2} \bar{\mathsf{z}}_{k}] (X)^{2} \right] \mathsf{E}_{\theta} \left[ \Pi_{\theta} [\lambda^{1 / 2} (\widehat{p} - p) | \lambda^{1 / 2} \bar{\mathsf{z}}_{k}] (X)^{2} \right] \leqslant \delta^{2} \mathsf{var}_{\theta} (\widehat{\psi}_{1})
\end{align*}
because rejection of $\mathsf{H}_{0, k, \mathsf{CS}} (\delta)$ implies rejection of $\mathsf{H}_{0, \mathsf{CS}} (\delta)$. $\mathsf{CSBias}_{\theta, k} (\widehat{\psi}_{1})^{2}$ can be unbiasedly estimated by the following infeasible fourth-order $U$-statistic:
\begin{equation*}
\widehat{\mathbb{IF}}_{44, \mathsf{CS}, k} (\Sigma_{k}^{-1}) \coloneqq \frac{(n - 4)!}{n!} \sum_{1 \leqslant i_{1} \neq i_{2} \neq i_{3} \neq i_{4} \leqslant n} [\mathcal{E}_{\widehat{b}, m_{2}} (\bar{\mathsf{z}}_{k}) (O)]^{\top}_{i_{1}} \Sigma_{k}^{-1} [\mathcal{E}_{\widehat{b}, m_{2}} (\bar{\mathsf{z}}_{k}) (O)]_{i_{2}} [\mathcal{E}_{\widehat{p}, m_{1}} (\bar{\mathsf{z}}_{k}) (O)]^{\top}_{i_{3}} \Sigma_{k}^{-1} [\mathcal{E}_{\widehat{p}, m_{1}} (\bar{\mathsf{z}}_{k}) (O)]_{i_{4}}.
\end{equation*}
Using $\widehat{\mathbb{IF}}_{44, \mathsf{CS}, k} (\Sigma_{k}^{-1})$ as a starting point, one can indeed follow the previous development in Section \ref{sec:main} to construct a feasible test of $\mathsf{H}_{0, k, \mathsf{CS}} (\delta)$. We decide not to further pursue this alternative strategy because it is much more complex to analyze and we leave it to future investigation.
\end{remark}

\section{What if nonstandard DML estimators are used in practice?}
\label{sec:nonstandard}

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

Several recent works \citep{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 \text{H\"{o}lder}{} smoothness \citep{newey2018cross, mcgrath2022undersmoothing} or (approximate) sparsity \citep{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 \text{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 \citep{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{sec:prelude} that $\mathsf{Bias}_{\theta} (\widehat{\psi}_{1})$ can be decomposed as follows:
\begin{equation*}
\mathsf{Bias}_{\theta} (\widehat{\psi}_{1}) \equiv \mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1}) + \mathsf{TB}_{\theta, k} (\widehat{\psi}_{1}).
\end{equation*}
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.

\subsection{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{thm:cs} 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{thm:soif_toif} 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:
\begin{align*}
\mathsf{EB}_{3, \theta, k} = \mathsf{E}_{\theta} [\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1}) - \mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})] \lesssim \mathbb{L}_{\theta, 2, \widehat{b}, k} \mathbb{L}_{\theta, 2, \widehat{p}, k} \frac{k \mathsf{log} k}{n}.
\end{align*}
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{prop:master} 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})$:
\begin{equation*}
\begin{split}
& \widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1}) \coloneqq \sum_{j = 2}^{m} \widehat{\mathbb{IF}}_{jj, k} (\widehat{\Sigma}_{k}^{-1}), \text{ where} \\
\small\widehat{\mathbb{IF}}_{jj, k} (\widehat{\Sigma}_{k}^{-1}) & \coloneqq (-1)^{j} \frac{(n - j)!}{n!} \left[ \mathcal{E}_{\widehat{b}, m_{2}} (\bar{\mathsf{z}}_{k}) (O) \right]^{\top}_{i_{1}} \left\{ \prod_{s = 3}^{j} \widehat{\Sigma}_{k}^{-1} \left( \left[ S_{bp} \bar{\mathsf{z}}_{k} (X) \bar{\mathsf{z}}_{k} (X)^{\top} \right]_{i_{s}} - \widehat{\Sigma}_{k} \right) \right\} \widehat{\Sigma}_{k}^{-1} \left[ \mathcal{E}_{\widehat{p}, m_{1}} (\bar{\mathsf{z}}_{k}) (O) \right]_{i_{2}}.
\end{split}
\end{equation*}
$\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{app:hoif})]. 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{thm:soif_toif}. Again, the bias and variance bounds under Condition \ref{cond:w} are proved in \citet{liu2017semiparametric}. Under Condition \ref{cond:sw}, see the proof of Lemma \ref{lem:hoif_var} in Appendix \ref{app:var}. Recall from Section \ref{sec:main} that Condition \ref{cond:sw} is weaker than Condition \ref{cond:w}.

\begin{theorem}\label{thm:hoif_stats}
Under the conditions as in Theorem \ref{thm:soif_toif} and $k \mathsf{log} k \vee (k m^{2}) = o (n)$:
\begin{enumerate}[label = (\arabic*)]
\item The bias and variance of $\widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1})$ as an estimator of $\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1})$ are of the following order:
\begin{align*}
\left\vert \mathsf{EB}_{\theta, m, k} \right\vert & \coloneqq \mathsf{E}_{\theta} \left[ \widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1}) - \mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1}) \right] \\
& \lesssim \mathbb{L}_{\theta, 2, \widehat{b}, k} \mathbb{L}_{\theta, 2, \widehat{p}, k} \mathbb{L}_{\theta, 2, \widehat{\Sigma}, k}^{\frac{m - 1}{2}} \lesssim \mathbb{L}_{\theta, 2, \widehat{b}, k} \mathbb{L}_{\theta, 2, \widehat{p}, k} \left( \frac{k \mathsf{log} k}{n} \right)^{\frac{m - 1}{2}}, \\
\mathsf{var}_{\theta} \left[ \widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1}) \right] & \lesssim \frac{1}{n} \left\{ \frac{k}{n} + \mathbb{L}_{\theta, 2, \widehat{b}, k} + \mathbb{L}_{\theta, 2, \widehat{p}, k} + k \mathbb{L}_{\theta, 2, \widehat{b}, k}^{2} \mathbb{L}_{\theta, 2, \widehat{p}, k}^{2} \right\}.
\end{align*}
If, however, Condition \ref{cond:sw} is strengthened to Condition \ref{cond:w},
\begin{align*}
\mathsf{var}_{\theta} \left[ \widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1}) \right] \lesssim \frac{1}{n} \left\{ \frac{k}{n} + \mathbb{L}_{\theta, 2, \widehat{b}, k} + \mathbb{L}_{\theta, 2, \widehat{p}, k} \right\}.
\end{align*}
\item Under $\mathsf{H}_{0, k} (\delta)$ and Condition \ref{cond:w},
\begin{equation*}
\frac{\widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1}) - \mathsf{E}_{\theta}\left[\widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1})\right]}{\mathsf{s.e.}_{\theta} [\widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1})]} \overset{d}{\rightarrow} N (0, 1).
\end{equation*}
When we estimate $\mathsf{s.e.}_{\theta} [\widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1})$ by its consistent bootstrap estimator $\widehat{\mathsf{s.e.}} [\widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1})$ (e.g. by \eqref{bvar:if2233} in Appendix \ref{app:bootstrap} with $\Sigma_{k}$ replaced by $\widehat{\Sigma}_{k}$), we also have
\begin{equation*}
\frac{\widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1}) - \mathsf{E}_{\theta}\left[\widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1})\right]}{\widehat{\mathsf{s.e.}} [\widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1})]} \overset{d}{\rightarrow} N (0, 1).
\end{equation*}
\end{enumerate}
\end{theorem}

\begin{remark}
\label{rem:comp}
First, note that under $\mathsf{H}_{0, k} (\delta)$ instead of $\mathsf{H}_{0, \mathsf{CS}} (\delta)$ as in the previous section, we only have asymptotic normality of $\widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1})$ under Condition \ref{cond:w}. Under Condition \ref{cond:sw}, if we were able to show $\mathbb{L}_{\theta, 4, \widehat{b}, k}^{2} \mathbb{L}_{\theta, 4, \widehat{p}, k}^{2} = O (1)$ when $\mathbb{L}_{\theta, 2, \widehat{b}}, \mathbb{L}_{\theta, 2, \widehat{p}}, \mathbb{L}_{\theta, \infty, \widehat{b}}, \mathbb{L}_{\theta, \infty, \widehat{p}}$ are all $O (1)$, then asymptotic normality would also hold. This remains an open question. For more detailed discussion, see Appendix \ref{app:var}.

Further note the following trade-off between the asymptotic statistical properties and computational cost:
\begin{enumerate}
\item As $m$ grows, the bias of $\widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1})$ decays to 0 at a rate monotonically increasing with $m$.
\item But larger $m$ incurs higher computational cost to compute $\widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1})$.
\end{enumerate}
We conjecture that such a trade-off is unavoidable but a rigorous proof is beyond the scope of this paper.

Finally, the extra condition $k m^{2} = o (n)$ imposes restrictions on both $k$ and $m$ and is a result of the variance bound of $\widehat{\mathbb{IF}}_{mm, k} (\widehat{\Sigma}_{k}^{-1})$ shown in Lemma \ref{lem:hoif_var}. For example, if one chooses $k = n / (\mathsf{log} n)^{2}$, then $m$ is at most $o (\mathsf{log} n)$.
\end{remark}

\subsection{The proposed assumption-lean test of $\mathsf{H}_{0, k} (\delta)$ and early-stopping}\leavevmode
\label{sec:hierarchy}

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$,
\begin{equation}
\label{eq:higher_test}
\widehat{\chi}_{m, k} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta) \coloneqq \mathbbm{1} \left\{ \frac{\vert \widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1}) \vert}{\widehat{\mathsf{s.e.}} [\widehat{\psi}_{1}]} - z_{\alpha^{\dag} / 2} \frac{\widehat{\mathsf{s.e.}} [\widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1})]}{\widehat{\mathsf{s.e.}} [\widehat{\psi}_{1}]} > \delta \right\}.
\end{equation}

Theorem \ref{thm:hoif_stats} immediately implies the following:

\begin{theorem}\leavevmode
\label{prop:hoif_test}
Under the conditions in Theorem \ref{thm:hoif_stats} with Condition \ref{cond:w}, $k \mathsf{log}(k) \vee (k m^{2}) = o(n)$ together with the additional restriction
\begin{equation}\label{eq:add}
\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1}) \neq o \left( \mathbb{L}_{\theta, 2, \widehat{b}, k} \mathbb{L}_{\theta, 2, \widehat{p}, k} \left( \frac{k \mathsf{log} (k)}{n} \right)^{\frac{m - 1}{2}} \right)
\end{equation}
for any given $\delta > 0$, suppose that $\frac{\vert \mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1}) \vert}{\mathsf{s.e.}_{\theta} [\widehat{\psi}_{1}]} = \gamma$ for some (sequence) $\gamma = \gamma (n)$ (where $\gamma (n)$ can diverge with $n$), then the rejection probability of $\widehat{\chi}_{m, k} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta)$ converges to
\begin{equation}
\label{rejection:3}
2 - \Phi \left( z_{\alpha^{\dag} / 2} - \lim_{n \rightarrow \infty} (\gamma - \delta) \frac{\mathsf{s.e.}_{\theta} [\widehat{\psi}_{1}]}{\mathsf{s.e.}_{\theta} [\widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1})]} \right) - \Phi \left( z_{\alpha^{\dag} / 2} + \lim_{n \rightarrow \infty} (\gamma + \delta) \frac{\mathsf{s.e.}_{\theta} [\widehat{\psi}_{1}]}{\mathsf{s.e.}_{\theta} [\widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1})]} \right)
\end{equation}
as $n \rightarrow \infty$. In particular,

\begin{enumerate}[label=(\arabic*)]
\item under $\mathsf{H}_{0, k} (\delta): \gamma \leqslant \delta$, $\widehat{\chi}_{m, k} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta)$ rejects the null with probability less than or equal to $\alpha^{\dag}$, as $n \rightarrow \infty$;

\item under the following alternative to $\mathsf{H}_{0, k} (\delta)$: $\gamma = \delta + c$, for any diverging sequence $c = c(n) \rightarrow \infty$, $\widehat{\chi}_{m, k} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta)$ rejects the null with probability converging to 1, as $n \rightarrow \infty$.
\end{enumerate}

\begin{enumerate}[label=(\arabic*')]
\setcounter{enumi}{2}
\item[(2')] If $\mathbb{L}_{\theta, 2, \widehat{b}, k}$ and $\mathbb{L}_{\theta, 2, \widehat{p}, k}$ converge to 0, under the following alternative to $\mathsf{H}_{0, k} (\delta)$: $\gamma = \delta + c$, for any fixed $c > 0$ or any diverging sequence $c = c(n) \rightarrow \infty$, $\widehat{\chi}_{m, k} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta)$ has rejection probability converging to 1, as $n \rightarrow \infty$.
\end{enumerate}
\end{theorem}

\begin{remark}
Note that the greater $m$ is, the less restrictive \eqref{eq:add} becomes. In fact we can let $m \equiv m (n) \rightarrow \infty$ as $n \rightarrow \infty$, by which \eqref{eq:add} converges to a null requirement.

However, as stated in Remark \ref{rem:comp}, although increasing $m$ relaxes the restriction \eqref{eq:add}, it incurs higher computational cost. In practice, one might want $m$ to be not too large. In view of such computational concern in practice, we propose the following ``early-stopping'' test that terminates at the smallest $m$ at which it fails to reject $\mathsf{H}_{0, k} (\delta)$ and reports ``failure to reject $\mathsf{H}_{0, k} (\delta)$''.

Finally, notice that we only state Theorem \ref{prop:hoif_test} under the stronger Condition \ref{cond:w}. If weakened to Condition \ref{cond:sw}, even under $\mathsf{H}_{0, k} (\delta)$, $k \mathbb{L}_{\theta, 2, \widehat{b}, k}^{2} \mathbb{L}_{\theta, 2, \widehat{p}, k}^{2}$ may still diverge, and the asymptotic normality of $\widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1})$ cannot be justified unless we were able to show $\mathbb{L}_{\theta, 4, \widehat{b}, k}^{2} \mathbb{L}_{\theta, 4, \widehat{p}, k}^{2}$ to be bounded, as discussed in Remark \ref{rem:comp} and Appendix \ref{app:var}. Thus we need to change the cutoff from $z_{\alpha^{\dag} / 2}$ in \eqref{eq:higher_test} to $(\alpha^{\dag})^{- 1 / 2}$ and justify the level of the test using Chebyshev inequality instead.


\end{remark}

Formally, for a fixed integer $M > 0$ that is determined by the analyst's computational budget, define
\begin{equation}
\label{eq:early_test}
\widehat{\chi}_{M, k}^{es} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta) \coloneqq \mathbbm{1} \left\{ \frac{\vert \widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1}) \vert}{\widehat{\mathsf{s.e.}} [\widehat{\psi}_{1}]} - z_{\alpha^{\dag} / 2} \frac{\widehat{\mathsf{s.e.}} [\widehat{\mathbb{IF}}_{22 \rightarrow mm, k} (\widehat{\Sigma}_{k}^{-1})]}{\widehat{\mathsf{s.e.}} [\widehat{\psi}_{1}]} > \delta: \forall \ m = 2, \ldots, M \right\}.
\end{equation}
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{prop:hoif_test}:

\begin{proposition}
\label{prop:early}
Assume all the conditions of Theorem \ref{prop:hoif_test}, but with \eqref{eq:add} replaced by
\begin{equation*}
\mathsf{Bias}_{\theta, k} (\widehat{\psi}_{1}) \neq o \left( \mathbb{L}_{\theta, 2, \widehat{b}, k} \mathbb{L}_{\theta, 2, \widehat{p}, k} \left( \frac{k \mathsf{log} (k)}{n} \right)^{\frac{M - 1}{2}} \right).
\end{equation*}
Under $\mathsf{H}_{0, k} (\delta)$:
\begin{align*}
\mathsf{P}_{\theta} \left( \widehat{\chi}_{M, k}^{es} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta) = 1 \right) \leqslant \alpha^{\dag}
\end{align*}
that is, $\widehat{\chi}_{M, k}^{es} (\widehat{\Sigma}_{k}^{-1}; z_{\alpha^{\dag} / 2}, \delta)$ is an asymptotically valid level-$\alpha^{\dag}$ test of $\mathsf{H}_{0, k} (\delta)$.
\end{proposition}

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

\section{Monte Carlo experiments}
\label{sec:simulation}
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{app:simulations}. 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{app:multiX}. We then simulate $Y$ and $A$ according to the following data generating mechanism:
\begin{equation*}
Y \sim b (X) + N (0, 1) = \sum_{j = 1}^{4} \tau_{b, j} h_{b} (X_{j}; 0.25) + N (0, 1)
\end{equation*}
and
\begin{equation*}
A \sim \mathsf{Bernoulli} \left( 1 / p (X) \equiv \mathsf{expit} \left\{ \sum_{j = 1}^{4} \tau_{p, j} h_{p} (X_{j}; 0.25) \right\} \right)
\end{equation*}
where $h_{b} (\cdot; 0.25)$ and $h_{p} (\cdot; 0.25)$ have the forms as defined in Appendix \ref{app:simulations} 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{tab:s1}. 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{app:simulations} 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
\begin{equation*}
\bar{\mathsf{z}}_{k} (X) = (\bar{\mathsf{z}}_{k'} (X_{1})^{\top}, \bar{\mathsf{z}}_{k'} (X_{2})^{\top}, \bar{\mathsf{z}}_{k'} (X_{3})^{\top}, \bar{\mathsf{z}}_{k'} (X_{4})^{\top})^{\top},
\end{equation*}
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.

\subsection{Finite sample performance of tests for $\mathsf{H}_{0, k} (\delta)$}
\label{sec:sim_bias}

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{fig:qq_ha} (simulation setup I) and \ref{fig:qq_h0} (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{app:bootstrap}. 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{tab:var_ha} (simulation setup I) and \ref{tab:var_h0} (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 \textit{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{tab:var_ha} and Table \ref{tab:var_h0}, 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{tab:main_ha} (simulation setup I) and \ref{tab:main_h0} (simulation setup II).

\begin{itemize}
\item In the upper panel of Table \ref{tab:main_ha} (simulation setup I), the first row and the third column, we display the MC bias of the DML estimator $\widehat{\psi}_{1}$ and the MC average of $\widehat{\mathsf{s.e.}}[\widehat{\psi}_{1}]$, which are $-34.26 \times 10^{-3}$ and $8.77 \times 10^{-3}$. Since the ratio between the bias and standard error is around 4, we expect the associated 90\% Wald CI $\widehat{\psi}_{1} \pm z_{0.05} \widehat{\mathsf{s.e.}}(\widehat{\psi}_{1})$ does not have the nominal coverage. This is indeed the case by reading from the first row, the second column of the upper panel of Table \ref{tab:main_ha}, showing the MC coverage probability for the 90\% Wald CI of $\widehat{\psi}_{1}$ is 0\%.

Similarly, in the upper panel of Table \ref{tab:main_h0} (simulation setup II), the first row and the third column, we display the MC bias of the DML estimator $\widehat{\psi}_{1}$ and the MC average of $\widehat{\mathsf{s.e.}}[\widehat{\psi}_{1}]$, which are $-8.88 \times 10^{-3}$ and $8.10 \times 10^{-3}$. This is as expected because the true nuisance functions $b$ and $p$ belong to $\text{H\"{o}lder} (0.6 + c)$ for some $c > 0$ and DML estimator (with $b$ and $p$ estimated in optimal rate in $L_{2}$ norm) is expected to have bias of $o (n^{-1/2})$ when the average smoothness between $b$ and $p$ is above 0.5 \citep{robins2009semiparametric}. Here the ratio between the bias and standard error is around 1 and we expect the associated 90\% Wald CI $\widehat{\psi}_{1} \pm z_{0.05} \widehat{\mathsf{s.e.}}(\widehat{\psi}_{1})$ is slightly undercovered. This is indeed the case by reading from the first row, the second column of the upper panel of Table \ref{tab:main_h0}, showing the MC coverage probability for the 90\% Wald CI of $\widehat{\psi}_{1}$ is 83\% (note that $\widehat{\psi}_{2, k} = \widehat{\psi}_{1}$ when $k = 0$).

\item In the second column of the upper panel of Table \ref{tab:main_ha}, we display the MC averages of $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$ and the MC averages of its estimated standard errors; in the middle column, we display those of $\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$; and in the lower panel, we display those of $\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$ after third-order bias correction. In the upper panel, we observe that increasing $k$ from 256 to 512 does improve the amount of bias recovered by $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$ (from $-18.76 \times 10^{-3}$ to $- 25.34 \times 10^{-3}$), but there is no obvious difference in the MC averages of $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$ between $k = 512$ and $k = 1024$ ($- 25.34 \times 10^{-3}$ and $- 25.43 \times 10^{-3}$, about 75\% of the total bias). The MC average of the estimated standard error increases with $k$, from $2.45 \times 10^{-3}$ at $k = 256$ to $3.43 \times 10^{-3}$ at $k = 1024$. This is consistent with the theoretical prediction that the variability of $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$ should grow linearly with $k$. In the middle, we observe very similar numerical results between $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$ and $\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$ for all the statistics that we are interested in. Interestingly, we did not see much improvement of using $\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$ instead of $\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$, when compared to $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$, at least in this example.

In the second column of the upper panel of Table \ref{tab:main_h0}, in the upper panel, we also observe that increasing $k$ from 256 to 512 also slightly improve the amount of bias recovered by $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$ (from $-4.54 \times 10^{-3}$ to $- 4.94 \times 10^{-3}$). The MC average of the estimated standard error increases with $k$, from $1.51 \times 10^{-3}$ at $k = 256$ to $2.43 \times 10^{-3}$ at $k = 1024$. In the middle panel, we observe very similar numerical results between $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$ and $\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$ for all the statistics that we are interested in. Interestingly, we again did not see much improvement of using $\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$ instead of $\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$, when compared to $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$, at least in this example.

\item In the third columns of Table \ref{tab:main_ha}, we display the MC coverage probabilities of the 90\% two-sided CIs of $\widehat{\psi}_{2, k} (\Sigma_{k}^{-1})$ (upper panel) and $\widehat{\psi}_{3, k} (\widehat{\Sigma}_{k}^{-1})$ (lower panel). Comparing the oracle bias-corrected estimator $\widehat{\psi}_{2, k} (\Sigma_{k}^{-1})$ vs. $\widehat{\psi}_{1}$, the coverage probability improves from 0\% to 33\% (82\%) at $k = 256$ (at $k = 1024$). Similar observations can be made for $\widehat{\psi}_{2, k} (\widehat{\Sigma}_{k}^{-1})$ and $\widehat{\psi}_{3, k} (\widehat{\Sigma}_{k}^{-1})$ as well when $\Sigma_{k}$ is unknown.

In the third columns of Table \ref{tab:main_h0}, we display the MC coverage probabilities of the 90\% two-sided CIs of $\widehat{\psi}_{2, k} (\Sigma_{k}^{-1})$ (upper panel) and $\widehat{\psi}_{3, k} (\widehat{\Sigma}_{k}^{-1})$ (lower panel). Comparing the oracle bias-corrected estimator $\widehat{\psi}_{2, k} (\Sigma_{k}^{-1})$ vs. $\widehat{\psi}_{1}$, the coverage probability improves from 83\% to 100\% at $k = 256$ and $k = 1024$. Similar observations can be made for $\widehat{\psi}_{2, k} (\widehat{\Sigma}_{k}^{-1})$ and $\widehat{\psi}_{3, k} (\widehat{\Sigma}_{k}^{-1})$ as well when $\Sigma_{k}$ is unknown.

\item In the fourth columns of Tables \ref{tab:main_ha} and \ref{tab:main_h0}, we display the MC biases of $\widehat{\psi}_{2, k} (\Sigma_{k}^{-1})$ (upper panel) and $\widehat{\psi}_{3, k} (\widehat{\Sigma}_{k}^{-1})$ (lower panel), together with their MC standard deviations (in the parentheses). As expected, the standard deviations of the estimators after bias correction are very similar to those of $\widehat{\psi}_{1}$. This indeed confirms that we are able to correct bias without drastically inflating the variance.

\item In the fifth column of Table \ref{tab:main_ha}, we display the MC rejection rates of the test statistic: the upper panel shows the MC rejection rates of the oracle test $\widehat{\chi}_{2, k} (\Sigma_{k}^{-1}; z_{0.10 / 2}, \delta = 3 / 4)$ and the lower panel shows the MC rejection rates of the test $\widehat{\chi}_{3, k} (\widehat{\Sigma}_{k}^{-1}; z_{0.10 / 2}, \delta = 3 / 4)$. These simulation results demonstrate that when a ``good'' dictionary (db6 father wavelets) is used, our proposed test does have power to reject the null hypothesis of actual interest $\mathsf{H}_{0} (\delta)$ that the ratio between the bias of $\widehat{\psi}_{1}$ and its standard error is lower than $\delta$. All the rejection rates in Table \ref{tab:main_ha} are 100\% when $\delta = 3 / 4$, regardless of whether we are using the oracle test or not. When $\delta = 2$, we obtained nontrivial rejection rates (data not shown). This is as expected because the ratio between $\widehat{\mathbb{IF}}_{22, k}$'s and $\widehat{\mathsf{s.e.}} (\widehat{\psi}_{1})$ is around 2 to 3. In the fifth column of Table \ref{tab:main_h0}, all the rejection rates are 0\% when $\delta = 3 / 4$. This is also as expected because the ratio between $\widehat{\mathbb{IF}}_{22, k}$'s and $\widehat{\mathsf{s.e.}} (\widehat{\psi}_{1}) $ is close to $1 / 2$.
\end{itemize}

\begin{remark}
Here we choose $k \lesssim n / (\mathsf{log} n)^{2}$, which is on the order of 1000 when $n = 50,000$. Larger $k$ could easily lead to numerical instability due to inverting a large-dimensional sample Gram matrix (also see \citet{liu2020nearly}). In a technical report by one of the authors \citep{liu2023hoif}, a new class of empirical HOIF estimators is proposed that overcomes the above issues by replacing $\widehat{\Sigma}_{k}$ computed from the training sample by that computed from the estimation sample. The theoretical results are somewhat more complicated to state than the estimators used in this paper. We will thus only refer interested readers to \citet{liu2023hoif}. In a separate manuscript \citep{wanis2023falsification}, we report the finite-sample and real-world data performance of these new estimators and the corresponding test statistics, together with a data-driven method for choosing $k$ in practice. We also refer readers to~\citet{breunig2020adaptive} and \citet{liu2021adaptive} for data-driven methods of choosing $k$ in a theoretically-oriented manner.
\end{remark}

\begin{table}[tbp]
\centering
\begin{tabular}{c|c|c|c}
\hline
$k$ & 256 & 512 & 1024 \\
\hline
$\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$ & \begin{tabular}{@{}c}
2.437 \\
2.268 (0.218)
\end{tabular} & \begin{tabular}{@{}c}
3.011 \\
2.716 (0.268)
\end{tabular} & \begin{tabular}{@{}c}
3.247 \\
2.841 (0.322)
\end{tabular} \\
\hline
$\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$ & \begin{tabular}{@{}c}
2.456 \\
2.271 (0.218)
\end{tabular} & \begin{tabular}{@{}c}
3.075 \\
2.778 (0.278)
\end{tabular} & \begin{tabular}{@{}c}
3.546 \\
3.032 (0.378)
\end{tabular} \\
\hline
$\widehat{\mathbb{IF}}_{33, k} (\widehat{\Sigma}_{k}^{-1})$ & \begin{tabular}{@{}c}
0.495 \\
0.603 (0.0598)
\end{tabular} & \begin{tabular}{@{}c}
0.814 \\
1.126 (0.126)
\end{tabular} & \begin{tabular}{@{}c}
1.518 \\
1.998 (0.285)
\end{tabular} \\
\hline
$\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$ & \begin{tabular}{@{}c}
2.282 \\
2.251 (0.227)
\end{tabular} & \begin{tabular}{@{}c}
2.745 \\
2.738 (0.290)
\end{tabular} & \begin{tabular}{@{}c}
2.743 \\
3.009 (0.427)
\end{tabular} \\
\hline
\end{tabular}
\caption{For data generating mechanism in simulation setup I: We reported the MC standard deviations $\times 10^{-3}$ (upper values in each cell), the MC averages of the estimated standard errors $\times 10^{-3}$ (lower values in each cell outside the parenthesis) and the MC standard deviations of the estimated standard errors $\times 10^{-3}$ (lower values in each cell inside the parenthesis) through nonparametric bootstrap resampling for $\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})$.}
\label{tab:var_ha}
\end{table}

\begin{table}[tbp]
\resizebox{\columnwidth}{!}{
\begin{tabular}{c|c|c|c|c}
\hline
k & $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1}) \times 10^{-3}$ & \shortstack{MC Coverage \\ ($\widehat{\psi}_{2, k} (\Sigma_{k}^{-1})$ 90\% Wald CI)} & $\mathsf{Bias} (\widehat{\psi}_{2, k} (\Sigma_{k}^{-1})) \times 10^{-3}$ & \shortstack{$\widehat{\chi}_{2, k} (\Sigma_{k}^{-1}; z_{0.10 / 2}, \delta = 3 / 4)$} \\
\hline
$0$ & $\mathsf{NA}$ ($\mathsf{NA}$) & 0\% & -34.26 (8.77) & $\mathsf{NA}$ \\
$256$ & -18.76 (2.27) & 31\% & -15.50 (8.60) & 100\% \\
$512$ & -25.34 (2.72) & 83\% & -8.92 (8.61) & 100\% \\
$1024$ & -25.43 (2.84) & 82\% & -8.83 (8.79) & 100\% \\
\hline
\hline
k & $\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1}) \times 10^{-3}$ & \shortstack{MC Coverage \\ ($\widehat{\psi}_{3, k} (\widehat{\Sigma}_{k}^{-1})$ 90\% Wald CI)} & $\mathsf{Bias} (\widehat{\psi}_{3, k} (\widehat{\Sigma}_{k}^{-1})) \times 10^{-3}$ & \shortstack{$\widehat{\chi}_{3, k} (\widehat{\Sigma}_{k}^{-1}; z_{0.10 / 2}, \delta = 3 / 4)$} \\
\hline
$256$ & -18.54 (2.25) & 30\% & -15.72 (8.60) & 100\% \\
$512$ & -24.65 (2.74) & 81\% & -9.61 (8.60) & 100\% \\
$1024$ & -23.82 (3.01) & 76\% & -10.44 (8.76) & 100\% \\
\hline
\end{tabular}}
\caption{For data generating mechanism in simulation setup I: We reported the MC averages of point estimates and standard errors (the first column in each panel) of $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$ (upper panel) and $\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$ (lower panel), together with the coverage probabilities of two-sided 90\% Wald CIs (the second column in each panel) of $\widehat{\psi}_{2, k} (\Sigma_{k}^{-1})$ (upper panel) and $\widehat{\psi}_{3, k} (\widehat{\Sigma}_{k}^{-1})$ (lower panel), the MC biases and MC standard deviations (the third column in each panel) of $\widehat{\psi}_{2, k} (\Sigma_{k}^{-1})$ (upper panel) and $\widehat{\psi}_{3, k} (\widehat{\Sigma}_{k}^{-1})$ (lower panel) and the empirical rejection rates (the fourth column in each panel) of $\widehat{\chi}_{2, k}^{(1)} (\Sigma_{k}^{-1}; z_{0.10 / 2}, \delta = 3 / 4)$ (upper panel) and $\widehat{\chi}_{3, k}^{(1)} (\widehat{\Sigma}_{k}^{-1}; z_{0.10 / 2}, \delta = 3 / 4)$ (lower panel). In the upper panel, for $k = 0$, $\widehat{\psi}_{2, k = 0} \equiv \widehat{\psi}_{1}$.}
\label{tab:main_ha}
\end{table}

\begin{table}[tbp]
\centering
\begin{tabular}{c|c|c|c}
\hline
$k$ & 256 & 512 & 1024 \\
\hline
$\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$ & \begin{tabular}{@{}c}
1.166 \\
1.209 (0.156)
\end{tabular} & \begin{tabular}{@{}c}
1.367 \\
1.408 (0.191)
\end{tabular} & \begin{tabular}{@{}c}
1.673 \\
1.584 (0.267)
\end{tabular} \\
\hline
$\widehat{\mathbb{IF}}_{22, k} (\widehat{\Sigma}_{k}^{-1})$ & \begin{tabular}{@{}c}
1.187 \\
1.219 (0.159)
\end{tabular} & \begin{tabular}{@{}c}
1.433 \\
1.439 (0.198)
\end{tabular} & \begin{tabular}{@{}c}
1.828 \\
1.697 (0.285)
\end{tabular} \\
\hline
$\widehat{\mathbb{IF}}_{33, k} (\widehat{\Sigma}_{k}^{-1})$ & \begin{tabular}{@{}c}
0.251 \\
0.331 (0.0451)
\end{tabular} & \begin{tabular}{@{}c}
0.417 \\
0.629 (0.0997)
\end{tabular} & \begin{tabular}{@{}c}
0.854 \\
1.251 (0.247)
\end{tabular} \\
\hline
$\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$ & \begin{tabular}{@{}c}
1.109 \\
1.186 (0.156)
\end{tabular} & \begin{tabular}{@{}c}
1.256 \\
1.346 (0.191)
\end{tabular} & \begin{tabular}{@{}c}
1.488 \\
1.560 (0.383)
\end{tabular} \\
\hline
\end{tabular}
\caption{For data generating mechanism in simulation setup II: We reported the MC standard deviations $\times 10^{-3}$ (upper values in each cell), the MC averages of the estimated standard errors $\times 10^{-3}$ (lower values in each cell outside the parenthesis) and the MC standard deviations of the estimated standard errors $\times 10^{-3}$ (lower values in each cell inside the parenthesis) through nonparametric bootstrap resampling for $\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})$.}
\label{tab:var_h0}
\end{table}

\begin{table}[tbp]
\resizebox{\columnwidth}{!}{
\begin{tabular}{c|c|c|c|c}
\hline
k & $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1}) \times 10^{-3}$ & \shortstack{MC Coverage \\ ($\widehat{\psi}_{2, k} (\Sigma_{k}^{-1})$ 90\% Wald CI)} & $\mathsf{Bias} (\widehat{\psi}_{2, k} (\Sigma_{k}^{-1})) \times 10^{-3}$ & \shortstack{$\widehat{\chi}_{2, k} (\Sigma_{k}^{-1}; z_{0.10 / 2}, \delta = 3 / 4)$} \\
\hline
$0$ & $\mathsf{NA}$ ($\mathsf{NA}$) & 83\% & -8.88 (8.10) & $\mathsf{NA}$ \\
$256$ & -4.54 (1.21) & 99\% & -4.34 (8.13) & 0\% \\
$512$ & -4.89 (1.41) & 99\% & -3.99 (8.21) & 0\% \\
$1024$ & -4.94 (1.58) & 100\% & -3.94 (8.35) & 0\% \\
\hline
\hline
k & $\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1}) \times 10^{-3}$ & \shortstack{MC Coverage \\ ($\widehat{\psi}_{3, k} (\widehat{\Sigma}_{k}^{-1})$ 90\% Wald CI)} & $\mathsf{Bias} (\widehat{\psi}_{3, k} (\widehat{\Sigma}_{k}^{-1})) \times 10^{-3}$ & \shortstack{$\widehat{\chi}_{3, k} (\widehat{\Sigma}_{k}^{-1}; z_{0.10 / 2}, \delta = 3 / 4)$} \\
\hline
$256$ & -4.49 (1.19) & 99\% & -4.39 (8.13) & 0\% \\
$512$ & -4.76 (1.35) & 99\% & -4.12 (8.20) & 0\% \\
$1024$ & -4.62 (1.56) & 100\% & -4.26 (8.32) & 0\% \\
\hline
\end{tabular}}
\caption{For data generating mechanism in simulation setup II: We reported the MC averages of point estimates and standard errors (the first column in each panel) of $\widehat{\mathbb{IF}}_{22, k} (\Sigma_{k}^{-1})$ (upper panel) and $\widehat{\mathbb{IF}}_{22 \rightarrow 33, k} (\widehat{\Sigma}_{k}^{-1})$ (lower panel), together with the coverage probabilities of two-sided 90\% Wald CIs (the second column in each panel) of $\widehat{\psi}_{2, k} (\Sigma_{k}^{-1})$ (upper panel) and $\widehat{\psi}_{3, k} (\widehat{\Sigma}_{k}^{-1})$ (lower panel), the MC biases and MC standard deviations (the third column in each panel) of $\widehat{\psi}_{2, k} (\Sigma_{k}^{-1})$ (upper panel) and $\widehat{\psi}_{3, k} (\widehat{\Sigma}_{k}^{-1})$ (lower panel) and the empirical rejection rates (the fourth column in each panel) of $\widehat{\chi}_{2, k}^{(1)} (\Sigma_{k}^{-1}; z_{0.10 / 2}, \delta = 3 / 4)$ (upper panel) and $\widehat{\chi}_{3, k}^{(1)} (\widehat{\Sigma}_{k}^{-1}; z_{0.10 / 2}, \delta = 3 / 4)$ (lower panel). In the upper panel, for $k = 0$, $\widehat{\psi}_{2, k = 0} \equiv \widehat{\psi}_{1}$.}
\label{tab:main_h0}
\end{table}

\section{Concluding remarks}
\label{sec:conclusion}
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{sec:main}). 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 \citep{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 \citep{chen2023efficient, kompa2022deep}. In Appendix \ref{app:prox_hoif}, we document the HOIFs for the running example $\psi (\theta) = - \mathsf{E}_{\theta} [Y (a = 1)]$ in Section \ref{sec:introduction} under the so-called ``proximal causal learning'' framework \citep{tchetgen2020introduction}. This framework has been shown to be closely related to other approaches dealing with endogeneity, e.g. synthetic controls \citep{abadie2010synthetic, shi2021theory} or quadratic functionals of nonparametric instrumental variable regression \citep{breunig2019simple}. Based on the form of these HOIFs, it is straightforward to generalize our method to scenarios under endogeneity. Second, as mentioned earlier in \textbf{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 \citep{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 \citep{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 \citet{cattaneo2018kernel, cattaneo2019two}, and \citet{koltchinskii2020estimation}. Finally, we are also considering to directly learn a data-driven representation $\widetilde{\bar{\mathsf{z}}}_{k}$ via the penultimate layer \citep{ansuini2019intrinsic, damian2022neural} or distillation \citep{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

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

\allowdisplaybreaks