New $\sqrt{n}$-consistent, numerically stable higher-order influence function estimators

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

100,496 characters · 21 sections · 118 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

New $n$-consistent, numerically stable higher-order influence function estimators

\affil[2]{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 and Shanghai Artificial Intelligence Laboratory, Shanghai, China} \affil[1]{Department of Statistics, University of Virginia, Charlottesville, VA, USA}

abstractHigher-Order Influence Functions (HOIFs) provide a unified theory for constructing rate-optimal estimators for a large class of low-dimensional (smooth) statistical functionals/parameters (and sometimes even infinite-dimensional functions) that arise in substantive fields including epidemiology, economics, and the social sciences. Since the introduction of HOIFs by robins2008higher or robins2016technical\footnote{robins2016technical is the complete version of robins2008higher, including more results and proofs. We therefore only refer to robins2016technical in the sequel.}, they have been viewed mostly as a theoretical benchmark rather than a useful tool for statistical practice. Works aimed to flip the script are scant, but a few recent papers liu2017semiparametric, liu2021assumption make some partial progress. In this paper, we take a fresh attempt at achieving this goal by constructing new, numerically stable HOIF estimators (or sHOIF estimators for short with “s” standing for “stable”) with provable statistical and computational guarantees. This new class of sHOIF estimators (up to the 2nd order) was foreshadowed in synthetic experiments conducted by liu2020nearly.

{ Keywords: Causal Inference, Functional Estimation, Higher-Order Influence Functions, Semiparametric Theory, Combinatorics}

Introduction

{\it Higher-Order Influence Functions} (HOIFs) robins2016technical are higher-order generalizations of the first-order influence functions (IFs), a staple in semiparametric statistical theory newey1990semiparametric, bickel1998efficient, van2002part. HOIFs are a powerful and unified approach to constructing minimax rate-optimal estimators for a class of statistical functionals/parameters (and sometimes even functions; see kennedy2022minimax) that arise in (bio)statistics, epidemiology, economics, and the social sciences. HOIF estimators originally proposed in robins2016technical, robins2017minimax\footnote{See robins2022corrigenda for corrections of the proofs in robins2017minimax.} remain the only known minimax rate-optimal estimators for statistical functionals/parameters with substantive interests in the above disciplines, including the Average Treatment Effect (ATE) under the strong ignorability assumption\footnote{In liu2021assumption, we derived the HOIFs for the ATE functional even when the strong ignorability assumption fails to hold, provided that we have access to valid proxies for both the treatment and outcome, following a series of works on proximal causal learning tchetgen2020introduction.} and the expected conditional covariance of two random variables $A$ and $Y$ given a third random variable $X$, even after highly active research by the statistics and econometrics communities in recent years newey2018cross, kennedy2020optimal, hirshberg2021augmented, yu2020treatment. More recent works kennedy2022minimax, bonvini2022fast also initiated the application of the HOIF machinery to the minimax optimal estimation of Conditional Average Treatment Effect (CATE) function or dose response curves. Their results lay important theoretical foundation for individualized decision making problems, e.g. personalized medicine. This is the first instance when HOIF estimators are also shown to be effective, at least in theory, for function estimation problems, or more precisely, “hybrid function and functional estimation problems”. Similar idea has also been applied to dose-response curve estimation bonvini2022fast. For an introductory level review of HOIFs, we refer the interested readers to van2014higher and Section 1 of liu2020rejoinder. A relatively more technical review of HOIFs is delegated to Section (ref).

Over the past decade, Robins and colleagues initiated the research program of establishing theoretical foundations for HOIFs and estimators based on HOIFs robins2004optimal, van2014higher, robins2016technical, robins2017minimax, liu2017semiparametric for a class of statistical functionals/parameters recently characterized in rotnitzky2021characterization, which are heretofore termed as {\it Doubly Robust Functionals} (DRF) in this paper. We adopt this terminology to reflect the fact that their nonparametric first-order IFs give rise to doubly robust estimators scharfstein1999adjusting, robins2001comments, chernozhukov2018double. This class of DRFs subsumes the class of functionals studied in robins2016technical and chernozhukov2018riesz. Under the standard H\"{o}lder-regularity assumptions on the nuisance parameters (abbreviated as H\"{o}lder nuisance models), robins2017minimax constructed minimax optimal but non-adaptive HOIF estimators for a sub-class of DRFs. liu2021adaptive constructed adaptive second-order IF estimators for DRFs using the celebrated Lepski\v{i}'s adaptation scheme lepskii1991problem, within a strict submodel of the H\"{o}lder nuisance models. But both estimators require estimating the density of the potentially high-dimensional covariates $X$, even in $\sqrt{n}$-estimable regimes. When the dimension $d$ of the covariates is only moderately large (e.g. $d = 10$), nonparametric density estimation is already a daunting computational and statistical task.

To overcome the above issues, liu2017semiparametric introduced {\it empirical} HOIF (eHOIF for short) estimators that obviate multi-dimensional density estimation by inverting the sample/empirical Gram matrix of vector-valued basis transformation of the covariates computed using a separate sample independent of the sample used to construct the estimator of the DRF. This sample-splitting strategy is adopted mainly for simplifying the mathematical analysis, leading to rather straightforward analysis of the statistical properties of the eHOIF estimators. In particular, the eHOIF estimators are still the only class of estimators that achieves $\sqrt{n}$-consistency and semiparametric efficiency for DRFs under the minimal H\"{o}lder-regularity assumptions robins2009semiparametric. These nice statistical properties of the eHOIF estimators also motivate the development of a class of assumption-lean hypothesis tests statistic that is designed to falsify if the standard $(1 - \alpha) \times 100\%$ Wald confidence interval of the DRF has the claimed coverage probability liu2020nearly, liu2021assumption. At this point, astute readers must wonder why we need a new class of empirical HOIF estimators at all, which is what this article is all about.

Motivation and main contributions

Despite the effort in liu2017semiparametric, from our past experience of using eHOIF estimators in practice liu2017semiparametric, liu2020nearly, liu2021assumption, wanis2023machine, several singular issues of their finite-sample performance were unveiled by large-scale simulation experiments\footnote{For interested readers, these simulation experiments have also been used to expose the gap between the (nonparametric) statistical theory deep neural networks (DNNs) and their practice in xu2022deepmed. One can access computer codes of generating such simulations \href{https://github.com/siqixu/DeepMed}{here}.}:

enumerate[label = (\roman*)] • Numerical instability: In liu2017semiparametric, although eHOIF estimators exhibit better finite-sample performance than the original HOIF estimators in robins2017minimax, the simulations were restricted to very low condition number $k / n$: e.g. $k \approx 1,000$ and $n \approx 10,000$. In the simulation studies of liu2020nearly, when $k$ gets near $n$, eHOIF estimators blow up numerically (see Section S3.1 of liu2020nearly) already at order two. What is more striking is that the eHOIF estimators at higher orders, though supposed to be correcting the bias, can only exacerbate the numeric blow-up. • Non-monotone bias reduction: Theoretical results in liu2017semiparametric hint that increasing the orders of the estimator should {\it in principle} reduce the bias. However, we found that this is not usually the case for eHOIF estimators in practice (e.g. see Section 5 of liu2021assumption). Interestingly, sHOIF estimators do not seem to suffer from this problem in simulations, elevating the theoretical results from {\it mere principles} closer to {\it empirical facts}; see liu2020nearly or wanis2023machine for simulations at orders 2 or 3.

Our contributions are three-fold.

itemize• Methodology and practical relevance: This article proposes a new class of numerically stable sHOIF estimators for DRFs, that overcomes the above two major limitations of eHOIF estimators. The stable Second-Order IF (SOIF) estimators first appeared in the simulation studies of liu2020nearly, but their statistical properties remain elusive. • Theory and the proof strategy: Obtaining a deeper theoretical underpinning of this phenomenon mandates meticulous calculations rather than crude upper bounds. This is the critical technical innovation vis-\`{a}-vis other HOIF-related works. In particular, we intensively use the following proof techniques: leave-out analysis, matrix-valued Taylor expansion, and combinatorial calculations (i.e. corollaries of the binomial identity). The proof strategy developed in this paper may be of independent interest. • Extensions of sHOIFs beyond ATE settings: We also generalize sHOIF estimators to all the DRFs, allowing us to handle more structural parameters in the current causal inference (or econometrics) literature.

Notation

Before proceeding, we gather some frequently used notation throughout the paper. We denote the observed data random vector as $O \in {\mathcal{O}}$, where ${\mathcal{O}}$ is its corresponding sample space. Let $\bar{{\mathsf{z}}}_{k} \coloneqq (z_{1}, \cdots, z_{k})^{\top}$ denote a collection of $k$ different functions, each of which has input domain ${\mathcal{X}}$. Fix some $\theta' \in \Theta$. ${\mathbb{E}}_{\theta'}$, $\mathsf{var}_{\theta'}$, and $\mathsf{cov}_{\theta'}$ are, respectively, the expectation, variance, and covariance operators under the probability law ${\mathbb{P}}_{\theta'}$. For any measurable function $h: {\mathcal{X}} \rightarrow {\mathbb{R}}$, let $\Vert h \Vert_{\infty} \coloneqq \mathrm{ess } \sup_{x \in {\mathcal{X}}} h (x)$ and $\Vert h \Vert_{\theta', p} \coloneqq \left\{ {\mathbb{E}}_{\theta'} \left[ h (X)^{p} \right] \right\}^{1 / p}$ for any $p \geq 1$. We adopt standard (stochastic) asymptotic notation $\lesssim$, $\gtrsim$, $\asymp$, $\gg$, $\ll$, $o (\cdot)$, $\omega (\cdot)$, $O (\cdot)$, $\Omega (\cdot)$, $o_{{\mathbb{P}}_{\theta'}} (\cdot)$, $O_{{\mathbb{P}}_{\theta'}} (\cdot)$. For any real-valued vector $\bar{v}$ and any $q \in {\mathbb{R}}$, let $v^{q}$ be the element-wise $q$-th power of $v$.

Furthermore, define $\Sigma_{\theta'} \coloneqq {\mathbb{E}}_{\theta'} [Q]$, where $Q \coloneqq A \bar{{\mathsf{z}}}_{k} (X) \bar{{\mathsf{z}}}_{k} (X)^{\top}$, as the ($A$-weighted) population Gram matrix of $\bar{{\mathsf{z}}}_{k} (X)$, till Section (ref), in which we generalize all our results from $\psi (\theta) \coloneqq {\mathbb{E}}_{\theta} [Y (a = 1)]$\footnote{We use the potential outcome notation without introducing it, which will not affect the understanding of the main theme of this work.}, the mean of outcome $Y$ in the treated group under strong ignorability, to all members of the DRFs. Similarly, define $\widehat{\Sigma} \coloneqq {\mathbb{P}}_{n} [Q_{k}] \equiv n^{-1} \sum_{i = 1}^{n} Q_{i}$ as the ($A$-weighted) sample Gram matrix of $\bar{{\mathsf{z}}}_{k} (X)$, again till Section (ref). Here ${\mathbb{P}}_{n} [\cdot]$ denotes the sample mean operator. To further lighten the notation, we let $Q_{I} \coloneqq \sum_{i \in I} Q_{i}$ for any multi-index set $I \subseteq [n]$. For convenience, we also denote multi-index set $\{i_{1}, i_{2}, \cdots, i_{j}\} \subseteq [n]$ as $\bar{i}_{j}$ for $j \leq n$. $\Omega_{\theta'} \equiv \Sigma_{\theta'}^{-1}$ and $\widehat{\Omega} \equiv \widehat{\Sigma}^{-1}$, when they exist, are respectively the inverse of the population and sample Gram matrices. The kernels constructed from $\bar{{\mathsf{z}}}_{k}$ are denoted as $K_{\theta', k} (x, x') \coloneqq \bar{{\mathsf{z}}}_{k} (x)^{\top} \Omega_{\theta'} \bar{{\mathsf{z}}}_{k} (x')$ and $\widehat{K}_{k} (x, x') = \bar{{\mathsf{z}}}_{k} (x)^{\top} \widehat{\Omega} \bar{{\mathsf{z}}}_{k} (x')$. Given a set of functions $\bar{{\mathsf{v}}}_{k}: {\mathcal{X}} \rightarrow {\mathbb{R}}^{k}$ and any $L_{2}$ function $h: {\mathcal{O}} \rightarrow {\mathbb{R}}$, $\Pi [h | \bar{{\mathsf{v}}}_{k}]$ denotes the linear projection operator of projecting $h$ onto the linear span of $\bar{{\mathsf{v}}}_{k}$: formally, $$ \Pi [h | \bar{{\mathsf{v}}}_{k}] (x) \coloneqq \bar{{\mathsf{v}}}_{k} (x)^{\top} {\mathbb{E}} [\bar{{\mathsf{v}}}_{k} (X) \bar{{\mathsf{v}}}_{k} (X)^{\top}]^{-1} {\mathbb{E}} [\bar{{\mathsf{v}}}_{k} (X) h (O)]. $$

We use ${\mathbb{P}}_{\theta}$ to denote the true data generating law, unless stated otherwise. When the reference measure is the true law ${\mathbb{P}}_{\theta}$, we often drop the dependence on $\theta$: for example, we write $\Vert h \Vert_{p} \equiv \Vert h \Vert_{\theta, p}$, $\Sigma \equiv \Sigma_{\theta}$, $\Omega \equiv \Omega_{\theta}$ and ${\mathbb{P}}$, ${\mathbb{E}}$, $\mathsf{var}$, $\mathsf{cov}$ correspond to ${\mathbb{P}}_{\theta}$, ${\mathbb{E}}_{\theta}$, $\mathsf{var}_{\theta}$, $\mathsf{cov}_{\theta}$. Note that $\Omega$ should not be confused with the asymptotic notation $\Omega (\cdot)$ and this will be clear from the context. A statistic is said to be “oracle” whenever it depends on some part(s) of the unknown {\it true} data generating law ${\mathbb{P}}_{\theta}$ (such as $\Omega$); otherwise it is said to be “feasible”. We also introduce $\mathsf{Diag}$ as the operator of extracting the diagonal elements of a matrix.

Finally, let ${\mathbb{U}}_{n, m} [\cdot]$ denote the $m$-th order $U$-statistic operator: for any function $h: {\mathbb{R}}^{m} \rightarrow {\mathbb{R}}$

align*[align* omitted — 186 chars of source]

When $m = 1$, ${\mathbb{U}}_{n, m} [\cdot]$ reduces to the sample mean operator ${\mathbb{P}}_{n} [\cdot]$. Similarly, let ${\mathbb{V}}_{n, m}$ be the corresponding $V$-statistic operator\footnote{Here we use the scaling $\frac{(n - m)!}{n!}$ instead of the more conventional $\frac{1}{n^{m}}$ for notational convenience.}:

align*[align* omitted — 178 chars of source]

Later in the paper, for $m \geq 2$, we will define “oracle” $m$-th order influence function estimators constructed using the dictionary $\bar{{\mathsf{z}}}_{k}$, denoted as $\widehat{\mathbb{IF}}_{m, m, k} (\Omega) \equiv \widehat{\mathbb{IF}}_{m, m, k} = {\mathbb{U}}_{n, m} [\widehat{\mathsf{IF}}_{m, m, k, \bar{i}_{m}}]$ with $U$-statistic kernel $\mathsf{IF}_{m, m, k, \bar{i}_{m}} \equiv \widehat{\mathsf{IF}}_{m, m, k, \bar{i}_{m}} (\Omega)$. Its stable feasible version is denoted by $\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega})$ with the corresponding kernel $\widehat{\mathsf{IF}}_{m, m, k, \bar{i}_{m}} (\widehat{\Omega})$.

The setup and a review of the theory of HOIFs

With the notation just introduced, we are poised to state the problem setup and briefly review the theory of HOIFs relevant for this paper, in particular the theory of eHOIFs.

Suppose that we are given $N$ i.i.d. observations $\{O_{i}\}_{i = 1}^{N} \sim {\mathbb{P}}_{\theta}$, where $\theta \in \Theta$ is the so-called nuisance parameter and $\Theta$ is its underlying parameter space. Let ${\mathcal{P}} \coloneqq \left\{ {\mathbb{P}}_{\theta}: \theta \in \Theta \right\}$ be the space of data generating probability measures. Our primary interest is to estimate and draw statistical inference on a smooth statistical functional $\psi (\theta): \rightarrow {\mathbb{R}}$, in the sense of van1991differentiable. We restrict $\psi (\theta)$ to be the DRFs defined in rotnitzky2021characterization. As mentioned, our running example is $\psi (\theta) = {\mathbb{E}}_{\theta} [Y (a = 1)]$ the mean of an outcome $Y$ in the treated group $A = 1$. Here the observed data specializes to $O = (X, A, Y)$: respectively the $d$-dimensional covariates belonging to a compact subset ${\mathcal{X}} \equiv [-B, B]^{d}$ of ${\mathbb{R}}^{d}$, the binary treatment assignment, and the bounded outcome variable. Under unconfoundedness assumption (that can be relaxed by using the HOIFs of $\psi (\theta)$ under the proximal causal inference setting liu2021assumption), $\psi (\theta)$ can be identified by either of the two statistical functionals of the observed data distribution:

equation[equation omitted — 126 chars of source]

where $a (x) \coloneqq \{{\mathbb{E}} [A | X = x]\}^{-1}$ and $b (x) \coloneqq {\mathbb{E}} [Y | X = x, A = 1]$ except Section (ref). For this functional $\psi (\theta)$, the nuisance parameter is $\theta \equiv (a, b, g)$ where $g (x)$ is the probability density/mass function of the covariates $X$ conditional on $A = 1$. Hence the nuisance parameter space $\Theta = {\mathcal{A}} \times {\mathcal{B}} \times {\mathcal{G}}$, where ${\mathcal{A}}, {\mathcal{B}}, {\mathcal{G}}$ are, respectively, the space where $a, b, g$ lie. We further divide the whole $N$ data points into two parts: one with sample size $n$, called the estimation sample, and the other with sample size $N - n$, called the nuisance sample used to estimate the nuisance parameter $\theta$. Throughout this paper, we condition on the nuisance sample data by treating it or any quantity computed from it as fixed.

For a smooth statistical functional $\psi (\theta): \Theta \rightarrow {\mathbb{R}}$ in the sense of van1991differentiable, its first-order influence function $\mathbb{IF}_{1} (\theta)$ is a mean-zero first-order $U$-statistic satisfying the following functional equation

align*[align* omitted — 178 chars of source]

where ${\mathbb{P}}_{\theta_{t}}$ is any parametric submodels in $\{{\mathbb{P}}_{\theta}, \theta \in \Theta\}$, such that when $t = 0$, ${\mathbb{P}}_{\theta_{t}} \equiv {\mathbb{P}}_{\theta}$, the true data generating law, and ${\mathbb{S}}_{1}$ is its first-order score vector, as defined in waterman1996projected; also see robins2016technical. Here $\mathbb{IF}_{1} (\theta)$ has the following form robins1994estimation:

equation[equation omitted — 216 chars of source]

Typically, classical semiparametric theory newey1990semiparametric, bickel1998efficient constructs semiparametric efficient first-order estimators $\widehat{\psi}_{1}$ of $\psi (\theta)$ based on its first-order influence function follows: $$ \widehat{\psi}_{1} = \frac{1}{n} \sum_{i = 1}^{n} A_{i} \widehat{a} (X_{i}) (Y_{i} - \widehat{b} (X_{i})) + \widehat{b} (X_{i}) $$ where $\widehat{a}, \widehat{b}$ are nuisance parameter estimates computed from the nuisance sample. In particular, $\widehat{\psi}_{1}$ has bias

equation[equation omitted — 369 chars of source]

Formally, $\mathsf{bias} (\widehat{\psi}_{1})$ is a {\it product of two nuisance estimation errors}\footnote{rotnitzky2021characterization actually define the general class of statistical functionals that permit doubly-robust estimators based on this second-order bias property; see Section (ref).}, and hence {\it doubly-robust} scharfstein1999adjusting.

Despite being doubly-robust, the veracity of inference based on first-order estimators like $\widehat{\psi}_{1}$ may nonetheless be questionable when the nuisance parameter $\theta$ is of high complexity: e.g. functions with low smoothness or without sparsity. For example, when $a, b$ belong to H\"{o}lder functions with smoothness $s_{a}, s_{b}$ and $g$ arbitrarily complex, by far no first-order estimators are known to be $\sqrt{n}$-consistency for estimating $\psi (\theta)$ throughout the entire range

equation[equation omitted — 94 chars of source]

but the eHOIF estimators of liu2017semiparametric or the original HOIF estimators of robins2016technical if additionally assuming $g$ to be H\"{o}lder with smoothness $s_{g} > 0$. In fact, robins2009semiparametric also showed that (ref) is the minimal condition for the existence of $\sqrt{n}$-consistent estimators of $\psi (\theta)$ under the H\"{o}lder nuisance modeling assumption. Outside (ref), $\psi (\theta)$ is non $\sqrt{n}$-estimable and the only known estimator with the optimal rate of convergence in minimax sense is again the HOIF estimator robins2016technical, robins2017minimax, robins2022corrigenda. When restricting to highly smooth $g$, liu2021adaptive construct minimax optimal and adaptive estimator of $\psi (\theta)$ by combining the HOIF estimators with the celebrated Lepskii's adaptation scheme lepskii1991problem.

This article is about the $\sqrt{n}$-estimable regime (ref), so we will focus our attention on the eHOIF estimators. First, we choose a set of $k$-dimensional functions $\bar{{\mathsf{z}}}_{k} \equiv (z_{1}, \cdots, z_{k})^{\top}: {\mathcal{X}} \rightarrow {\mathbb{R}}^{k}$ satisfying certain regularity conditions to be given later in Section (ref). The Second-Order Influence Function (SOIF) estimator of $\psi (\theta)$ is the following second-order $U$-statistic:

equation[equation omitted — 338 chars of source]

and

align*[align* omitted — 381 chars of source]

Based on the definition of HOIFs robins2016technical, $- \widehat{\mathbb{IF}}_{2, 2, k}$ is in fact the SOIF of $\mathsf{bias} (\widehat{\psi}_{1})$\footnote{The difference in the signs in $\widehat{\mathbb{IF}}_{2, 2, k}$ between here and robins2016technical is non-essential.}. A more intuitively appealing explanation goes as follows: $- \widehat{\mathbb{IF}}_{2, 2, k}$ is an unbiased estimator of the following quantity:

equation[equation omitted — 477 chars of source]

which is simply replacing the estimation errors $\widehat{a} / a - 1$ and $b - \widehat{b}$ in (ref) by $$ \Pi \left[ \left. \frac{\widehat{a}}{a} - 1 \right\vert \bar{{\mathsf{z}}}_{k} \right] \text{ and } \Pi \left[ \left. b - \widehat{b} \right\vert A \bar{{\mathsf{z}}}_{k} \right]. $$ Hence $\widehat{\mathbb{IF}}_{2, 2, k}$ can be interpreted as a bias correction term that {\it partially} debiases $\mathsf{bias} (\widehat{\psi}_{1})$.

However, evaluating $\Omega$ in practice relies on the knowledge of $g$, which is generally unknown to the analyst. The initial attempt by robins2016technical and robins2017minimax was to estimate $g$ from the nuisance sample by $\widehat{g}$, leading to statistical properties affected by $g - \widehat{g}$ and thus complexity-reducing assumptions on ${\mathcal{G}} \ni g$. To completely resolve this reliance, liu2017semiparametric choose to estimate $\Omega$ by its empirical analogue using the {\it nuisance sample}, denoted as $\widehat{\Omega}_{\mathrm{nuis}} = \widehat{\Sigma}_{\mathrm{nuis}}^{-1}$. The resulting estimated kernel is denoted as $\widehat{K}_{k}^{\mathrm{nuis}} (x, x')$, similar to $\widehat{K}_{k}$ defined in Section (ref). Then the empirical SOIF (eSOIF) estimator $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}_{\mathrm{nuis}})$ of $\mathsf{bias}_{k} (\widehat{\psi}_{1})$ is

align*[align* omitted — 194 chars of source]

which, unlike $\widehat{\mathbb{IF}}_{2, 2, k}$, incurs a kernel estimation bias

align*[align* omitted — 525 chars of source]

shown to be of order at most $\sqrt{k \log k / n}$ in liu2017semiparametric. To further reduce the kernel estimation bias, one can consider the following $m$-th order eHOIF estimator, which is an $m$-th order $U$-statistic:

align*[align* omitted — 421 chars of source]

and

align*[align* omitted — 395 chars of source]

liu2017semiparametric showed that the kernel estimation bias of $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}_{\mathrm{nuis}})$ is of order at most $(k \log k / n)^{m / 2}$ and variance of order at most $1 / n \vee k / n^{2}$. Hence by taking $m \asymp \sqrt{\log n}$ and $k \asymp n / \log^{c} n$ for some absolute constant $c > 0$, we could estimate $\mathsf{bias}_{k} (\widehat{\psi}_{1})$ with essentially no bias without inflating the order of the variance of $\widehat{\psi}_{1}$. Furthermore, under H\"{o}lder nuisance models on ${\mathcal{A}} \times {\mathcal{B}}$, liu2017semiparametric demonstrate that the sHOIF estimator $\widehat{\psi}_{1} + \widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}_{\mathrm{nuis}})$, with said choices of $m$ and $k$, is $\sqrt{n}$-consistent in (ref) and semiparametric efficient in the interior of (ref) under some additional mild assumptions. In this paper, the sHOIF estimators to be introduced in Section (ref) simply replace $\widehat{\Sigma}_{\mathrm{nuis}}$ and $\widehat{\Omega}_{\mathrm{nuis}}$ in the eHOIF estimators by $\widehat{\Sigma}$ and $\widehat{\Omega}$, the empirical analogues of $\Sigma$ and $\Omega$ computed from the estimation sample. One can easily see that, due to the correlation between $\widehat{\Omega}$ and the estimation sample, the analysis of the statistical properties of sHOIF estimators becomes significantly more challenging.

Plan

The rest of the paper is organized as follows. Section (ref) defines the stable Second-Order IF (sSOIF) estimators and studies their statistical and numerical properties as a warm-up. Section (ref) presents the full version of sHOIF estimators, together with their statistical, numerical, and computational properties. We then apply sHOIF estimators and their statistical properties to two concrete problems Section (ref): one is to show that sHOIF estimators for $\psi (\theta)$ achieve semiparametric efficiency under the minimal conditions within the classical H\"{o}lder nuisance models; the other is to use sHOIF estimators to test if the nominal $(1 - \alpha) \times 100\%$ Wald confidence interval centered at the first-order DML estimator has the claimed coverage, a novel assumption-lean statistical procedure recently proposed in liu2020nearly, and further developed in liu2021assumption. To demonstrate the generality of sHOIF estimators, Section (ref) extends results heretofore in several directions. Finally, Section (ref) concludes the paper and discusses several open problems and possible future directions. Appendix contains technical details that provide insights on the proof strategy. The remaining technical details are deferred to Supplementary Materials shoif_supp.

Assumptions and warm-up: Stable second-order influence function estimators

In this section, we disclose the main assumptions, accompanied with an illustration of the main results using the stable second-order influence function (sSOIF) estimator $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ as a warm-up of what follows.

The assumptions below are imposed throughout the paper unless stated otherwise.

assumption[Conditions on initial first-step nuisance parameter estimates.] Nuisance parameter estimators $\widehat{a}$ and $\widehat{b}$ are attained from a separate independent frozen nuisance sample. For simplicity, we assume this sample to also have size $n$. $\widehat{a}$ and $\widehat{b}$ further satisfy the following properties until otherwise noticed: \begin{enumerate}[label = (\roman*)] • $\Vert \widehat{a} - a \Vert_{2} = o (1)$ and $\Vert \widehat{b} - b \Vert_{2} = o (1)$, i.e. both nuisance parameter estimators are $L_{2}$-consistent; • $\Vert a \Vert_{\infty}$, $\Vert b \Vert_{\infty}$, $\Vert \widehat{a} \Vert_{\infty}$ and $\Vert \widehat{b} \Vert_{\infty}$ are bounded by some absolute constant $B > 0$. • In the case of $\psi (\theta) = {\mathbb{E}} [Y (1)]$ under strong ignorability, we additionally need $1 / a$ and $1 / \widehat{a}$ to be bounded between $(c, 1 - c)$ for some absolute constant $0 < c < 0.5$. \end{enumerate}

assumption[Conditions on $\bar{{\mathsf{z}}}_{k}$ related quantities.] The following are assumed on the basis functions $\bar{{\mathsf{z}}}_{k}$ and the corresponding (inverse) Gram matrices $\Sigma, \widehat{\Sigma}, \Omega, \widehat{\Omega}$ and projection kernels $K_{k}$ and $\widehat{K}_{k}$: \begin{enumerate}[label = (\roman*)] • There exists an absolute constant $B > 0$ such that $\sup_{x \in {\mathcal{X}}} K_{k} (x, x) \leq B k$ and $\sup_{x \in {\mathcal{X}}} \widehat{K}_{k} (x, x) \leq B k$; • Both $\Sigma$ and $\widehat{\Sigma}$ have bounded spectra; • The projection kernel satisfies the following $L_{\infty}$-stability condition: for any measurable function $h: {\mathcal{X}} \rightarrow {\mathbb{R}}$, \begin{equation} \left\Vert \Pi \left[ h | \bar{{\mathsf{z}}}_{k} \right] (\cdot) \right\Vert_{\infty} \lesssim \Vert h \Vert_{\infty}. \end{equation} \end{enumerate}

remark[Comments on Assumptions (ref) and (ref)] \leavevmode \begin{enumerate}[label = (\roman*)] • Given Assumption (ref)(ii), there is no loss of generality by assuming $\Sigma \equiv \Omega \equiv {\mathbb{I}}$, the identity matrix of the same size as $\Sigma$ or $\Omega$. We make such a simplification throughout the paper, unless stated otherwise. • The assumptions on the nuisance parameters and their estimators in Assumption (ref) are quite mild. In particular, we do not assume $\widehat{a}$, $\widehat{b}$ converge to $a$, $b$ at any algebraic rate in $L_{2}$-norm. In fact, if content with $\sqrt{n}$-consistency instead of semiparametric efficiency, $\Vert \widehat{a} - a \Vert_{2} = o (1)$ and $\Vert \widehat{b} - b \Vert_{2} = o (1)$ can be even relaxed to $\Vert \widehat{a} - a \Vert_{2} = O (1)$ and $\Vert \widehat{b} - b \Vert_{2} = O (1)$; see liu2017semiparametric. • Assumption (ref) on the dictionary $\bar{{\mathsf{z}}}_{k}$ also appeared in robins2017minimax, liu2017semiparametric, liu2020nearly, liu2021assumption; also see comments in liu2020rejoinder. The $L_{\infty}$-stability condition (iii) have been established for Cohen-Daubechies-Vial wavelets, B-splines, and local polynomial partition series belloni2015some. It is possible to relax such a condition to a high-probability version, which we decide not to further pursue in this paper. \end{enumerate}

The following result on the sSOIF estimator $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ is a special case of Theorem (ref) to be revealed in Section (ref).

proposition[Bias and variance bounds of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$.] Under Assumptions (ref) -- (ref), with $k = o (n)$, one has the following: \begin{enumerate}[label = {\normalfont(\roman*)}] • The kernel estimation bias of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ satisfies \begin{equation} \begin{split} & \mathsf{kern-bias}_{2, k} (\widehat{\psi}_{1}) \coloneqq {\mathbb{E}} \left[ \widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}) \right] - \mathsf{bias}_{\theta, k} (\widehat{\psi}_{1}) \equiv {\mathbb{E}} \left[ \widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}) - \widehat{\mathbb{IF}}_{2, 2, k} \right] \\ & \lesssim \frac{k}{n} \left\{ \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{2} \Vert \widehat{b} - b \Vert_{2} + \left\Vert \frac{\widehat{a} - a}{a} \right\Vert_{2} \Vert \widehat{b} - b \Vert_{2} + \left( \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{2} \left\Vert \widehat{b} - b \right\Vert_{\infty} \wedge \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{\infty} \left\Vert \widehat{b} - b \right\Vert_{2} \right) \right\}. \end{split} \end{equation} • The variance of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ satisfies \begin{equation} \begin{split} \mathsf{var} \left[ \widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}) \right] \lesssim \frac{1}{n} \left\{ \frac{k}{n} + \left( \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{2} \left\Vert \widehat{b} - b \right\Vert_{\infty} \wedge \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{\infty} \left\Vert \widehat{b} - b \right\Vert_{2} \right) \right\}. \end{split} \end{equation} \end{enumerate}

remarkThe dependence on the condition number $k / n$ in the kernel estimation bias upper bound of the eSOIF estimator $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}_{\mathrm{nuis}})$ in liu2017semiparametric ($\sqrt{k \log k / n}$) is worse than that of the sSOIF estimator $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ reported here ($k / n$).

\allowdisplaybreaks

Proof sketch of Proposition (ref)

Kernel estimation bias bound

$\mathsf{kern\mbox{-}bias}_{2, k} (\widehat{\psi}_{1})$ can be controlled by repeatedly using the matrix identity $(A - B)^{-1} - A^{-1} = - A^{-1} B (A - B)^{-1}$ with $A \equiv \widehat{\Sigma}^{\dag} \coloneqq n^{-1} \sum_{i = 3}^{n} Q_{i}$ and $B = n^{-1} Q_{1, 2}$:

align*[align* omitted — 863 chars of source]

By choosing $J \asymp \log n$, the second term of the above display can be shown to be $o (n^{- 1 / 2})$.

For the first term, we only look at $j = 1, 2$ in the main text and the remaining analysis is a special case of the proof of Theorem (ref) in Appendix (ref).

For $j = 1$, we have \allowdisplaybreaks

align[align omitted — 1,682 chars of source]

where the last line follows from triangle inequality, Cauchy-Schwarz inequality and Assumptions (ref), (ref)(i) and (ref)(ii).

For $j = 2$, we have

align[align omitted — 1,175 chars of source]

Since $(\mathrm{II})$ is dominated by the term for $j = 1$, we only need to further analyze $(\mathrm{I})$ and $(\mathrm{III})$. $\mathrm{(III)}$ can be bounded by

equation[equation omitted — 528 chars of source]

where the first two terms are due to the first three terms in the (non-commutative) expansion of

equation[equation omitted — 120 chars of source]

and the third term comes from the last term in the above expansion. The appearance of the estimation error in $L_{\infty}$-norm is due to the opposite order of sample points indexed by $1$ and $2$ between the “meat” $Q_{2} Q_{1}$ and the “bread slices” of the “sandwich” structure $\bar{{\mathsf{z}}}_{k} (X_{1})^{\top} [\cdots] \bar{{\mathsf{z}}}_{k} (X_{2})$.

For $(\mathrm{I})$, we need to expand $({\mathbb{I}} - \widehat{\Sigma}^{\dag})^{2}$.

align[align omitted — 1,820 chars of source]

It is straightforward to see the first term in the last equality of the above display is dominated by the term for $j = 1$, whereas the second term can be shown to be bounded by

align*[align* omitted — 123 chars of source]

Taken together, the terms for $j = 1$ and $j = 2$ give the desired bound for $\mathsf{kern\mbox{-}bias}_{2, k} (\widehat{\psi}_{1})$ in (ref). It remains to prove the terms for $j \geq 3$ are of smaller order, which is deferred to Appendix (ref) for the general case. For $j \geq 3$, the corresponding term is of order

align*[align* omitted — 146 chars of source]

remarkNow is a perfect time to compare how the analysis of the kernel estimation bias of sSOIF differs from that of eSOIF of liu2017semiparametric. The only difference between the eSOIF and sSOIF estimators are the samples used to estimate $\Omega$. Using the nuisance sample instead, the (conditional) kernel estimation bias of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}_{\mathrm{nuis}})$ conditioning on the nuisance sample data is \begin{align*} & \ {\mathbb{E}} \left[ \widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}_{\mathrm{nuis}}) - \widehat{\mathbb{IF}}_{2, 2, k} \right] \\ = & \ {\mathbb{E}} \left[ (A_{1} \widehat{a} (X_{1}) - 1) \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \right] \left( \widehat{\Omega}_{\mathrm{nuis}} - {\mathbb{I}} \right) {\mathbb{E}} \left[ A_{2} \bar{{\mathsf{z}}}_{k} (X_{2}) (Y_{2} - \widehat{b} (X_{2})) \right]. \end{align*} From this, we can conclude \begin{align*} \left\vert {\mathbb{E}} \left[ \widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}_{\mathrm{nuis}}) - \widehat{\mathbb{IF}}_{2, 2, k} \right] \right\vert \lesssim \left( \frac{k \log k}{n} \right)^{1 / 2} \left\Vert \frac{\widehat{a}}{a} - 1 \right\Vert_{2} \Vert \widehat{b} - b \Vert_{2} \end{align*} by using matrix Bernstein or Khintchine inequality rudelson1999random, bandeira2021matrix; also see couillet2022random. liu2017semiparametric further show that \begin{align*} \left\vert {\mathbb{E}} \left[ \widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}_{\mathrm{nuis}}) - \widehat{\mathbb{IF}}_{2, 2, k} \right] \right\vert \lesssim \left( \frac{k \log k}{n} \right)^{m / 2} \left\Vert \frac{\widehat{a}}{a} - 1 \right\Vert_{2} \Vert \widehat{b} - b \Vert_{2}. \end{align*} However, as pointed out in Section (ref), the finite-sample performance of eHOIF estimators is not well-reflected by these upper bounds, prompting the need of developing sHOIF estimators.

Variance bound

The variance bound is technically involved. The missing steps can be found in Appendix (ref). The key step is to show

equation[equation omitted — 442 chars of source]

To prove (ref), it is sufficient to exhibit

equation[equation omitted — 629 chars of source]

and

equation[equation omitted — 724 chars of source]

Recall $\widehat{\Omega} = \widehat{\Sigma}^{-1}$ and $\widehat{\Sigma} = n^{-1} \sum_{i = 1}^{n} Q_{i}$. We introduce independent “ghost copies” $Q_{1}', Q_{2}'$ of $Q_{1}, Q_{2}$ and denote $\widehat{\Omega}' = (\widehat{\Sigma}')^{-1}$ and $\widehat{\Sigma}' = n^{-1} \sum_{i = 3}^{n} Q_{i} + n^{-1} (Q_{1}' + Q_{2}')$. Then (ref) is equivalent to

equation[equation omitted — 323 chars of source]

Let $\bar{\Omega} \equiv \bar{\Sigma}$ and $\bar{\Sigma} \equiv n^{-1} \sum_{i = 5}^{n} Q_{i}$. Repeating the matrix identity $(A + B)^{-1} - A^{-1} = - A^{-1} B (A + B)^{-1}$ on (ref) by setting $A = \bar{\Sigma}$ and $B = n^{-1} (Q_{1, 2} + Q_{3, 4})$ or $B = n^{-1} (Q_{1, 2}' + Q_{3, 4})$, we have

align*[align* omitted — 1,005 chars of source]

Let $J \asymp \log n$. The second term of the above display can be shown to be $o (1 / n)$. Proceeding to the first term, it is easy to see from Assumption (ref)(iii) that the term corresponding to $j = 1$:

align*[align* omitted — 1,341 chars of source]

Similarly, the terms corresponding to $j \geq 2$ can be shown to be $O \left( \frac{1}{n} \left( \frac{2 k}{n} \right)^{j - 1} \right)$, a consequence of Lemma (ref) below. Note that the extra factor $2$ appears because there are $O (2^{j})$ terms in total by expanding out $\left( \frac{Q_{1, 2} + Q_{3, 4}}{n} \right)^{j}$ and $\left( \frac{Q_{1, 2}' + Q_{3, 4}}{n} \right)^{j}$.

lemmaGiven a positive integer $j$. Given any pair of integers $j_{1} \geq 0, j_{2} > 0$ such that $j_{1} + j_{2} = j$, for any subset $\mathfrak{j} \subseteq \{0, 1\}^{j}$ of the $j$-dimensional Boolean hypercube with $\Vert \mathfrak{j} \Vert_{1} = j_{1}$, we have \begin{equation} {\mathbb{E}} \left[ \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \bar{{\mathsf{z}}}_{k} (X_{2}) \bar{{\mathsf{z}}}_{k} (X_{3})^{\top} \left( \prod_{\ell = 1}^{j} Q_{1, 2}^{\mathfrak{j}_{\ell}} Q_{3, 4}^{(1 - \mathfrak{j}_{\ell})} \right) \bar{{\mathsf{z}}}_{k} (X_{4}) \right] \lesssim k^{j - 1}. \end{equation} However, if $j_{1} = 0$ and $j_{2} = j$, we have \begin{equation} {\mathbb{E}} \left[ \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \bar{{\mathsf{z}}}_{k} (X_{2}) \bar{{\mathsf{z}}}_{k} (X_{3})^{\top} Q_{3, 4}^{j} \bar{{\mathsf{z}}}_{k} (X_{4}) \right] \lesssim k^{j}. \end{equation}

The proof of Lemma (ref) can be found in Appendix (ref). Finally, we defer the proof of (ref) to the online supplements, which can be proved in a similar fashion.

Numerical stability and time complexity of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$

As discussed in Section (ref), the key motivation for proposing sHOIF estimators is the numerical instability observed for eHOIF estimators. As a warm-up, we rigorously prove the numerical stability and calculate the time complexity of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ in this section. The reason why $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ can be numerically unstable is that when $k$ is near $n$, it is highly likely $\lambda_{\min} (\widehat{\Sigma}) \approx 0$ and hence $\lambda_{\max} (\widehat{\Omega})$ is close to infinity. But:

proposition$\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ does not depend on the eigenvalues of $\widehat{\Omega}$.

For ease of exposition, in what follows we let

itemize• $\bar{{\mathsf{Z}}}_{n, k} \coloneqq \left( \bar{{\mathsf{z}}}_{k} (X_{1}), \cdots, \bar{{\mathsf{z}}}_{k} (X_{n}) \right)^{\top}$ as the $n \times k$-matrix of the dictionary vectors for all $n$ samples; • $\bar{{\mathsf{Z}}}_{n, k}^{A} \coloneqq \left( A_{1} \bar{{\mathsf{z}}}_{k} (X_{1}), \cdots, A_{n} \bar{{\mathsf{z}}}_{k} (X_{n}) \right)^{\top}$ as the $n \times k$-matrix of the $A$-weighted dictionary vectors for all $n$ samples; • $\bm{{\mathcal{E}}}_{n, a} (\widehat{a}) \coloneqq \left( {\mathcal{E}}_{a} (\widehat{a}, O_{1}), \cdots, {\mathcal{E}}_{a} (\widehat{a}, O_{n}) \right)^{\top}$ and $\bm{{\mathcal{E}}}_{n, b} (\widehat{b}) \coloneqq \left( {\mathcal{E}}_{b} (\widehat{b}, O_{1}), \cdots, {\mathcal{E}}_{b} (\widehat{b}, O_{n}) \right)^{\top}$.

proofWe can rewrite $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ as \begin{equation} \widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}) = \frac{1}{n - 1} \bm{{\mathcal{E}}}_{n, a} (\widehat{a})^{\top} \left[ {\mathbb{I}} - \mathsf{Diag} \right] \left\{ \bar{{\mathsf{Z}}}_{n, k} \left( \bar{{\mathsf{Z}}}_{n, k}^{\top} \bar{{\mathsf{Z}}}_{n, k}^{A} \right)^{-1} \bar{{\mathsf{Z}}}_{n, k}^{A \top} \right\} \bm{{\mathcal{E}}}_{n, b} (\widehat{b}). \end{equation} Now apply Singular Value Decomposition (SVD) on the matrices $\bar{{\mathsf{Z}}}_{n, k}^{A}$: \begin{align*} \bar{{\mathsf{Z}}}_{n, k}^{A} = U^{A} \mathsf{Diag} (D^{A}) V^{A \top}. \end{align*} Then \begin{align*} \bar{{\mathsf{Z}}}_{n, k} \left( \bar{{\mathsf{Z}}}_{n, k}^{\top} \bar{{\mathsf{Z}}}_{n, k}^{A} \right)^{-1} \bar{{\mathsf{Z}}}_{n, k}^{A \top} = \bar{{\mathsf{Z}}}_{n, k}^{A} \left( \bar{{\mathsf{Z}}}_{n, k}^{A \top} \bar{{\mathsf{Z}}}_{n, k}^{A} \right)^{-1} \bar{{\mathsf{Z}}}_{n, k}^{A \top} = U^{A} U^{A \top}. \end{align*} So \begin{align*} (ref) = \frac{1}{n - 1} \bm{{\mathcal{E}}}_{n, a} (\widehat{a})^{\top} \left[ {\mathbb{I}} - \mathsf{Diag} \right] \left\{ U^{A} U^{A \top} \right\} \bm{{\mathcal{E}}}_{n, b} (\widehat{b}), \end{align*} which is completely independent of the eigenvalues of $\widehat{\Omega}$ ($(D^{A})^{2}$ up to constant).

Hence it is not surprising that $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ is numerically stable even when $k \rightarrow n$.

remarkIn a sense, $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ can be viewed as a {\it self-normalized} version of $\widehat{\mathbb{IF}}_{2, 2, k}$. It is generally expected that self-normalized statistics could have better statistical properties than the non-self-normalized ones pena2008self. However, whether the perspective of self-normalization is useful for establishing statistical properties of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ is still unclear to us and is worth pursuing as a research problem.

Furthermore, not only does the alternative formula (ref) of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ directly imply its numerical stability, but also it hints at the complexity of computing $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$. Barring the time complexity of SVD ($O (n k^{2})$), the time complexity of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ scales with $k n$ at a linear instead of a quadratic rate. This can be seen from (ref), in which only two vector-matrix products are involved, each taking $O (n k)$ operations. Thus we have

propositionThe time complexity of computing $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ is $O (n k^{2})$, dominated by that of SVD.

remarkAnother alternative way of arriving at the above conclusion is to observe that the $U$-statistic kernel of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$, denoted as $\widehat{\mathsf{IF}}_{2, 2, k, \bar{i}_{2}} (\widehat{\Omega})$, is separable, in the following sense: there exists a pair (but not necessarily a unique pair) of functions $h_{1}, h_{2}$ such that \begin{align*} \widehat{\mathsf{IF}}_{2, 2, k, \bar{i}_{2}} (\widehat{\Omega}) \equiv h_{1} (O_{i_{1}}, \widehat{\Omega}) \cdot h_{2} (O_{i_{2}}, \widehat{\Omega}). \end{align*}

The hierarchy of sHOIF estimators

As indicated in Section (ref), the $m$-th order sHOIF estimator $\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega})$ takes the same form as the $m$-th order eHOIF estimator $\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega}_{\mathrm{nuis}})$, with the sole difference that $\widehat{\Omega}_{\mathrm{nuis}}$ is replaced by $\widehat{\Omega}$. Formally, the $m$-th order sHOIF and the corresponding $m$-th order estimator of $\psi (\theta)$ read as follows:

equation[equation omitted — 717 chars of source]

remarkThe above sHOIF statistics are the same as the eHOIF statistics except that $\widehat{\Omega}$ is constructed from the estimation sample instead of the nuisance sample.

In this section, we first explain heuristically why $\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega})$ is enough to correct for the kernel estimation bias (see Section (ref)), after which the statistical, numerical and computational properties of $\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega})$ are stated formally.

Heuristic explanation

In what follows we explain heuristically why the kernel estimation bias $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ can be further corrected by adding:

align*[align* omitted — 358 chars of source]

and

align*[align* omitted — 402 chars of source]

For short, we define $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (3, 3), k} (\widehat{\Omega}) \coloneqq \widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega}) + \widehat{\mathbb{IF}}_{3, 3, k} (\widehat{\Omega})$ and $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (4, 4), k} (\widehat{\Omega}) \coloneqq \sum_{j = 2}^{4} \widehat{\mathbb{IF}}_{j, j, k} (\widehat{\Omega})$. Also note that the majority of this section is written for mathematical rigor.

Simple algebra gives

align*[align* omitted — 1,080 chars of source]

and \allowdisplaybreaks

align*[align* omitted — 1,866 chars of source]

First, observe that the expectation of the oracle version of $\widetilde{\widehat{\mathbb{IF}}}_{3, 3, k} (\widehat{\Omega})$

align*[align* omitted — 295 chars of source]

exactly cancels (ref), the leading-order part of the kernel estimation bias of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$ corresponding to $j = 1$.

Next, observe that the expectation of the oracle version of $\widetilde{\widehat{\mathbb{IF}}}_{4, 4, k} (\widehat{\Omega})$ is

align*[align* omitted — 562 chars of source]

which again cancels the kernel estimation bias of $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (3, 3), k} (\widehat{\Omega})$ truncated at $j = 2$, dominated by

equation[equation omitted — 499 chars of source]

which can be derived from (ref), (ref) and the kernel estimation bias of $\widehat{\mathbb{IF}}_{3, 3, k} (\widehat{\Omega})$ truncated at level $j = 1$; see Appendix (ref) for a more detailed calculation. Hence $\widehat{\mathbb{IF}}_{(3, 3) \rightarrow (4, 4), k} (\widehat{\Omega})$ further reduces the kernel estimation bias of $\widehat{\mathbb{IF}}_{2, 2, k} (\widehat{\Omega})$.

Characterization of the bias and variance of the sHOIF estimators

We now state the main theoretical result of this paper.

theoremUnder Assumptions (ref) -- (ref), with $k \lesssim \frac{n}{\log^{2} n}$ and $m \gtrsim \sqrt{\log n}$, one has the following: \begin{enumerate}[label = {\normalfont(\roman*)}] • The kernel estimation bias of $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega})$ satisfies \begin{equation} \begin{split} & \mathsf{kern-bias}_{m, k} (\widehat{\psi}_{1}) \coloneqq {\mathbb{E}} \left[ \widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}) \right] - \mathsf{bias}_{\theta, k} (\widehat{\psi}_{1}) \equiv {\mathbb{E}} \left[ \widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}) - \widehat{\mathbb{IF}}_{2, 2, k} \right] \\ & \lesssim \left( \frac{k m}{n} \right)^{\lceil \frac{\lceil \frac{m - 1}{2} \rceil - 1}{2} \rceil \vee 1} \left\{ \begin{array}{c} \left\Vert \dfrac{\widehat{a} - 1}{a} \right\Vert_{2} \Vert \widehat{b} - b \Vert_{2} + \left\Vert \dfrac{\widehat{a} - a}{a} \right\Vert_{2} \Vert \widehat{b} - b \Vert_{2} \\ + \left( \left\Vert \dfrac{\widehat{a} - 1}{a} \right\Vert_{2} \left\Vert \widehat{b} - b \right\Vert_{\infty} \wedge \left\Vert \dfrac{\widehat{a} - 1}{a} \right\Vert_{\infty} \left\Vert \widehat{b} - b \right\Vert_{2} \right) \end{array} \right\}. \end{split} \end{equation} • For $m \geq 2$, the variance of $\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega})$ satisfies \begin{equation} \begin{split} \mathsf{var} \left[ \widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega}) \right] \lesssim \frac{1}{n} \left\{ \frac{k}{n} + \left( \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{2} \left\Vert \widehat{b} - b \right\Vert_{\infty} \wedge \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{2} \left\Vert \widehat{b} - b \right\Vert_{\infty} \right) \right\}. \end{split} \end{equation} And thus the variance of $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega})$ satisfies \begin{equation} \begin{split} \mathsf{var} \left[ \widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}) \right] \lesssim \frac{1}{n} \left\{ \frac{k}{n} + \left( \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{2} \left\Vert \widehat{b} - b \right\Vert_{\infty} \wedge \left\Vert \frac{\widehat{a} - 1}{a} \right\Vert_{2} \left\Vert \widehat{b} - b \right\Vert_{\infty} \right) \right\}. \end{split} \end{equation} \end{enumerate}

The proof of the above theorem can be found in Appendix (ref) (for kernel estimation bias bound) and the online supplements (for variance bound).

remark[Asymptotic normality and the bootstrap approximation] As shown in liu2020nearly, the asymptotic normality of the oracle statistic $\frac{\widehat{\mathbb{IF}}_{2, 2, k} - \mathsf{bias}_{\theta, k} (\widehat{\psi}_{1})}{\mathsf{se}_{\theta} (\widehat{\mathbb{IF}}_{2, 2, k})}$ follows from Theorem 1 of bhattacharya1992class whence $1 \ll k \ll n^{2}$. Thus to show CLT of $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega})$ for any $m \geq 2$, it is sufficient to demonstrate under what conditions $\mathsf{kern\mbox{-}bias}_{m, k} (\widehat{\Omega}) \ll \mathsf{se}_{\theta} (\widehat{\mathbb{IF}}_{2, 2, k})$. Bootstrap approximation (and its rate) of the distribution of $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega})$, or even of $\widehat{\mathbb{IF}}_{2, 2, k}$ is still an important open problem, though liu2021assumption have made some partial progress. A more thorough study of the conditions under which central limit theorem (CLT) or bootstrap approximation holds is beyond the scope of this paper.

Numerical stability and time complexity of sHOIF estimators

In what follows we consider the numerical and computational properties of sHOIF estimators, which extends the results in Section (ref) to higher-order. The first result in this section, Theorem (ref), earmarks the “stability” of sHOIF estimators in terms of their independence of the eigenvalues of $\widehat{\Omega}$, the root cause of the instability of eHOIF estimators.

theorem$\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega})$ does not depend on the eigenvalues of $\widehat{\Omega}$.

proofThe proof resembles the proof of Proposition (ref) closely by realizing that, for any $i, j \in [n]$, \begin{align*} A_{i} \bar{{\mathsf{z}}}_{k} (X_{i})^{\top} \widehat{\Omega} \bar{{\mathsf{z}}}_{k} (X_{j}) = U^{A}_{i, \bullet} U^{A \top} U U_{j, \bullet}, \end{align*} which is completely independent of the eigenvalues of $\widehat{\Omega}$.

Hence sHOIF estimators do not suffer from any numerical instability resulted from the large condition number of the sample Gram matrix when we let $k$ near $n$ in practice.

remarkTheorem (ref) also suggests a better way to compute sHOIF estimators. Instead of computing the sample Gram matrix $\widehat{\Sigma}$ and its inverse $\widehat{\Omega}$ using numerical methods, we should instead perform SVD on the basis matrices $\bar{{\mathsf{Z}}}_{n, k}$ and $\bar{{\mathsf{Z}}}_{n, k}^{A}$ and then compute $\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega})$. In fact, the upcoming R package wanis2023machine for computing HOIF related statistics exactly uses this strategy.

Since sHOIF estimators are numerically stable and thus are potentially useful tools for statistical practice liu2020nearly, wanis2023machine, it is worth discussing the computational complexity of sHOIF estimators for general order $m$ as well.

theoremThe time complexity of computing $\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega})$ is $O (\max \{(n k)^{\lceil (m - 1) / 2 \rceil}, n k^{2}\})$.

proofSimilar to the proof of Proposition (ref), we need to rewrite $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega})$ in the form of a linear combination of $V$-statistics. Without loss of generality, we take ${\mathcal{E}}_{a} (\widehat{a}; O) \equiv {\mathcal{E}}_{b} (\widehat{b}; O) \equiv 1$. But let us first represent $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega})$ as the following series: \begin{align*} \widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}) \equiv \sum_{j = 1}^{m} (-1)^{j} \binom{m - 1}{j - 1} {\mathbb{U}}_{n, j} \left[ \bar{{\mathsf{z}}}_{k} (X_{1})^{\top} \widehat{\Omega} \cdot \prod_{s = 3}^{m} \left( Q_{s} \widehat{\Omega} \right) \cdot \bar{{\mathsf{z}}}_{k} (X_{2}) \right]. \end{align*} Note that the number of summations in $\widehat{\mathbb{IF}}_{m, m, k} (\widehat{\Omega})$ is \begin{align*} n (n - 1) \cdots (n - m + 1) = \sum_{j = 1}^{m} (-1)^{m - j} s (m, j) n^{j} \end{align*} where $s (m, j)$ are unsigned Stirling numbers of the first kind, or the number of permutations on $m$ elements with $j$ cycles. Accordingly one can write an $m$-th order $U$-statistic into a linear combination of $V$-statistics from order $1$ to order $m$, with the number of $j$-th order $V$-statistics, for $j = 1, \cdots, m$, equal to $s (m, j)$. The proof is completed by leveraging the special structure of the $U$-statistic kernel for sHOIF estimators.

remarkConsidering Theorem (ref) and Theorem (ref) in tandem, there is a clear statistical-computational trade-off. However, whether or not such statistical-computational trade-off is an emanation of possibly intrinsic computational hardness of estimating certain smooth statistical functionals is still an open problem Finally, we briefly comment on our philosophical stance on the usefulness of sHOIF estimators. sHOIF estimators are effectively infinite-order $U$-statistics, so given the current computing devices, there is no doubt that practitioners are not using sHOIF estimators in practice in near term. This is “conditional” on the availability of hardware. The numerical stability or lack thereof, however, is an issue regardless of the availability of more powerful computing resources.

remarkTheorem (ref) also applies to eHOIF estimators liu2017semiparametric and the original HOIF estimators of robins2016technical, that needs an estimate of the density of the covariates $X$, if the time for density estimation is not counted.

Applications of the statistical properties of sHOIF estimators

Semiparametric efficiency under minimal H\"{o}lder assumptions on the nuisance functions

In nonparametric statistics, the optimality of a statistical procedure is often evaluated under the H\"{o}lder nuisance models.

The above calculations culminate into the following theorem, which is the second main result of this paper.

theoremIf ${\mathcal{A}} \times {\mathcal{B}} \subseteq {\mathcal{H}} (s_{a}, {\mathcal{X}}) \times {\mathcal{H}} (s_{b}, {\mathcal{X}})$ with $(s_{a} + s_{b}) / 2 \geq d / 4$, and choosing $m \asymp \sqrt{\log n}$ and $k \lesssim n / \log (n)^{2}$, \begin{equation} \sqrt{n} \left( \widehat{\psi}_{m, k} - \psi (\theta) \right) \overset{{\mathcal{L}}}{\rightarrow} {\mathcal{N}} (0, {\mathbb{E}} [\mathsf{IF}_{1} (\theta)^{2}]) \end{equation} where ${\mathbb{E}} [\mathsf{IF}_{1} (\theta)^{2}]$ is the semiparametric efficiency bound of $\psi (\theta)$.

remarkAccording to the lower bound of robins2009semiparametric under the H\"{o}lder nuisance model, $(s_{a} + s_{b}) / 2 \geq d / 4$ is the minimal condition for the existence of a semiparametric efficient estimator of $\psi (\theta)$. It is not unreasonable to expect that this minimal condition also holds for most, if not all, of the DRFs.

Implications on the assumption-free bias testing procedure of liu2020nearly and liu2021assumption

In light of the growing interest in understanding the performance of deep-learning-based causal inference farrell2021deep, chen2020causal and the gap between these theoretical results and empirical performance xu2022deepmed, liu2020nearly proposed the following oracle assumption-free valid nominal $\alpha$-level test statistic:

equation[equation omitted — 325 chars of source]

for the following null hypothesis:

equation[equation omitted — 171 chars of source]

where

equation[equation omitted — 248 chars of source]

liu2021assumption in turn constructed a feasible assumption-lean valid nominal $\alpha$-level test statistic

equation[equation omitted — 451 chars of source]

and the following higher-order test statistic based on eHOIF estimators:

equation[equation omitted — 448 chars of source]

liu2021assumption showed that all the standard errors in the above test statistics can be estimated consistently by certain bootstrapping procedure. More importantly, they proved the following.

propositionLet \begin{align*} \mathsf{cs-bias}_{k} (\widehat{\psi}_{1}) \coloneqq \left\{ {\mathbb{E}} \left[ \Pi [(A \widehat{a} - 1) | \bar{{\mathsf{z}}}_{k}] (X)^{2} \right] {\mathbb{E}} \left[ \Pi [A (\widehat{b} - b) | \bar{{\mathsf{z}}}_{k}] (X)^{2} \right] \right\}^{1 / 2}. \end{align*} Under the assumptions of Theorem (ref), with $k \lesssim n / (\log n)^{2}$ and the following extra condition: \begin{equation} \vert \mathsf{kern-bias}_{3, k} (\widehat{\psi}_{1}) \vert \not\ll \mathsf{cs-bias}_{k} (\widehat{\psi}_{1}) \end{equation} then $\widehat{\chi}_{3, k} (\widehat{\Omega}_{tr})$ is a valid nominal $\alpha$-level test of ${\mathsf{H}}_{0} (\delta)$ (ref). The extra condition (ref) can be relaxed to \begin{equation} \vert \mathsf{kern-bias}_{m, k} (\widehat{\psi}_{1}) \vert \not\ll \mathsf{cs-bias}_{k} (\widehat{\psi}_{1}) \left( \frac{k \log k}{n} \right)^{\frac{m - 1}{2}} \end{equation} if one uses $\widehat{\chi}_{m, k} (\widehat{\Omega}_{tr})$ instead of $\widehat{\chi}_{3, k} (\widehat{\Omega}_{tr})$.

We can similarly define the following sHOIF-based test statistics: for $m \geq 2$,

align*[align* omitted — 402 chars of source]

Then as an immediate corollary of Theorem (ref), we have

theoremUnder the assumptions of Theorem (ref), with $k \lesssim n / (\log n)^{2}$ and a different relaxed extra condition from (ref): \begin{equation} \vert \mathsf{kern-bias}_{m, k} (\widehat{\psi}_{1}) \vert \not\ll \mathsf{cs-bias}_{k} (\widehat{\psi}_{1}) \left( \frac{k}{n} \right)^{m - 1} \end{equation} then $\widehat{\chi}_{m, k} (\widehat{\Omega})$ is a valid nominal $\alpha$-level test of ${\mathsf{H}}_{0} (\delta)$ (ref).

Given the above theoretical guarantees, and further considering that the sHOIF estimators and tests have better finite-sample performance than the corresponding eHOIF estimators and tests, we recommend using $\widehat{\chi}_{m, k} (\widehat{\Omega})$ in practice. For more examples of its application, see wanis2023machine.

Further extensions of sHOIF estimators

Generalization to the entire class of DRFs

In this subsection, we briefly comment on how our results can be generalized to the entire class of DRFs characterized in rotnitzky2021characterization. The class of DRFs includes many other functionals that arise in substantive studies in (bio)statistics, epidemiology, economics, and social sciences, including:

itemize• the expected conditional variance, which is useful for constructing confidence/predictive sets robins2006adaptive; • the expected conditional covariance, which is useful for both causal inference and conditional independence testing shah2020hardness; • average causal effect of continuous treatment, which is important for treatment allocations ai2021unified, bonvini2022fast.

rotnitzky2021characterization defined the class of DRFs as follows:

definition[Doubly Robust Functionals; Definition 1 of rotnitzky2021characterization] $\psi (\theta)$ is a doubly robust functional if, for each $\theta \in \Theta$ there exists $a: x \mapsto a (x) \in {\mathcal{A}}$ and $b: x \mapsto b (x) \in {\mathcal{B}}$ such that (i) $\theta = (a, b, \theta \setminus \{b, p\})$ and $\Theta = {\mathcal{A}} \times {\mathcal{B}} \times \Theta \setminus \{{\mathcal{A}}, {\mathcal{B}}\}$ and (ii) for any $\theta, \theta' \in \Theta$ \begin{equation} \psi (\theta) - \psi (\theta^{\prime}) + {\mathbb{E}} \left[ \mathsf{IF}_{1} (\theta^{\prime}) \right] = {\mathbb{E}} \left[ S (a (X) - a^{\prime} (X)) (b (X) - b^{\prime} (X)) \right] \end{equation} where $S \equiv s (O)$ with $o \mapsto s (o)$ a known function that does not depend on $\theta$ or $\theta'$ satisfying either ${\mathbb{P}}_{\theta} (S \geq 0) = 1$ or ${\mathbb{P}}_{\theta} (S \leq 0) = 1$. We also denote $\lambda (x) \coloneqq {\mathbb{E}} [S | X = x]$. Then the first-order influence function of $\psi (\theta)$ has the following form: given $\theta' \equiv (a', b', \theta' \setminus \{a', b'\})^{\top} \in \Theta$, \begin{equation} \mathsf{IF}_{1} (\theta') \equiv S a' (X) b' (X) + m_{a} (O, a') + m_{b} (O, b') + S_{0} \end{equation} where $S_{0}$ is some known statistic that does not depend on $a'$ and $b'$, and $h \mapsto m_{a} (o, h)$ for $h \in {\mathcal{A}}$ and $h \mapsto m_{b} (o, h)$ for $h \in {\mathcal{B}}$ are two known linear maps satisfying \begin{align*} & {\mathbb{E}} \left[ S h (X) b (X) + m_{a} (O, h) \right] \equiv 0, \ \forall \ h \in {\mathcal{A}}, \\ & {\mathbb{E}} \left[ S a (X) h (X) + m_{b} (O, h) \right] \equiv 0, \ \forall \ h \in {\mathcal{B}}. \end{align*} As a result, $\psi (\theta) \equiv {\mathbb{E}} [m_{a} (O, a) + S_{0}] \equiv {\mathbb{E}} [m_{b} (O, b) + S_{0}]$.

remarkFor $\psi (\theta) \equiv {\mathbb{E}} [Y (1)]$ under strong ignorability, $S$, $a$, $b$, $m_{a} (O, a)$, $m_{b} (O, b)$, and $S_{0}$ correspond to $- A$, $\{{\mathbb{E}} [A | X]\}^{-1}$, ${\mathbb{E}} [Y | X, A = 1]$, $A Y a (X)$, $b (X)$, and $0$, respectively. Thus $\lambda (X) = {\mathbb{E}} [- A | X] = - \frac{1}{a (X)}$. We also have \begin{align*} & {\mathbb{E}} [{\mathcal{E}}_{a} (\widehat{a}; O) \bar{{\mathsf{z}}}_{k} (X)] = {\mathbb{E}} \left[ \frac{1}{a (X)} (\widehat{a} (X) - a (X)) \bar{{\mathsf{z}}}_{k} (X) \right], \\ & {\mathbb{E}} [{\mathcal{E}}_{b} (\widehat{b}; O) \bar{{\mathsf{z}}}_{k} (X)] = {\mathbb{E}} \left[ \frac{1}{a (X)} (\widehat{b} (X) - b (X)) \bar{{\mathsf{z}}}_{k} (X) \right]. \end{align*}

We have the following notation correspondence that maps the results for $\psi (\theta) \equiv {\mathbb{E}} [Y (1)]$ under strong ignorability to any DRF $\psi (\theta)$:

itemize• ${\mathcal{E}}_{a} (\widehat{a}; O) \bar{{\mathsf{z}}}_{k} (X) \Rightarrow {\mathcal{E}}_{a} (\widehat{a}, \bar{{\mathsf{z}}}_{k}; O)$ and ${\mathcal{E}}_{b} (\widehat{b}; O) \bar{{\mathsf{z}}}_{k} (X) \Rightarrow {\mathcal{E}}_{b} (\widehat{b}, \bar{{\mathsf{z}}}_{k}; O)$ where ${\mathcal{E}}_{a} (\widehat{a}, \bar{{\mathsf{z}}}_{k}; O)$ and ${\mathcal{E}}_{b} (\widehat{b}, \bar{{\mathsf{z}}}_{k}; O)$ satisfy \begin{align*} & {\mathbb{E}} [{\mathcal{E}}_{a} (\widehat{a}, \bar{{\mathsf{z}}}_{k}; O)] = {\mathbb{E}} \left[ \lambda (X) (\widehat{a} (X) - a (X)) \bar{{\mathsf{z}}}_{k} (X) \right], \\ & {\mathbb{E}} [{\mathcal{E}}_{b} (\widehat{b}, \bar{{\mathsf{z}}}_{k}; O)] = {\mathbb{E}} \left[ \lambda (X) (\widehat{b} (X) - b (X)) \bar{{\mathsf{z}}}_{k} (X) \right]. \end{align*} • $\Sigma = {\mathbb{E}} [A \bar{{\mathsf{z}}}_{k} (X) \bar{{\mathsf{z}}}_{k} (X)^{\top}]$ $\Rightarrow$ $\Sigma = {\mathbb{E}} [S \bar{{\mathsf{z}}}_{k} (X) \bar{{\mathsf{z}}}_{k} (X)^{\top}]$ and $\widehat{\Sigma} = {\mathbb{P}}_{n} [A \bar{{\mathsf{z}}}_{k} (X) \bar{{\mathsf{z}}}_{k} (X)^{\top}]$ $\Rightarrow$ $\widehat{\Sigma} = {\mathbb{P}}_{n} [S \bar{{\mathsf{z}}}_{k} (X) \bar{{\mathsf{z}}}_{k} (X)^{\top}]$.

With the above mappings, all the theoretical results for $\psi (\theta) \equiv {\mathbb{E}} [Y (1)]$ developed herein can be applied to those for an arbitrary DRF $\psi (\theta)$ {\it mutatis mutandis}.

A special case: the expected conditional covariance

Before concluding our paper, we further study the implications of the sHOIF theory developed so far for a special cases of DRFs: the expected conditional covariance between two random variables $A$ and $Y$ given a third random variable $X$, $\psi \coloneqq {\mathbb{E}} [\mathsf{cov} (A, Y | X)]$. When $A = Y$ almost surely, $\psi$ reduces to the expected conditional variance of $A$ given $X$, $\psi \coloneqq {\mathbb{E}} [\mathsf{cov} (A | X)]$. For differences between these two parameters, see an extended discussion in liu2020nearly.

The main feature that distinguishes $\psi$ from many other DRFs is $S = 1$, which leads to its SOIF:

align*[align* omitted — 284 chars of source]

in which $\bar{\Omega} \equiv \{{\mathbb{E}} [S \bar{{\mathsf{z}}}_{k} (X) \bar{{\mathsf{z}}}_{k} (X)^{\top}]\}^{-1} \equiv \{{\mathbb{E}} [\bar{{\mathsf{z}}}_{k} (X) \bar{{\mathsf{z}}}_{k} (X)^{\top}]\}^{-1}$ only depends on the distribution of $X$. Also, $\widehat{a}$ and $\widehat{b}$ are nuisance estimates of $a (x) = {\mathbb{E}} [A | X = x]$ and $b (x) = {\mathbb{E}} [Y | X = x]$ in this context. This leads to the following improved kernel estimation bias bound:

corollaryUnder Assumptions (ref) -- (ref), with $k \lesssim \frac{n}{\log^{2} n}$ and $m \gtrsim \sqrt{\log n}$, one has the following: The kernel estimation bias of $\widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega})$ satisfies \begin{equation} \begin{split} & \mathsf{kern-bias}_{m, k} (\widehat{\psi}_{1}) \coloneqq {\mathbb{E}} \left[ \widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}) \right] - \mathsf{bias}_{\theta, k} (\widehat{\psi}_{1}) \equiv {\mathbb{E}} \left[ \widehat{\mathbb{IF}}_{(2, 2) \rightarrow (m, m), k} (\widehat{\Omega}) - \widehat{\mathbb{IF}}_{2, 2, k} \right] \\ & \lesssim \left( \frac{k m}{n} \right)^{\lceil \frac{\lceil \frac{m - 1}{2} \rceil - 1}{2} \rceil \vee 1} \left\{ \left\Vert \widehat{a} - a \right\Vert_{2} \Vert \widehat{b} - b \Vert_{2} + \left( \left\Vert \widehat{a} - a \right\Vert_{2} \left\Vert \widehat{b} - b \right\Vert_{\infty} \wedge \left\Vert \widehat{a} - a \right\Vert_{\infty} \left\Vert \widehat{b} - b \right\Vert_{2} \right) \right\}. \end{split} \end{equation}

Note that the variance bound is improved in a similar manner and is omitted here.

Discussion

In this paper, we propose a novel class of HOIF estimators, stable HOIF (sHOIF) estimators, for the doubly robust functionals (DRFs) characterized in rotnitzky2021characterization. They are semiparametric efficient under the minimal H\"{o}lder-smoothness condition $\frac{s_{a} + s_{b}}{2} > \frac{d}{4}$ of robins2009semiparametric, allowing the dimension $k$ of the basis function diverging at a rate just slower than the sample size $n$. As can be seen from Theorem (ref), sHOIF estimators have improved rate of convergence than eHOIF developed in liu2017semiparametric. More importantly, as well documented in the simulation studies of liu2020nearly and wanis2023machine, the sHOIF estimators also have significantly better finite-sample performance over existing higher-order estimators in practice, making them more amenable for tasks such as testing if the bias of a first-order DML estimator $\widehat{\psi}_{1}$ of a causal effect is dominated by its standard error liu2020nearly, liu2021assumption, wanis2023machine. Finally, we end our paper by mentioning several future research directions:

enumerate[label = (\arabic*)] • It will be interesting to study if one can extend the idea of sHOIF estimators to the non-$\sqrt{n}$-estimable regimes by, for instance, estimating $\Omega$ via some shrinkage or regularized algorithms. As conjectured in robins2016technical, the minimax convergence rate of the functionals studied in this paper may depend on the regularity of the density of the covariates $X$. Hence it is expected that the shrinkage or regularization also depends on the density of $X$. Simulation studies in liu2020nearly suggest the nonlinear shrinkage covariance matrix estimators of ledoit2012nonlinear could be a viable option. Preliminary simulation studies in liu2020nearly and wanis2023machine suggest that the performance of these shrinkage covariance matrix estimators does degrade with the smoothness of the design density. • As pointed out in kennedy2022minimax, their Second-Order R-Learner (SORL) for CATE also involves inverting large Gram matrices of certain basis functions (in which they use the Legendre polynomials) under additional complexity-reducing assumptions on the covariates $X$. It will be interesting to investigate if the sHOIF estimators can be generalized to the CATE estimation problems and stabilize their SORL or even HORL estimators. • Another important open problem was also mentioned in van2014higher, liu2020nearly, liu2021assumption. To define HOIFs for DRFs, one needs to choose a set of $k$-dimensional basis functions $\bar{{\mathsf{z}}}_{k}$ or an approximation kernel $K$ of the Kronecker delta function, ideally in prior to the data analysis. However, such a strategy seems to go against the current data analytic paradigm, which strongly advocates learning representations (e.g. in the form of bases or kernels) adaptively from data rather than choosing some fixed bases/frames {\it a priori}. Prominent examples include DNNs, autoencoders, and GANs. It is thus interesting to construct HOIF estimators along different basis directions and then select one or aggregate all, guided by certain optimality criterion. We leave this important and difficult problem to future endeavor. • It will be interesting to also derive HOIFs and sHOIFs for identifiable causal effect functionals in graphical models with latent variables bhattacharya2022semiparametric and implicitly defined functionals robins2016technical, ai2021unified in general semiparametric regression problems for improved quality of estimation and statistical inference, which however requires extension of the current work to $U$-processes, a much more difficult research problem that we are working on in a separate paper.