EconBase
← Back to paper

Debiased Nonparametric Regression for Statistical Inference and Distributionally Robustness

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.

30,101 characters · 16 sections · 28 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.

Debiased Nonparametric Regression for Statistical Inference and Distributionally Robustness

abstractThis study proposes a debiasing method for smooth nonparametric estimators. While machine learning techniques such as random forests and neural networks have demonstrated strong predictive performance, their theoretical properties remain relatively underexplored. In particular, many modern algorithms lack guarantees of pointwise and uniform risk convergence, as well as asymptotic normality. These properties are essential for statistical inference and robust estimation and have been well-established for classical methods such as Nadaraya-Watson regression. To ensure these properties for various nonparametric regression estimators, we introduce a model-free debiasing method. By incorporating a correction term that estimates the conditional expected residual of the original estimator, or equivalently, its estimation error, into the initial nonparametric regression estimator, we obtain a debiased estimator that satisfies pointwise and uniform risk convergence, along with asymptotic normality, under mild smoothness conditions. These properties facilitate statistical inference and enhance robustness to covariate shift, making the method broadly applicable to a wide range of nonparametric regression problems.

Introduction

This study investigates the problem of nonparametric regression. While modern machine learning methods enable the estimation of complex regression functions, they often lack the theoretical guarantees that traditional estimators possess, such as pointwise and uniform mean squared error (MSE) convergence, as well as asymptotic normality. Asymptotic normality is crucial for statistical inference, whereas uniform risk convergence plays an important role in prediction under covariate shift. In this study, we propose a debiasing method for regression estimators and show that, under mild smoothness conditions, they achieve asymptotic normality and uniform convergence.

Here, we formulate the problem. Let $X \in {\mathcal{X}}$ be a $d$-dimensional covariate and $Y \in {\mathcal{Y}}$ be a target variable, where ${\mathcal{X}} \subseteq {\mathbb{R}}^d$ and ${\mathcal{Y}} \subset {\mathbb{R}}$ represent the covariate and target spaces, respectively. Let $P$ denote the joint distribution of $(X, Y)$, and define the regression function under $P$ as \[ f_P(X) = {\mathbb{E}}_P\big[Y \mid X\big], \] where ${\mathbb{E}}_P\big[\cdot\big]$ represents the expectation operator with respect to $P$. We assume access to observations $\big\{(X_i, Y_i)\big\}_{i=1}^n$, where $(X_i, Y_i)$ are independent and identically distributed samples from a true distribution $P_0$. The true regression function is denoted by $f_0 = f_{P_0}$. Our objective is to estimate $f_0$ from the observations $\big\{(X_i, Y_i)\big\}_{i=1}^n$, aiming for pointwise and uniform risk convergence, as well as asymptotic normality.

\paragraph{Notation.} For a real-valued vector $z$, let $\|z\|$ denote its Euclidean norm. For a measurable function $h\colon {\mathcal{X}} \to {\mathbb{R}}$, let $\|h\|_\infty \coloneqq \sup_{x \in {\mathcal{X}}} |h(x)|$ denote the sup-norm, and $\|h\|_2 \coloneqq \sqrt{{\mathbb{E}}\left[h(X)^2\right]}$ denote the $L^2$ norm.

Content of this study

In Section (ref), we define our debiased estimator for nonparametric regression. Our proposed estimator consists of the following three steps: (1) estimate the regression function $f_0$ using a nonparametric regression method; (2) estimate the expected conditional residual, or equivalently, the estimation error of the first-stage estimator, using local polynomial regression; and (3) add the second-stage estimator to the first-stage estimator.

In Section (ref), we present an example of our proposed debiased estimator, applying local polynomial regression for conditional expected residual estimation.

In Section (ref), we establish that our debiased nonparametric regression estimator achieves pointwise and uniform MSE convergence if either (a) the first-stage nonparametric regression estimator and the true regression function, or (b) their difference, belongs to the H\"older class. Using pointwise MSE convergence, we also prove asymptotic normality. Furthermore, our estimator possesses the doubly robust property: if either the first-stage nonparametric regression estimator or the second-stage conditional expected residual estimator is consistent, then the resulting debiased estimator is also consistent.

The key takeaway is that our debiased estimator ensures desirable theoretical properties, typically associated with classical methods, for a wide range of modern first-stage estimators. As long as certain smoothness conditions are met, our method effectively corrects bias in a broad class of initial estimators, equipping the resulting estimator with strong theoretical guarantees.

Related work

Doubly robust and debiased estimation has attracted significant attention across various fields, including statistics tsiatis2007semiparametric, economics ChernozhukovVictor2018Dmlf, epidemiology BangRobins2005, and machine learning Kallus2020double. With the advancement of machine learning regression models, debiasing methods based on doubly robust estimators have been extensively studied ChernozhukovVictor2018Dmlf, particularly in conjunction with sample splitting techniques klaassen1987,vanderVaart2002,ZhengWenjing2011CTME. While much of the existing literature focuses on parametric or semiparametric models, this study extends these approaches to nonparametric regression. The application of debiased machine learning to nonparametric regression is discussed in kennedy2023semiparametricdoublyrobusttargeted, though their focus differs from ours.

Nonparametric regression in machine learning presents several challenges:

itemize• Pointwise and uniform MSE convergence is not established for many estimators, such as neural networks and random forests. • Asymptotic normality remains unproven for numerous machine learning estimators.

Pointwise and uniform MSE convergence and asymptotic normality are critical in statistical analysis, motivating significant efforts to establish these properties. For instance, Wager2018 and Mourtada2020 demonstrate asymptotic normality and optimality, respectively, for their modified random forests, but not for the original random forest of Breiman2001.

Pointwise and uniform convergence is particularly relevant for robustness under distributional shifts, such as covariate shift Shimodaira2000. These convergence properties are rarely established for modern nonparametric estimators, whereas classical methods like local linear and series regression possess them. This gap arises in part due to the data-adaptive nature of modern methods and the frequent use of empirical process arguments, which often focus on population MSE, where expectations are taken over the covariate distribution rather than on pointwise or sup-norm MSE. SchmidtHieber2024 addresses this concern by demonstrating restricted uniform optimality for neural network regression in the one-dimensional covariate setting, but this approach does not generalize easily.

Our method shares some motivations with debiasing approaches in high-dimensional regression. For example, the Lasso estimator, while widely used, introduces a bias that decreases with sample size at a slower rate than $\sqrt{n}$, thereby hindering asymptotic normality tibshirani96regression,Buhlmann2011. To mitigate this, vandeGeer2014 proposes incorporating bias-correction terms, leading to improved asymptotic properties such as normality and efficiency Jankova2018.

This study generalizes our previous work on doubly robust methods for nonparametric regression discontinuity design kato2024doublyrobustregressiondiscontinuity. Concurrently with our work, chernozhukov2024conditionalinfluencefunctions and vanderlaan2025automaticdebiasedmachinelearning also explore double machine learning in the context of nonparametric regression, though their emphases differ substantially from ours. For example, chernozhukov2024conditionalinfluencefunctions develop conditional influence functions, which are related to the framework of Ichimura2022.

The debiased estimator

In this section, we define our debiased estimator in a general way. First, we randomly split the observations ${\mathcal{D}}$ into two datasets, ${\mathcal{D}}^{(1)}$ and ${\mathcal{D}}^{(2)}$, such that ${\mathcal{D}} = {\mathcal{D}}^{(1)} \cup {\mathcal{D}}^{(2)}$. For simplicity, let $n$ be even and set $m = n/2$. Then, for each $\ell \in \{1, 2\}$, we define \[ {\mathcal{D}}^{(\ell)} \coloneqq \big\{(X_i, Y_i)\big\}_{i \in {\mathcal{I}}^{(\ell)}}, \] where ${\mathcal{I}}^{(\ell)}$ is the index set of ${\mathcal{D}}^{(\ell)}$, i.e., ${\mathcal{I}}^{(1)} \cup {\mathcal{I}}^{(2)} = \{1, 2, \dots, n\}$.

We consider the following three-stage estimation for each point $x_0 \in {\mathcal{X}}$ of interest:

description• Estimate $f_0$ using any smooth model, and denote this estimator by $\widehat{f}_n$. • Estimate the conditional expected residual of the first-stage estimator, ${\mathbb{E}}\big[Y - f(X)\mid X = x_0\big]$, or equivalently, the estimation error $f_0(x_0) - \widehat{f}_n(x_0)$. Denote this estimator by $\widehat{b}_n(x_0)$. • Sum the first- and second-stage estimators to obtain \begin{align} &\widetilde{f}_n(x_0) \coloneqq \widehat{b}_n(x_0) + \widehat{f}_n(x_0), \end{align} which serves as a debiased estimator of $f_0$.

The term debiased estimator refers to our objective of correcting bias arising from the first-stage estimator. In many nonparametric and machine learning methods, bias can be problematic. By subtracting the estimated residual, ${\mathbb{E}}\big[Y - \widehat{f}_n(X) \mid X = x_0\big]$, we compensate for the difference between $f_0$ and $\widehat{f}_n$, resulting in an estimator with improved convergence and inferential properties.

Our estimator ((ref)) is closely related to the influence function in conditional mean estimation kennedy2023semiparametricdoublyrobusttargeted. It is also connected to the Neyman orthogonal score ChernozhukovVictor2018Dmlf, which is essentially equivalent to the canonical gradient in one-step bias correction Schuler2024. Notably, in conditional mean estimation, constructing the influence function is more intricate due to the conditioning. For an influence function that involves parameters represented via the conditional expected value, see Ichimura2022. Independently and concurrently, chernozhukov2024conditionalinfluencefunctions investigates a related topic with different motivations and methods.

The debiased nonparametric regression with local polynomial conditional expected residual estimation

In this section, we present an example of our debiased estimator. We employ local polynomial regression to estimate the conditional expected residual. For simplicity, we focus on the case where the covariate is one-dimensional ($d = 1$) and ${\mathcal{X}} = [0, 1]$.

First-stage nonparametric regression

For the first-stage estimator $\widehat{f}_n$, any regression method can be used, provided that smoothness conditions are met. As we show later, our key assumption is that the difference $f_0 - \widehat{f}_n$ is smooth. While various notions of smoothness exist, for simplicity, we focus on smoothness in the sense of the H\"older class (Definition (ref)). Notably, neither $f_0$ nor $\widehat{f}_n$ needs to be smooth individually, as long as their difference $f_0 - \widehat{f}_n$ is smooth.

We emphasize that our analysis does not impose any specific requirements on $\widehat{f}_n$ beyond smoothness. The first-stage estimator $\widehat{f}_n$ may converge to $f_0$ at an arbitrarily slow rate or even be inconsistent. Nevertheless, as long as the smoothness conditions hold, our asymptotic theoretical results remain valid. However, to achieve strong finite-sample performance, a high-quality first-stage estimator is typically preferred. In practice, methods such as random forests, neural networks, or other advanced machine learning techniques can be employed, provided they satisfy the smoothness conditions to a reasonable extent.

Second-stage conditional expected residual estimation

In the second-stage estimation of the conditional expected residual, we employ local polynomial regression. We define the $\ell$-th order polynomial basis as follows: \[ \rho(u) \coloneqq

pmatrix[pmatrix omitted — 59 chars of source]

^\top. \] In the theoretical analysis, the value of $\ell$ is chosen based on the smoothness of $f_0 - \widehat{f}_n$.

Using this polynomial basis, we define the local polynomial regression estimator as follows:

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

where $K\colon {\mathcal{X}} \to {\mathbb{R}}$ is a kernel function defined by \[ K(u) \coloneqq K_h(u) \coloneqq \mathbbm{1}\left[\big|u\big| \leq h\right]. \]

The conditional expected residual estimator is then given by the first component of $\widehat{\beta}_n(x)$: \[ \widehat{b}_n(x) = \widehat{\beta}_{0, n}(x). \]

Convergence analysis

This section presents a convergence analysis of $\widetilde{f}_n$ defined in Section (ref). We begin by defining the H\"older class.

definition[H\"older class] Given an interval ${\mathcal{T}} \subset {\mathbb{R}}$, and positive constants $\beta$ and $L$, the H\"older class $\Sigma(\beta, L)$ on ${\mathcal{T}}$ is the set of $\ell = \floor{\beta}$ times differentiable functions $f \colon {\mathcal{T}} \to {\mathbb{R}}$ whose $\ell$-times derivative $f^{(\ell)}$ satisfies \[ \Big|f^{(\ell)}(x) - f^{(\ell)}(z)\Big|\leq L \big|x - z\big|^{\beta - \ell}\quad \forall x, z \in {\mathcal{T}}. \]

Closed-form solution

The debiased estimator with local polynomial regression admits a closed-form solution. We define

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

Then, the second-stage estimator is given by \[ \widehat{b}_n(x_0) = \frac{1}{m}\sum_{i\in{\mathcal{I}}^{(2)}} \left(Y_i - \widehat{f}_n(X_i)\right) \widehat{w}_{h}(X_i, x_0). \] Thus, the debiased estimator is expressed as \[ \widetilde{f}_n(x_0) = \frac{1}{m}\sum_{i\in{\mathcal{I}}^{(2)}} \left(Y_i - \widehat{f}_n(X_i)\right) \widehat{w}_{h}(X_i, x_0) + \widehat{f}_n(x_0). \]

Bias and variance decomposition

We begin by analyzing the bias and variance of $\widetilde{f}_n$, defined as

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

The bias term quantifies how closely the estimator $\widetilde{f}_n$ centers around the true regression function $f_0$ at $x_0$, while the variance term measures the fluctuation of $\widetilde{f}_n$ around its own mean. Balancing bias and variance is a fundamental concern in statistical estimation. Our approach, which employs local polynomial smoothing in the second stage, aims to systematically reduce the bias introduced by the first-stage estimator while controlling variance through appropriate bandwidth selection.

To facilitate our analysis, we introduce the population version of $\widehat{{\mathcal{B}}}_{n}(x_0)$: \[ \widetilde{{\mathcal{B}}}_{n}(x_0) \coloneqq \int \rho\left(\frac{x - x_0}{h_n}\right) \rho\left(\frac{x - x_0}{h_n}\right)^\top K\left(\frac{x - x_0}{h_n}\right) {\mathrm{d}} F(x). \]

We then establish bounds for the bias and variance terms as follows:

theorem[Bias and variance decomposition] Let $s, L, C, C_1, C_2, C_3 > 0$ be constants independent of $f_0$ and $n$. Let $h_n$ be the bandwidth of the local polynomial estimator. For every $x \in {\mathcal{X}}$, the following hold: \begin{itemize} • The eigenvalues of $\widetilde{{\mathcal{B}}}_{n}(x)$ are uniformly bounded above and below away from zero. • There exists a constant $C > 0$ such that \[ \int \mathbbm{1}\big[| x - x_0 | \leq h_n\big]{\mathrm{d}} P(x) \leq C h_n, \] where $P(x)$ denotes the distribution of $x$. \end{itemize} If $f_0 - \widehat{f}_n$ belongs to the H\"older class $\Sigma(s, L)$ almost surely as $n\to \infty$, then for any $\varepsilon > 0$ and for all $x_0 \in {\mathcal{X}}$, there exists $n_0 > 0$ such that for all $n \geq n_0$, with probability at least $1 - \varepsilon$, it holds that \begin{align*} \bigl|\mathrm{Bias}(x_0)\bigr| &\leq C_1h^s_n,\\ \mathrm{Variance}(x_0) &\leq \frac{C_2}{nh_n} + C_3h^{2s}_n. \end{align*}

The smoothness condition on $f_0 - \widehat{f}_n$ can be satisfied in various ways; see Section (ref) for further discussion. The above result implies that as long as the smoothness of $f_0 - \widehat{f}_n$ is well-controlled and the bandwidth $h_n$ is chosen appropriately to shrink with $n$, both bias and variance remain bounded and ultimately vanish as $n$ grows.

Pointwise convergence of the MSE

From Theorem (ref), the bounds on bias and variance imply the convergence of the MSE, noting that \[ {\mathbb{E}}\left[\left(\widehat{f}_n(x_0) - f_0(x_0)\right)^2\right] = \mathrm{Bias}(x_0)^2 + \mathrm{Variance}(x_0). \] To minimize the MSE, we set \[ h_n = \alpha n^{-\frac{1}{2s + 1}}, \] where $\alpha$ is a constant independent of $f_0$ and $n$, $s$ is the smoothness parameter, and $d$ is the dimension of $X$. We then establish the pointwise MSE convergence result.

theorem[Pointwise MSE convergence] Let $s, L, C, \alpha > 0$ be constants independent of $f_0$ and $n$. Let \[ h_n = \alpha n^{-\frac{1}{2s + 1}} \] be the bandwidth of the local polynomial estimator. For every $x_0 \in {\mathcal{X}}$, the following hold: \begin{itemize} • The eigenvalues of $\widetilde{{\mathcal{B}}}_{n}(x)$ are uniformly bounded above and below away from zero. • There exists a constant $C > 0$ such that \[ \int \mathbbm{1}\big[| x - x_0 | \leq h_n\big]{\mathrm{d}} P(x) \leq C h_n, \] where $P(x)$ denotes the distribution of $x$. \end{itemize} If $f_0 - \widehat{f}_n$ belongs to the H\"older class $\Sigma(s, L)$ almost surely as $n\to \infty$, then for all $x_0 \in {\mathcal{X}}$, it holds that \[ {\mathbb{E}}\left[\left(\widetilde{f}_n(x_0) - f_0(x_0)\right)^2\right] = O\left(n^{-\frac{s}{2s + 1}}\right), \] as $n\to\infty$.

Theorem (ref) can be extended to establish the convergence of the MSE over the distribution of $X$.

corollary[MSE over the distribution of $X$] Suppose the conditions in Theorem (ref) hold. Then, we have \[ {\mathbb{E}}\left[\left(\widetilde{f}_n(X) - f_0(X)\right)^2\right] = O\left(n^{-\frac{2s}{2s + 1}}\right), \] where the expectation is taken over both $\widetilde{f}_n$ and $X$, as $n \to \infty$.

Pointwise asymptotic normality

We now establish the asymptotic normality of $\widetilde{f}(x_0)$. This result follows directly from Theorem (ref).

theorem[Asymptotic normality] Let $s, L, C > 0$ be constants independent of $f_0$ and $n$. Let $h_n$ be the bandwidth of the local polynomial estimator such that $h_n \to 0$, $nh_n \to \infty$, and $n h^{2s + 1}_n \to 0$ hold as $n\to \infty$. For every $x \in {\mathcal{X}}$, the following hold: \begin{itemize} • The eigenvalues of $\widetilde{{\mathcal{B}}}_{n}(x)$ are uniformly bounded above and below away from zero. • There exists a constant $C > 0$ such that \[ \int \mathbbm{1}\big[| x - x_0 | \leq h\big]{\mathrm{d}} P(x) \leq C h_n, \] where $P(x)$ is the distribution of $x$. \end{itemize} If $f_0 - \widehat{f}_n$ belongs to the H\"older class $\Sigma(s, L)$ almost surely as $n\to \infty$, then for every $x_0 \in {\mathcal{X}}$, it holds that \[ \sqrt{nh_n}\left(\widetilde{f}(x_0) - f_0(x_0)\right) \xrightarrow{{\mathrm{d}}} {\mathcal{N}}\big(0, V(x_0)\big), \] as $n \to \infty$, where $V(x_0) \geq 0$ is a constant depending on $x_0$ but independent of $n$.

Here, note that the following holds from the bias and variance decomposition: \[ V(x_0) \coloneqq nh_n\mathrm{Variance}(x_0) \leq C_2 + o(h_n). \]

The asymptotic normality result facilitates statistical inference, such as constructing confidence intervals or hypothesis tests for $f_0(x_0)$. In practice, knowing that $\sqrt{nh_n}(\widetilde{f}(x_0) - f_0(x_0))$ converges in distribution to a normal allows for classical inferential techniques, provided that the variance $V(x_0)$ can be consistently estimated. However, constructing a confidence interval can be challenging, because $s$ (the smoothness parameter) is often unknown and affects the choice of $h_n$. In real applications, data-driven methods, including cross-validation, can be employed to select $h_n$, and bootstrap or other resampling approaches can help in variance estimation.

Uniform convergence

Finally, we establish the uniform (sup-norm) convergence of $\widetilde{f}_n(x_0)$. Unlike mean squared convergence, which depends on the distribution of $X$, uniform convergence bounds the supremum of the estimation error over all $x \in {\mathcal{X}}$. This property is particularly valuable for addressing the covariate shift problem Shimodaira2000,SchmidtHieber2024, as uniform convergence remains unaffected by changes in the distribution of $X$ between training and test datasets.

theorem[Uniform convergence] Let $s, L_K, L, C, \alpha > 0$ be constants independent of $f_0$ and $n$. Let \[ h_n \coloneqq \alpha\left(\frac{\log(n)}{n}\right)^{\frac{1}{2s + 1}} \] be the bandwidth of the local polynomial estimator. For every $x_0 \in {\mathcal{X}}$, the following hold: \begin{itemize} • The eigenvalues of $\widetilde{{\mathcal{B}}}_{n}(x)$ are uniformly bounded above and below away from zero. • There exists a constant $C > 0$ such that \[ \int \mathbbm{1}\big[| x - x_0 | \leq h_n\big]{\mathrm{d}} P(x) \leq C h_n, \] where $P(x)$ denotes the distribution of $x$. • The noise $\varepsilon_i$ is sub-Gaussian, satisfying \[ {\mathbb{E}}\left[\exp(\lambda \varepsilon_i)\right] \leq \exp\left(K^2\lambda^2\right) \quad \forall \lambda \in {\mathbb{R}}, \] for some constant $K > 0$. • The kernel function $K$ is Lipschitz: $K \in \Sigma(1, L_K)$ on ${\mathcal{X}}$. \end{itemize} If $f_0 - \widehat{f}_n$ belongs to the H\"older class $\Sigma(s, L)$ almost surely as $n\to \infty$, then the following holds: \[ {\mathbb{E}}\left[\left\|\widetilde{f}_n - f_0\right\|^2_\infty\right] = O\left(\frac{\log(n)}{n}\right)^{\frac{2s}{2s + 1}}. \]

Although uniform convergence ensures robust performance under distribution shifts, it has not been established for many nonparametric regression estimators derived from modern machine learning methods. One reason for this gap is their data-adaptive nature and the reliance on empirical process techniques, which primarily focus on population risk measures whose expectation is taken over the covariate distribution rather than on sup-norm bounds. For example, SchmidtHieber2024 analyzes a neural network-based estimator in a one-dimensional covariate setting, demonstrating restricted uniform optimality. In contrast, our approach is model-free and accommodates multi-dimensional covariates, offering broader applicability and greater robustness.

Double robustness

Our proposed debiased estimator also possesses the property of double robustness: if either $\widehat{b}_n$ or $\widehat{f}_n$ is consistent, then $\widetilde{f}_n$ remains consistent.

Recall that \[ \widetilde{f}_n(x_0) = \frac{1}{m}\sum_{i \in {\mathcal{I}}^{(2)}} \left(Y_i - \widehat{f}_n(X_i)\right)\widehat{w}_h(X_i, x_0) + \widehat{f}_n(x_0). \] Suppose $\widehat{f}_n(x_0)$ converges in probability to $f^\dagger(x_0)$ and $\widehat{b}_n(x_0)$ converges in probability to $b^\dagger(x_0)$. We analyze two cases:

itemize• If $f^\dagger(x_0) = f_0(x_0)$ but $b^\dagger(x_0) \neq f_0(x_0) - f^\dagger(x_0)$, then \[ \widetilde{f}_n(x_0) = \frac{1}{m}\sum_{i \in {\mathcal{I}}^{(2)}}\Big(Y_i - f_0(X_i)\Big)w^\dagger(X_i, x_0) + f_0(x_0) + o_p(1), \] where $w^\dagger(X_i, x_0)$ depends only on $X_i$ and $x_0$. Since the first term converges in probability to zero, it follows that $\widetilde{f}_n(x_0) \xrightarrow{{\mathrm{p}}} f_0(x_0)$. • If $f^\dagger(x_0) \neq f_0(x_0)$ but $b^\dagger(x_0) = f_0(x_0) - f^\dagger(x_0)$, then \[ \widetilde{f}_n(x_0) = \left(f_0(x_0) - f^\dagger(x_0)\right) + f^\dagger(x_0) + o_p(1) = f_0(x_0) + o_p(1). \]

In both cases, $\widetilde{f}_n(x_0) \xrightarrow{{\mathrm{p}}} f_0(x_0)$. We summarize this result as follows.

theorem[Double robustness] For any $x_0 \in {\mathcal{X}}$, if either $\widehat{b}_n(x_0)$ or $\widehat{f}_n(x_0)$ is consistent, then \[ \widetilde{f}_n(x_0) \xrightarrow{{\mathrm{p}}} f_0(x_0)\qquad \text{as } n\to \infty. \]

Double robustness thus provides an additional safeguard against model misspecification or suboptimal performance of one of the estimators. In practice, this property is particularly valuable when the first-stage model $\widehat{f}_n$ fails to accurately capture the true function, but the second-stage residual model $\widehat{b}_n$ remains reliable—or vice versa.

Smoothness is all you need

As demonstrated above, our method relies solely on the smoothness of $\widehat{f} - f_0$. This assumption holds when both $f_0$ and $\widehat{f}$ are individually smooth, but it also remains valid in cases where neither $f_0$ nor $\widehat{f}$ is smooth, as long as their difference is smooth. Neural networks satisfy this condition under the analysis in SchmidtHieber2020.

Notably, our results require only the smoothness of the first-stage estimator (i.e., the difference $\widehat{f}_n - f_0$). This assumption significantly broadens the applicability of our method. Even when $f_0$ and $\widehat{f}_n$ are not individually smooth, our approach can still ensure MSE convergence, asymptotic normality, and uniform convergence, provided that the difference is smooth.

From a different perspective, this flexibility allows one to use state-of-the-art or specialized methods in the first stage without compromising the theoretical guarantees of the final estimator, as long as the first-stage error exhibits a level of smoothness akin to the H\"older condition. This feature is particularly useful when the true regression function $f_0$ has a complex or partially known structure, yet the difference $f_0 - \widehat{f}_n$ remains more regular due to design choices or inductive biases in $\widehat{f}_n$.

Conclusion

This study developed a debiased estimator for nonparametric regression. For any smooth regression estimator, our debiased approach ensures pointwise and uniform MSE convergence, asymptotic normality, and double robustness. The asymptotic normality facilitates statistical inference, while uniform convergence provides distributional robustness for nonparametric regression. Additionally, if either the first-stage or the second-stage nonparametric estimator is consistent, then the resulting regression estimator remains consistent, exhibiting the doubly robust property.

Notably, our estimator requires only the smoothness of the first-stage regression estimator. Specifically, if either (a) both the first-stage estimator $\widehat{f}$ and the true regression function $f_0$ are smooth, or (b) the difference $\widehat{f} - f_0$ is smooth, then the desirable properties of the debiased estimator hold. Thus, our method imposes minimal structural constraints while guaranteeing strong theoretical results, making it an attractive and flexible approach for a wide range of nonparametric regression problems.