# Debiased inference for a covariate-adjusted regression function

Kenta Takatsu  Department of Statistics and Data Science, Carnegie Mellon
University  Ted Westling  Department of Mathematics and Statistics, University
of Massachusetts Amherst

###### Abstract

In this article, we study nonparametric inference for a covariate-adjusted
regression function. This parameter captures the average association between a
continuous exposure and an outcome after adjusting for other covariates. Under
certain causal conditions, it also corresponds to the average outcome had all
units been assigned to a specific exposure level, known as the causal dose-
response curve. We propose a debiased local linear estimator of the covariate-
adjusted regression function and demonstrate that our estimator converges
pointwise to a mean-zero normal limit distribution. We use this result to
construct asymptotically valid confidence intervals for function values and
differences thereof. In addition, we use approximation results for the
distribution of the supremum of an empirical process to construct
asymptotically valid uniform confidence bands. Our methods do not require
undersmoothing, permit the use of data-adaptive estimators of nuisance
functions, and our estimator attains the optimal rate of convergence for a
twice differentiable regression function. We illustrate the practical
performance of our estimator using numerical studies and an analysis of the
effect of air pollution exposure on cardiovascular mortality.

###### Contents

  1. 1 Introduction
    1. 1.1 Motivation and literature review
    2. 1.2 Contribution and organization of the article
  2. 2 Proposed estimator
    1. 2.1 Notation and statistical setting
    2. 2.2 Parameter of interest and its interpretation
    3. 2.3 Debiased local linear estimator
    4. 2.4 Local parameter and its efficient influence function
    5. 2.5 Bandwidth selection
  3. 3 Asymptotic properties of the proposed methods
    1. 3.1 Pointwise convergence in distribution
    2. 3.2 Inference on causal effects
    3. 3.3 Uniform inference
  4. 4 Numerical studies
    1. 4.1 Study design
    2. 4.2 Results of numerical studies
  5. 5 Effect of fine particulate matter on cardiovascular mortality
  6. 6 Concluding remarks
  7. A Guide to notation
  8. B Conditions
  9. C Efficient influence function of the local parameter
  10. D Decomposition of the estimator
  11. E Proof of Theorems
  12. F Analysis of the leading term
  13. G Analysis of remainder term $R_{n,h,b,a_{0},1}$
  14. H Analysis of remainder terms $R_{n,h,b,a_{0},2}$ and $R_{n,h,b,a_{0},3}$
  15. I Analysis of remainder term $R_{n,h,b,a_{0},4}$
  16. J Analysis of remainder term $R_{n,h,b,a_{0},5}$
  17. K Analysis of remainder term $R_{n,h,b,a_{0},6}$
  18. L Analysis of the covariance estimator
  19. M Lemmas supporting uniform result
  20. N Additional results from numerical studies

##  1 Introduction

###  1.1 Motivation and literature review

In this article, we study nonparametric inference for a covariate-adjusted
regression function, which is also known as a G-computed regression function.
This statistical problem arises in the context of observational studies where
interest focuses on the causal effect of a continuous exposure, such as the
dose of a drug, the amount of air pollution exposure, or the amount of a
biochemical in the bloodstream. However, a covariate-adjusted regression
function may also be of interest outside of causal contexts as a one-
dimensional marginal summary of a multivariate regression function. Despite
the simple formulation and many applications of this parameter, a method
achieving valid nonparametric inference without undersmoothing is not yet
available. In this paper, we close this gap by introducing a novel
nonparametric, doubly-robust estimator, pointwise confidence intervals, and
uniform confidence bands for the covariate-adjusted regression function that
do not require undersmoothing.

One setting where a covariate-adjusted regression function arises is causal
inference with a continuous exposure. The gold standard for assessing the
causal effect of a treatment or exposure on an outcome is a randomized
experiment, where units in the population are assigned values of the exposure
by a random process known to the researchers. Frequently, however, such an
experiment is infeasible, unethical, or cost-prohibitive. For example, it is
unethical to purposefully expose people to a chemical or pollutant known to
have negative health effects. In such cases, researchers may instead wish to
assess the causal effect of the exposure using observational data in which the
exposure varies according to an unknown mechanism. Recovering a causal effect
with observational data is more challenging due to potential common causes of
the exposure and outcome. However, if all common causes are recorded in the
data, then causal effects can be recovered by appropriately adjusting for
them. Specifically, the average outcome had all units been assigned a specific
exposure value coincides with the covariate-adjusted regression function,
which is known as the G-formula or G-computation in causal inference (Robins,,
1986; Gill and Robins,, 2001).

The covariate-adjusted regression function is also of interest outside of
causal contexts. The regression (i.e., the conditional expectation) of an
outcome on a vector of covariates can be difficult to visualize and summarize
when there are more than two regressors and when nonparametric methods are
used. The covariate-adjusted regression function summarizes the adjusted
association between a single continuous covariate and the outcome by averaging
the regression function over all other covariates for each value of the
covariate of interest. This is related to the use of _marginal effects_ to
summarize the results of nonlinear regression models (Mize et al.,, 2019).
Furthermore, since the covariate-adjusted regression is univariate, it is easy
to visualize.

Several nonparametric methods exist for this problem. Kennedy et al., (2016)
proposed an estimator based on local linear smoothing. Theirs was the first
_doubly robust_ estimator for this parameter, meaning that their estimator is
consistent if either of two nuisance estimators is consistent. van der Laan et
al., (2018) and Colangelo and Lee, (2020) also proposed doubly-robust
estimators based on smoothing, and used cross-fitting to remove empirical
process conditions for nuisance estimators. Semenova and Chernozhukov, (2021)
proposed an estimator based on a series expansion. Finally, Westling and
Carone, (2020) and Westling et al., (2020) explored inference under a
monotonicity assumption. Additional literature related to the covariated-
adjusted regression function includes Robins, (1986, 2000) and Zhang et al.,
(2016), who proposed plug-in estimators based on a parametric outcome
regression model; Hirano and Imbens, (2004) and Imai and Van Dyk, (2004), who
proposed estimators based on a parametric propensity model; Neugebauer and van
der Laan, (2007), who studied inference for the best projection onto a working
parametric model; Rubin and van der Laan, (2006) and Díaz and van der Laan,
(2013), who proposed data-adaptive methods; Bonvini and Kennedy, (2022), who
proposed a higher-order estimator; and Westling, (2022) and Weng et al.,
(2022), who proposed tests of the null hypothesis that the function is flat.

As is the case for an ordinary regression function, nonparametric inference
for the covariate-adjusted regression function is a challenging task. The bias
of smoothing-based methods is not asymptotically negligible when the bandwidth
is chosen to optimize the rate of convergence of the estimator. This bias
complicates the task of obtaining valid inference. Some authors have argued
for interpreting confidence sets constructed based on the resulting limit
theorems as valid for a smoothed parameter (Wasserman,, 2006; Kennedy et al.,,
2016). An alternative approach is to choose the bandwidth to go to zero faster
than the optimal rate to guarantee that the bias goes to zero faster than the
variance, which is called undersmoothing (van der Laan et al.,, 2018;
Colangelo and Lee,, 2020; Semenova and Chernozhukov,, 2021). While this
approach theoretically allows valid inference, it yields a suboptimal rate of
convergence for the estimator. Furthermore, there is little guidance about how
to select a bandwidth for undersmoothing in practice beyond ad-hoc methods.

Recently, Calonico et al., (2018) proposed a method of bias correction for
density and regression functions estimated using kernel smoothing. Bias
correction in the context of kernel smoothing is challenging because the bias
depends on the second or higher derivative of the true function, which is more
difficult to estimate than the function itself (Wasserman,, 2006; Hall,,
1992). However, Calonico et al., (2018) demonstrated that, via a careful
choice of the bandwidth parameter of the bias estimator, it is possible to
effectively debias kernel smoothing-based estimators. These estimators permit
asymptotically valid inference without undersmoothing, and unlike
undersmoothing, retain the optimal rate of convergence relative to the assumed
smoothness of the true function.

###  1.2 Contribution and organization of the article

In this article, we contribute to the existing literature in the following
ways: (1) we propose a novel debiased estimator of the covariate-adjusted
regression function motivated by the approach proposed by Calonico et al.,
(2018); (2) we propose methods of pointwise and uniform inference and provide
conditions under which our methods are asymptotically valid; and (3) we
illustrate the practical performance of our proposed methods using numerical
studies and an analysis of the effect of air quality on health. To the best of
our knowledge, ours are the first asymptotically valid methods of pointwise
and uniform inference for the covariate-adjusted regression function without
undersmoothing. We note that our problem is substantively different from that
of Calonico et al., (2018) due to the presence of nuisance parameters, which
introduces remainder terms and technical considerations not present when
estimating density and regression functions. We elucidate these differences
more below.

The remainder of this article is organized as follows. In Section 2, we define
our statistical setting and proposed estimator. In Section 3, we present our
approach to inference and our theoretical results. In Section 4, we
demonstrate the empirical performance of our proposed methods using numerical
studies, and in Section 5, we use our methods to assess the effect of air
pollution on cardiovascular mortality. Section 6 presents brief concluding
remarks. Proofs of all theorems are provided in Supplementary Material.

##  2 Proposed estimator

###  2.1 Notation and statistical setting

We now define the statistical setting we will work in and notation we will
use. We consider a univariate outcome $Y\in\mathcal{Y}\subset\mathbb{R}$, a
univariate exposure $A\in\mathcal{A}\subset\mathbb{R}$, and a vector of
covariates $W\in\mathcal{W}\subseteq\mathbb{R}^{d}$. We assume that the
distribution of $A$ possesses a Lebesgue density. We define the observed data
unit $O\coloneqq(Y,A,W)$, which takes values in
$\mathcal{O}\coloneqq\mathcal{Y}\times\mathcal{A}\times\mathcal{W}$. We then
observe $n$ independent and identically distributed observations
$O_{1},\dots,O_{n}$ from an unknown probability distribution $P_{0}$. We
denote by $\mathbb{P}_{n}$ the empirical distribution function of the observed
data. We index objects by $P$ when they depend on a generic distribution $P$
over the observed data unit $O$, and we use subscript $0$ as short-hand for
the true distribution $P_{0}$. For instance, we denote the expectation under
$P_{0}$ by $E_{0}$. For a distribution $P$ on $\mathcal{O}$, we denote by
$F_{P}$ the marginal distribution and $f_{P}$ the Lebesgue density of $A$
under $P$. We denote by $Q_{P}$ the marginal distribution of $W$ under $P$. We
denote by $F_{n}$ and $Q_{n}$ the empirical distributions of
$A_{1},\dotsc,A_{n}$ and $W_{1},\dotsc,W_{n}$, respectively. We let
$\mu_{P}(a,w)\coloneqq E_{P}(Y\mid A=a,W=w)$ denote the _outcome regression_
function and $g_{P}(a,w)\coloneqq\left[\frac{\partial}{\partial a}P(A\leq
a\mid W=w)\right]/f_{P}(a)$ denote the _standardized propensity_ function. For
a probability measure $P$ and $P$-integrable function $h$, we define
$Ph\coloneqq\int h\,dP$. For $q\geq 1$, we denote by
$\|h\|_{P,q}\coloneqq(P|h|^{q})^{1/q}$ the $L_{q}(P)$ norm of $h$. We also
define $e_{j}$ as the ($j+1$)-vector of zero’s with one at the $j$th position
for any integer $j\geq 1$.

###  2.2 Parameter of interest and its interpretation

Our parameter of interest is the covariate-adjusted regression function
$a\mapsto\theta_{0}(a)$ defined as

$\theta_{0}(a)\coloneqq E_{0}\left\\{E_{0}(Y\mid
A=a,W)\right\\}=E_{0}\left\\{\mu_{0}(a,W)\right\\}=\int\mu_{0}(a,w)\,dQ_{0}(w).$

Under certain conditions, $\theta_{0}$ has a causal interpretation.
Informally, if the covariates $W$ contain all common causes of $A$ and $Y$,
exposure value $A=a$ is in the support of the conditional distribution of $A$
given $W=w$ for almost all $w$, and one unit’s exposure doesn’t impact other
units’ outcomes, then $\theta_{0}(a)$ can be interpreted as the average
outcome value had all units in the population been assigned to exposure value
$A=a$. This latter parameter is called the _causal dose-response curve_ or the
_average dose-response function_. We refer readers interested in a formal
statement of these conditions and result and additional discussion of this
causal parameter to Robins, (1986), Gill and Robins, (2001), Kennedy et al.,
(2016), Westling et al., (2020), and Westling, (2022).

As mentioned in the introduction, $\theta_{0}$ is also of interest outside of
causal settings. The outcome regression function $\mu_{0}$ is the expected
outcome given exposure and covariates. Hence, for fixed $a$, $\mu_{0}(a,W)$ is
a random variable representing the expected outcome value given $A=a$ across
the distribution of the covariates $W$ in the population, and $\theta_{0}(a)$
is the mean of this variable. The curve $a\mapsto\theta_{0}(a)$ depicts how
this average conditional mean changes with $a$. Hence, $\theta_{0}$ is a
marginal summary of the multivariate regression function $\mu_{0}$.

###  2.3 Debiased local linear estimator

We now define our estimator of the covariate-adjusted regression function. We
begin with a review of the local linear method proposed by Kennedy et al.,
(2016) and its key properties. For an outcome regression function $\mu$, a
standardized propensity $g$, and a covariate distribution $Q$, we define the
_pseudo-outcome_ mapping

$\xi_{\mu,g,Q}:(y,a,w)\mapsto\frac{y-\mu(a,w)}{g(a,w)}+\int\mu(a,\bar{w})\,dQ(\bar{w}).$
(1)

Theorem 1 of Kennedy et al., (2016) showed that this mapping possesses a
double-robust property: it holds that $E_{0}[\xi_{\mu,g,Q_{0}}(Y,A,W)\mid
A=a_{0}]=\theta_{0}(a_{0})$ if either $\mu=\mu_{0}$ or $g=g_{0}$. Kennedy et
al., (2016) thus proposed first constructing estimators $\mu_{n}$ and $g_{n}$
of $\mu_{0}$ and $g_{0}$, respectively, using the empirical distribution
$Q_{n}$ of the observed covariates to estimate $Q_{0}$, then regressing the
estimated pseudo-outcomes
$\xi_{n}(Y_{1},A_{1},W_{1}),\dotsc,\xi_{n}(Y_{n},A_{n},W_{n})$ on the observed
exposures $A_{1},\dotsc,A_{n}$ using local linear regression, where
$\xi_{n}\coloneqq\xi_{\mu_{n},g_{n},Q_{n}}$. Specifically, let $K$ be a kernel
function (i.e., a symmetric density on $\mathbb{R}$), $h>0$ be a bandwidth,
and $K_{h,a_{0}}(a)\coloneqq K([a-a_{0}]/h)/h$. The local linear estimator at
a point $a_{0}$ is then defined as the intercept of the weighted ordinary
least squares regression of the estimated pseudo-outcomes on the observed
exposures with weights $K_{h,a_{0}}(A_{1}),\dotsc,K_{h,a_{0}}(A_{n})$.
Mathematically, the local linear estimator can be written as

$\displaystyle\theta_{n,h}^{LL}(a_{0})\coloneqq
e_{1}^{T}\mathbf{D}_{n,h,a_{0},1}^{-1}\mathbb{P}_{n}\left(w_{h,a_{0},1}K_{h,a_{0}}\xi_{n}\right),$
(2)

where
$w_{h,a_{0},j}(a)\coloneqq\left(1,[a-a_{0}]/h,\dotsc,[a-a_{0}]^{j}/h^{j}\right)^{T}$
and
$\mathbf{D}_{n,h,a_{0},j}\coloneqq\mathbb{P}_{n}(w_{h,a_{0},j}w^{T}_{h,a_{0},j}K_{h,a_{0}})$
for any integer $j\geq 1$. Other authors have used a similar approach as
Kennedy et al., (2016), but replaced the local linear regression step with an
alternative nonparametric regression (Westling et al.,, 2020; Semenova and
Chernozhukov,, 2021; Bonvini and Kennedy,, 2022).

As discussed in the introduction, standard approaches to nonparametric
regression, including local linear regression, do not yield valid inference
when the bandwidth is chosen to minimize mean squared error because the bias
of the resulting estimator is of the same order as its standard deviation.
Specifically, in their Theorem 3, Kennedy et al., (2016) showed that under
suitable conditions,
$(nh)^{1/2}\big{[}\theta_{n,h}^{LL}(a_{0})-\theta_{0}(a_{0})-\tfrac{1}{2}h^{2}c_{2}\theta_{0}^{(2)}(a_{0})\big{]}$
converges to a mean-zero normal distribution, where $c_{2}\coloneqq\int
u^{2}K(u)\,du$ and $\theta_{0}^{(2)}$ is the second derivative of
$\theta_{0}$. To minimize mean squared error, the bandwidth $h$ should be
chosen to balance bias squared and variance, which means choosing $h$ such
that $(nh)^{1/2}h^{2}=(nh^{5})^{1/2}$ converges to a positive, finite
constant, or equivalently $h$ proportional to $n^{-1/5}$. Hence, if
$\theta_{0}^{(2)}(a_{0})\neq 0$, then
$(nh)^{1/2}\big{[}\theta_{n,h}^{LL}(a_{0})-\theta_{0}(a_{0})\big{]}$ converges
to a normal distribution with non-zero mean, implying that confidence
intervals centered around $\theta_{n,h}^{LL}(a_{0})$ will be asymptotically
anti-conservative.

We propose debiasing the local linear estimator by subtracting an estimator of
the bias in the spirit of Calonico et al., (2018). Our debiased estimator will
be of the form
$\theta_{n,h,b}(a_{0}):=\theta_{n,h}^{LL}(a_{0})-\tfrac{1}{2}h^{2}c_{2}\theta_{n,b}^{(2)}(a_{0})$,
where $\theta_{n,b}^{(2)}$ is a third-order local polynomial estimator of
$\theta_{0}^{(2)}$ with bandwidth $b>0$ (Fan and Gijbels,, 1996). Using the
notation defined above, we have
$\theta_{n,b}^{(2)}(a_{0})=2b^{-2}e_{3}^{T}\mathbf{D}_{n,b,a_{0},3}^{-1}\mathbb{P}_{n}(w_{b,a_{0},3}K_{b,a_{0}}\xi_{n})$
so that we can write
$\theta_{n,h,b}(a_{0})=\mathbb{P}_{n}(\Gamma_{n,h,b,a_{0}}\xi_{n})$ for

$\displaystyle\Gamma_{n,h,b,a_{0}}(a)\coloneqq
e_{1}^{T}\mathbf{D}_{n,h,a_{0},1}^{-1}w_{h,a_{0},1}(a)K_{h,a_{0}}(a)-e_{3}^{T}c_{2}(h/b)^{2}\mathbf{D}_{n,b,a_{0},3}^{-1}w_{b,a_{0},3}(a)K_{b,a_{0}}(a).$
(3)

To summarize, for given bandwidths $h$ and $b$ and kernel $K$, our proposed
estimator is constructed in two steps: (1) construct estimators $\mu_{n}$ and
$g_{n}$ of $\mu_{0}$ and $g_{0}$ respectively, and (2) compute the plug-in
estimates of pseudo-outcomes
$\xi_{n}(Y_{i},A_{i},W_{i})=\xi_{\mu_{n},g_{n},Q_{n}}(Y_{i},A_{i},W_{i})$ for
$i=1,\dotsc,n$, and regress them on the observed exposures using the bias-
corrected local linear estimator
$\theta_{n,h,b}(a_{0})=\mathbb{P}_{n}(\Gamma_{n,h,b,a_{0}}\xi_{n})$ for each
$a_{0}$. In Section 3, we provide conditions on the true data-generating
mechanism and on $h$, $b$, $K$, $\mu_{n}$, and $g_{n}$, as well as practical
guidance for selecting or estimating these quantities.

It may seem that basing the second derivative estimator $\theta_{n,b}^{(2)}$
on a third-order local polynomial estimator would require additional
smoothness of $\theta_{0}$, hence violating our stated goal of obtaining the
optimal rate of convergence relative to the assumed smoothness of
$\theta_{0}$. This is not the case. As an intuitive explanation, we decompose
$(nh)^{1/2}\left[\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})\right]$ as

$(nh)^{1/2}\left[\theta_{n,h}^{LL}(a_{0})-\theta_{0}(a_{0})-\tfrac{1}{2}h^{2}c_{2}\theta_{0}^{(2)}(a_{0})\right]+\tfrac{1}{2}c_{2}(nh^{5})^{1/2}\left[\theta_{n,b}^{(2)}(a_{0})-\theta_{0}^{(2)}(a_{0})\right].$

As shown in Kennedy et al., (2016), the first term converges in distribution
to a mean-zero limit. Thus, if $nh^{5}=O(1)$ and
$\theta_{n,b}^{(2)}(a_{0})\operatorname*{\stackrel{{\scriptstyle
P}}{{\longrightarrow}}}\theta_{0}^{(2)}(a_{0})$, then
$(nh)^{1/2}\left[\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})\right]$ converges to
this same limit. However, perhaps surprisingly, $\theta_{n,b}^{(2)}(a_{0})$
need not be consistent for $\theta_{0}^{(2)}(a_{0})$ to achieve good inference
using $\theta_{n,h,b}(a_{0})$. Roughly speaking, if
$\theta_{n,b}^{(2)}(a_{0})-\theta_{0}^{(2)}(a_{0})$ merely converges in
distribution to a mean-zero limit, then so does
$(nh)^{1/2}\left[\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})\right]$. The variance
of $\theta_{n,b}^{(2)}(a_{0})$ is proportional to $(nb^{5})^{-1}$, so if
$h/b\longrightarrow\tau>0$, then the bias correction _will not_ be
asymptotically negligible; it will contribute to the asymptotic variance of
the estimator as in Calonico et al., (2018). If in addition $h\propto
n^{-1/5}$, then the variance of $\theta_{n,b}^{(2)}(a_{0})$ does not go to
zero so it is not consistent.

Even if the bias correction is asymptotically negligible, and especially if it
is not, accounting for its finite-sample variability is important for
achieving good finite-sample inference. As discussed more below, our variance
estimator will account for the variability of $\theta_{n,b}^{(2)}(a_{0})$, and
in particular, we can still achieve valid inference when the variance of
$\theta_{n,b}^{(2)}(a_{0})$ is not going to zero. We will show that
$\theta_{n,b}^{(2)}(a_{0})-\theta_{0}^{(2)}(a_{0})$ converges to a mean-zero
limit if $b\longrightarrow 0$ and $\theta_{0}^{(2)}$ is continuous at $a_{0}$
(and additional technical conditions unrelated to smoothness of $\theta_{0}$
hold).

###  2.4 Local parameter and its efficient influence function

We now provide an alternative motivation for our proposed estimator, which
also motivates our approach to inference. We show that $\theta_{n,h,b}(a_{0})$
can be considered as a one-step estimator of a debiased smoothed parameter.

The parameter mapping $P\mapsto\theta_{P}(a_{0})$ is not _pathwise
differentiable_ relative to a nonparametric model, meaning that it is not
smooth enough as a function of the distribution $P$ to permit $n^{-1/2}$-rate
estimation (Bickel,, 1982; Pfanzagl and Wefelmeyer,, 1985; Bickel and
Klaassen,, 1993). One way to develop inference methods for such a parameter is
to consider a sequence of smoothed parameters approaching the parameter of
interest, each of which is pathwise differentiable (van der Laan et al.,,
2018). Our debiased local linear estimator can be viewed through this lens.
For any integer $j\geq 1$, we define $\mathbf{D}_{P,h,a_{0},j}\coloneqq
P(w_{h,a_{0},j}w_{h,a_{0},j}^{T}K_{h,a_{0}})$. We then define the _debiased
local parameter_ $\theta_{P,h,b}(a_{0})\coloneqq
P(\Gamma_{P,h,b,a_{0}}\theta_{P})$, where

$\displaystyle\Gamma_{P,h,b,a_{0}}(a)\coloneqq
e_{1}^{T}\mathbf{D}_{P,h,a_{0},1}^{-1}w_{h,a_{0},1}(a)K_{h,a_{0}}(a)-e_{3}^{T}c_{2}(h/b)^{2}\mathbf{D}_{P,b,a_{0},3}^{-1}w_{b,a_{0},3}(a)K_{b,a_{0}}(a).$
(4)

The asymptotic properties of $\theta_{n,h,b}(a_{0})$ can be understood through
the decomposition

$(nh)^{1/2}\left[\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})\right]=(nh)^{1/2}\left[\theta_{n,h,b}(a_{0})-\theta_{0,h,b}(a_{0})\right]+(nh)^{1/2}\left[\theta_{0,h,b}(a_{0})-\theta_{0}(a_{0})\right].$

In Supplementary Material, we show that, if $\theta_{0}$ is twice continuously
differentiable in a neighborhood of $a_{0}$,
$h/b\longrightarrow\tau\in[0,\infty)$, and additional mild conditions hold,
then $\theta_{0,h,b}(a_{0})-\theta_{0}(a_{0})=o(h^{2})$ as $h\longrightarrow
0$. Hence, if $nh^{5}=O(1)$, then
$\theta_{0,h,b}(a_{0})-\theta_{0}(a_{0})=o(\\{nh\\}^{-1/2})$, and so the
first-order asymptotic properties of
$(nh)^{1/2}\left[\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})\right]$ are
determined by
$(nh)^{1/2}\left[\theta_{n,h,b}(a_{0})-\theta_{0,h,b}(a_{0})\right]$. This
expression can be studied using semiparametric efficiency theory. The first
step in doing so is to derive the efficient influence function of
$\theta_{0,h,b}(a_{0})$. This is the subject of the following lemma.

######  Lemma 1 (Efficient influence function).

For each $h$ and $b>0$ and $a_{0}\in\mathcal{A}$,
$P\mapsto\theta_{P,h,b}(a_{0})$ is a pathwise differentiable parameter with
respect to the model $\mathcal{M}$ consisting of $P$ such that (1)
$E_{P}[Y^{2}]<\infty$ and (2) there exists $\kappa>0$ and
$\mathcal{W}_{0}\subseteq\mathcal{W}$ with $Q_{P}(\mathcal{W}_{0})=1$ such
that $g_{P}(a,w)\geq\kappa$ for all $(a,w)$ such that
$|a-a_{0}|\leq\max\\{h,b\\}$ and $w\in\mathcal{W}_{0}$, and the efficient
influence function of $\theta_{P,h,b}(a_{0})$ relative to this model is

$\displaystyle\phi^{*}_{P,h,b,a_{0}}:(y,a,w)$
$\displaystyle\mapsto\Gamma_{P,h,b,a_{0}}(a)\xi_{\mu_{P},g_{P},Q_{P}}(y,a,w)-\gamma_{P,h,b,a_{0}}(a)$
$\displaystyle\qquad+\int\Gamma_{P,h,b,a_{0}}(\bar{a})\left\\{\mu_{P}(\bar{a},w)-\theta_{P}(\bar{a})\right\\}\,dF_{P}(\bar{a}),\text{
where}$ (5) $\displaystyle\gamma_{P,h,b,a_{0}}(a)$
$\displaystyle:=e_{1}^{T}\mathbf{D}^{-1}_{P,h,a_{0},1}w_{h,a_{0},1}(a)K_{h,a_{0}}(a)w^{T}_{h,a_{0},1}(a)\mathbf{D}^{-1}_{P,h,a_{0},1}P\left(w_{h,a_{0},1}K_{h,a_{0}}\theta_{P}\right)$
$\displaystyle\qquad-
e_{3}^{T}c_{2}(h/b)^{2}\mathbf{D}_{P,b,a_{0},3}^{-1}w_{b,a_{0},3}(a)K_{b,a_{0}}(a)w^{T}_{b,a_{0},3}(a)\mathbf{D}^{-1}_{P,b,a_{0},3}P\left(w_{b,a_{0},3}K_{b,a_{0}}\theta_{P}\right).$

The proof of Lemma 1 and all other results are provided in Supplementary
Material.

We now let $\phi_{n,h,b,a_{0}}^{*}$ be the estimated efficient influence
function obtained by replacing $\mu_{P}$ and $g_{P}$ with estimators $\mu_{n}$
and $g_{n}$, $Q_{P}$ and $F_{P}$ with the empirical distributions $Q_{n}$ and
$F_{n}$, and $\mathbf{D}_{P,h,a_{0},j}$ with $\mathbf{D}_{n,h,a_{0},j}$. Due
to the cancellation of terms in the influence function, it then holds that
$\theta_{n,h,b}(a_{0})=\int\mu_{n}(a_{0},w)\,dQ_{n}(w)+\mathbb{P}_{n}\phi_{n,h,b,a_{0}}^{*}$;
i.e., the debiased estimator $\theta_{n,h,b}(a_{0})$ is a one-step estimator
of the debiased smoothed parameter $\theta_{0,h,b}(a_{0})$.

The representation of our debiased estimator as a one-step estimator of a
smoothed parameter plays an important role in motivating our approach to
inference and our asymptotic results. We would typically expect that one-step
estimators of pathwise differentiable parameters are asymptotically linear
under appropriate conditions. Similarly, in Theorem 1, we will see that a
finite-sample version of asymptotic linearity holds for our estimator:
$\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})=\mathbb{P}_{n}\phi_{\infty,h,b,a_{0}}^{*}+o_{p}(\\{nh\\}^{-1/2})$.
As with asymptotic linearity, this representation is useful because it reduces
the derivation of further asymptotic properties to the study of the linear
term $\mathbb{P}_{n}\phi_{\infty,h,b,a_{0}}^{*}$. Furthermore, it suggests the
variance estimator
$\sigma_{n,h,b}^{2}(a_{0}):=h\mathbb{P}_{n}(\phi_{\infty,h,b,a_{0}}^{*})^{2}$.
Importantly, this variance estimator accounts for the contribution of the bias
estimator, and in particular it is a consistent estimator of the asymptotic
variance even when the bias estimator contributes to the asymptotic variance.

###  2.5 Bandwidth selection

As with the debiased density and regression estimators proposed by Calonico et
al., (2018), our estimator requires the choice of two bandwidths. Many
bandwidth selection methods for local polynomial regression can be adapted to
our setting. Here, we briefly discuss several strategies for data-driven
bandwidth selection. In Section 4, we compare the empirical behavior of the
three methods outlined here.

First, we can choose the bandwidths to minimize a cross-validated estimator of
the integrated mean squared error (IMSE) of $\theta_{n,h,b}$, as Kennedy et
al., (2016) did for the local linear estimator. Specifically, since
$\theta_{n,h,b}(a_{0})$ can be written as a linear smoother of $\xi_{n}$, a
computationally efficient leave-one-out cross-validated estimator of the IMSE
of $\theta_{n,h,b}$ is given by

$\displaystyle\mathrm{IMSE}_{\mathrm{cv}}(h,b):=\frac{1}{n}\sum_{i=1}^{n}\left\\{\frac{\xi_{n}(Y_{i},A_{i},W_{i})-\theta_{n,h,b}(A_{i})}{1-\Gamma_{n,h,b,A_{i}}(A_{i})/n}\right\\}^{2}.$
(6)

We refer the reader to Chapter 5.3, and specifically Theorem 5.34, of
Wasserman, (2006) for additional discussion of this formula. We then define
$(h_{\text{cv}},b_{\text{cv}}):=\operatorname*{arg\,min}_{h,b}\mathrm{IMSE}_{\mathrm{CV}}(h,b)$.
Alternatively, we can fix $b$ as a function of $h$ via
$h/b=\tau\in(0,\infty)$. For example, we can fix $b=h$ so that $\tau=1$. We
can then optimize the estimated IMSE over $h$ alone, i.e.
$h_{\text{cv},\tau}:=\operatorname*{arg\,min}_{h}\mathrm{IMSE}_{\mathrm{CV}}(h,h/\tau)$.
Fixing $\tau$ removes the need to select $b$, which reduces the search to a
one-dimensional space, and guarantees that $h/b=O(1)$, which is required by
our conditions below. Furthermore, in the setting of regression estimation,
Calonico et al., (2018) found that fixing $\tau>0$ to a positive constant
yields improved coverage accuracy. However, the choice of $\tau$ is somewhat
arbitrary, and as we will see in Theorem 1, $\tau>0$ also yields an estimator
with larger asymptotic variance than that of Kennedy et al., (2016).

As a second bandwidth selection procedure, we will consider an adaptation of
the plug-in method proposed by Calonico et al., (2018) for their debiased
local linear estimator of a regression function. The method works by
minimizing an estimator of the IMSE of the local linear estimator
$\theta_{n,h}^{LL}$ over $h$. We approximate the large-sample IMSE of
$\theta_{n,h}^{LL}$ with respect to a probability measure $\omega$ as
$\mathrm{IMSE}_{\mathrm{plug-
in}}(h):=h^{4}\int\hat{B}^{LL}(a)^{2}d\omega(a)+(nh)^{-1}\int\hat{V}^{LL}(a)\,d\omega(a)$,
where

$\displaystyle\hat{B}^{LL}(a_{0})$
$\displaystyle:=e_{1}^{T}\mathbf{D}_{n,h_{1},a_{0},1}^{-1}\frac{1}{n}\sum_{i=1}^{n}w_{h_{1},a_{0},1}(A_{i})K_{h_{1},a_{0}}(A_{i})\tfrac{1}{2}\hat{\theta}_{h_{1}}^{(2)}(a_{0})\left(\frac{A_{i}-a_{0}}{h_{1}}\right)^{2},\text{
and}$ $\displaystyle\hat{V}^{LL}(a_{0})$
$\displaystyle:=e_{1}^{T}\mathbf{D}_{n,h_{1},a_{0},1}^{-1}\left\\{\frac{1}{n}\sum_{i=1}^{n}w_{h_{1},a_{0},1}(A_{i})K_{h_{1},a_{0}}(A_{i})^{2}\hat{\sigma}^{2}(A_{i})w^{T}_{h_{1},a_{0},1}(A_{i})\right\\}\mathbf{D}_{n,h_{1},a_{0},1}^{-1}e_{1}.$

Here, $h_{1}$ is a _pilot bandwidth_ , $\hat{\theta}_{h_{1}}^{(2)}$ is an
estimator of $\theta_{0}^{(2)}$, and $\hat{\sigma}^{2}(a)$ is a nearest-
neighbors estimator of the conditional variance of $\xi_{n}$ given $A=a$.
Since by design neither $\hat{B}^{LL}$ nor $\hat{V}^{LL}$ depends on $h$, the
bandwidth $h_{\mathrm{plug-in}}$ minimizing $\mathrm{IMSE}_{\mathrm{plug-in}}$
is given explicitly by

$h_{\mathrm{plug-
in}}:=n^{-1/5}\left(\frac{\int\hat{V}^{LL}\,d\omega}{4\int[\hat{B}^{LL}]^{2}d\omega}\right)^{1/5}.$

Finally, the bandwidth $b_{\text{plug-in}}$ of the bias correction is defined
as $b_{\text{plug-in}}=h_{\text{plug-in}}/\tau$, where $\tau$ is user-
specified. A benefit of this method is that it does not require numerical
optimization. A second benefit is that the measure $\omega$ can be chosen
based on the range over which it is of interest to estimate $\theta_{0}$.

##  3 Asymptotic properties of the proposed methods

###  3.1 Pointwise convergence in distribution

In this section, we study the asymptotic properties of our proposed estimator,
and use these properties to derive approaches to pointwise and uniform
inference. We begin by showing that
$(nh)^{1/2}[\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})]$ converges in
distribution to a normal limit for fixed $a_{0}$. Using this result, we will
then show that the $(1-\alpha)$-level Wald-style confidence interval

$\theta_{n,h,b}(a_{0})\pm(nh)^{-1/2}q_{1-\alpha/2}\sigma_{n,h,b}(a_{0}),$ (7)

for $q_{p}$ the $p$th quantile of a standard normal distribution has
asymptotic coverage level $1-\alpha$ of $\theta_{0}(a_{0})$ without
undersmoothing.

We begin by introducing technical conditions we will rely upon. We discuss
these conditions following the statement of Theorem 1. Our first two
conditions concern the kernel function and bandwidths, which will be required
in all of our results.

  1. (A1)

The kernel $K$ is a mean-zero, symmetric, nonnegative, and Lipschitz
continuous density function with support contained in $[-1,1]$. Additionally,
$K$ belongs to the linear span of the functions whose subgraph can be
represented as a finite number of Boolean operations among sets of the form
$\\{(s,u)\in\mathbb{R}\times\mathbb{R}:p(s,u)\leq\varphi(u)\\}$ where $p$ is a
polynomial and $\varphi$ is an arbitrary real function.

  2. (A2)

As $n\longrightarrow\infty$, the bandwidths $h=h_{n}$ and $b=b_{n}$ satisfy
$h_{n}\longrightarrow 0$, $nh_{n}\longrightarrow\infty$, $b_{n}\longrightarrow
0$, and $\tau_{n}:=h_{n}/b_{n}\longrightarrow\tau\in[0,\infty)$.

For some results, we will impose additional restrictions on the rate that $h$
approaches zero.

The next condition restricts the uniform entropy of the class of functions in
which the nuisance estimators are assumed to reside. We briefly define uniform
entropy; we refer the reader to van der Vaart and Wellner, (1996) for
additional details. For a class of functions $\mathcal{F}$, a probability
measure $Q$ and any $\varepsilon>0$, the $\varepsilon$ covering number
$N(\varepsilon,\mathcal{F},L_{2}(Q))$ of $\mathcal{F}$ relative to the
$L_{2}(Q)$ metric is defined as the minimal number of $L_{2}(Q)$ balls of
radius less than or equal to $\varepsilon$ needed to cover $\mathcal{F}$. The
uniform $\varepsilon$-entropy of $\mathcal{F}$ is defined as $\sup_{Q}\log
N(\varepsilon,\mathcal{F},L_{2}(Q))$, where the supremum is taken over all
probability measures. We now state the following condition.

  1. (A3)

There exist classes of functions $\mathcal{F}_{\mu}$ and $\mathcal{F}_{g}$
such that almost surely for all $n$ large enough,
$\mu_{0},\mu_{n}\in\mathcal{F}_{\mu}$, $g_{0},g_{n}\in\mathcal{F}_{g}$, and
for some constants $C_{j}\in(0,\infty)$, $V_{\mu}\in(0,1)$, and
$V_{g}\in(0,2)$:

    1. (a)

$|\mu|\leq C_{1}$ for all $\mu\in\mathcal{F}_{\mu}$, and $C_{2}\leq g\leq
C_{3}$ for all $g\in\mathcal{F}_{g}$; and

    2. (b)

$\sup_{Q}\log N(\varepsilon,\mathcal{F}_{\mu},L_{2}(Q))\leq
C_{4}\varepsilon^{-V_{\mu}}$ and $\sup_{Q}\log
N(\varepsilon,\mathcal{F}_{g},L_{2}(Q))\leq C_{5}\varepsilon^{-V_{g}}$ for all
$\varepsilon>0$.

Next, we control the behavior of limiting functions to which nuisance
estimators converge. We define the following pseudo-distance for any
$P_{0}$-square integrable functions
$\gamma_{1},\gamma_{2}:\mathcal{A}\times\mathcal{W}\to\mathbb{R}$,
$\mathcal{A}_{0}\subset\mathcal{A}$, and
$\mathcal{S}\subseteq\mathcal{A}\times\mathcal{W}$:

$\displaystyle
d(\gamma_{1},\gamma_{2};\mathcal{A}_{0},\mathcal{S}):=\sup_{a\in\mathcal{A}_{0}}\left\\{E_{0}[I_{\mathcal{S}}(a,W)\\{\gamma_{1}(a,W)-\gamma_{2}(a,W)\\}^{2}]\right\\}^{1/2}$

We also define $B_{\varepsilon}(a_{0})$ as the closed ball of radius
$\varepsilon$ centered at $a_{0}$. We then state the final two conditions
concerning the rate of convergence of the nuisance estimators and properties
of the true distribution $P_{0}$. These conditions are specific to a value
$a_{0}$ because they will be used in the pointwise result.

  1. (A4)

There exist $\mu_{\infty}\in\mathcal{F}_{\mu}$ and
$g_{\infty}\in\mathcal{F}_{g}$, $\delta_{1}>0$, and subsets
$\mathcal{S}_{1},\mathcal{S}_{2}$ and $\mathcal{S}_{3}$ of
$B_{\delta_{1}}(a_{0})\times\mathcal{W}$ such that
$\mathcal{S}_{1}\cup\mathcal{S}_{2}\cup\mathcal{S}_{3}=B_{\delta_{1}}(a_{0})\times\mathcal{W}$
and:

    1. (a)

$\mu_{\infty}(a,w)=\mu_{0}(a,w)$ for all
$(a,w)\in\mathcal{S}_{1}\cup\mathcal{S}_{3}$ and $g_{\infty}(a,w)=g_{0}(a,w)$
for all $(a,w)\in\mathcal{S}_{2}\cup\mathcal{S}_{3}$;

    2. (b)

$d(\mu_{n},\mu_{\infty};B_{\delta_{1}}(a_{0}),\mathcal{S}_{1})=o_{p}\left(\\{nh\\}^{-1/2}\right)$,
and $d(g_{n},g_{\infty};B_{\delta_{1}}(a_{0}),\mathcal{S}_{1})=o_{p}(1)$;

    3. (c)

$d(g_{n},g_{\infty};B_{\delta_{1}}(a_{0}),\mathcal{S}_{2})\\}=o_{p}\left(\\{nh\\}^{-1/2}\right)$,
and $d(\mu_{n},\mu_{\infty};B_{\delta_{1}}(a_{0}),\mathcal{S}_{1})=o_{p}(1)$;
and

    4. (d)

$d(\mu_{n},\mu_{\infty};B_{\delta_{1}}(a_{0}),\mathcal{S}_{3})d(g_{n},g_{\infty};B_{\delta_{1}}(a_{0}),\mathcal{S}_{3})=o_{p}\left(\\{nh\\}^{-1/2}\right)$.

  2. (A5)

It holds that:

    1. (a)

$\theta_{0}$ is twice continuously differentiable on $B_{\delta_{1}}(a_{0})$;

    2. (b)

$f_{0}$ is positive and Lipschitz continuous on $B_{\delta_{1}}(a_{0})$;

    3. (c)

there exist $\delta_{2}>0$ and $C_{6}<\infty$ such that
$E_{0}[|Y|^{2+\delta_{2}}\mid A=a,W=w]\leq C_{6}$ for all $a\in
B_{\delta_{1}}(a_{0})$ and $P_{0}$-almost every $w$ and
$E_{0}[|Y|^{4}]<\infty$; and

    4. (d)

$a\mapsto\sigma_{0}^{2}(a)\coloneqq
E_{0}\left[\\{\xi_{\infty}(Y,A,W)-\theta_{0}(A)\\}^{2}\mid A=a\right]$ is
bounded and continuous on $B_{\delta_{1}}(a_{0})$, where
$\xi_{\infty}\coloneqq\xi_{\mu_{\infty},g_{\infty},Q_{0}}$ is the limiting
pseudo-outcome.

Finally, we define $\phi_{\infty,h,b,a_{0}}^{*}$ as the limiting efficient
influence function,

$\displaystyle\phi_{\infty,h,b,a_{0}}^{*}(y,a,w)$
$\displaystyle:=\Gamma_{0,h,b,a_{0}}(a)\xi_{\infty}(y,a,w)-\gamma_{0,h,b,a_{0}}(a)$
$\displaystyle\qquad+\int\Gamma_{0,h,b,a_{0}}(\bar{a})\left\\{\mu_{\infty}(\bar{a},w)-\int\mu_{\infty}(\bar{a},\bar{w})\,dQ_{0}(\bar{w})\right\\}\,dF_{0}(\bar{a})$

and we define $\sigma_{n,h,b}^{2}(a_{0})\coloneqq
h\mathbb{P}_{n}(\phi_{n,h,b,a_{0}}^{*})^{2}$ as an estimator of the variance
of $\theta_{n,h,b}(a_{0})$. Under the five conditions defined above, we have
the following result concerning the pointwise asymptotics of our estimator.

######  Theorem 1.

If (A1)–(A5) hold, then
$\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})=\mathbb{P}_{n}\phi_{\infty,h,b,a_{0}}^{*}+o_{p}\left(\\{nh\\}^{-1/2}+h^{2}\right),$
and
$(nh)^{1/2}\mathbb{P}_{n}\phi_{\infty,h,b,a_{0}}^{*}\operatorname*{\stackrel{{\scriptstyle
d}}{{\longrightarrow}}}N\left(0,V_{K,\tau}f_{0}(a_{0})^{-1}\sigma_{0}^{2}(a_{0})\right)$,
where

$\displaystyle V_{K,\tau}$
$\displaystyle=\int\left\\{K(u)-\tau^{3}c_{2}\frac{(\tau
u)^{2}-c_{2}}{c_{4}-c_{2}^{2}}K(\tau u)\right\\}^{2}\,du$
$\displaystyle=c_{0}^{*}-2\tau^{3}c_{2}\frac{\tau^{2}c_{2,\tau}^{*}-c_{2}c_{0,\tau}^{*}}{c_{4}-c_{2}^{2}}+\tau^{5}c_{2}^{2}\frac{c_{4}^{*}-2c_{2}c_{2}^{*}+c_{2}^{2}c_{0}^{*}}{\left(c_{4}-c_{2}\right)^{2}}$

for $c_{j}\coloneqq\int u^{j}K(u)\,du$, $c^{*}_{j}\coloneqq\int
u^{j}K^{2}(u)\,du$ and $c^{*}_{j,\tau}\coloneqq\int u^{j}K(u)K(\tau u)\,du$.
Hence, if $nh^{5}=O(1)$, then
$(nh)^{1/2}\left[\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})\right]$ converges in
distribution to this same limit. Furthermore,
$(nh)^{1/2}\left[\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})\right]/\sigma_{n,h,b}(a_{0})$
converges in distribution to $N(0,1)$.

The two most crucial features of Theorem 1 are that the estimator is centered
around $\theta_{0}(a_{0})$, and that the conditions permit the bandwidth to be
selected at the optimal rate for estimation. In particular, the final
statement of Theorem 1 implies that the confidence interval defined in (7) has
asymptotic coverage $1-\alpha$ of $\theta_{0}(a_{0})$ without undersmoothing.
The first statement of Theorem 1 resembles asymptotic linearity but differs in
that the influence function $\phi_{\infty,h,b,a_{0}}^{*}$ changes with $n$
through $h$ and $b$. Nevertheless, the result is useful for suggesting a
natural variance estimator and in revealing the process driving the first-
order behavior of the estimator.

The limiting variance of our estimator is the same as that of the local linear
estimator proposed by Kennedy et al., (2016) up to the constant $V_{K,\tau}$.
When $\tau=0$, so that $h/b\longrightarrow 0$, $V_{k,0}=\int K^{2}=c_{2}^{*}$,
which is the same as the constant in Kennedy et al., (2016). Hence, in this
case, the bias correction has no impact on the asymptotic variance of the
estimator. When $\tau>0$, the asymptotic variance of our estimator is a
constant factor larger than that of the local linear estimator. When $\tau=1$,
the constant simplifies to

$V_{K,1}=\frac{c_{0}^{*}c_{4}^{2}-2c_{2}c_{4}c_{2}^{*}+c_{2}^{2}c_{4}^{*}}{(c_{4}-c_{2}^{2})^{2}}.$

Therefore, if $\tau>0$, our debiased estimator asymptotically reduces bias at
the expense of variance. For the Epanechnikov kernel, $V_{K,1}=1.25$, while
$V_{K,0}=0.6$. Hence, debiasing approximately doubles the asymptotic variance
in this case. However, our variance estimator is consistent even when the bias
estimator contributes to the asymptotic variance.

We note that we can decompose the linear term as

$\displaystyle\mathbb{P}_{n}\phi_{\infty,h,b,a_{0}}^{*}$
$\displaystyle=\mathbb{P}_{n}\left\\{\Gamma_{0,h,b,a_{0}}(\xi_{\infty}-\theta_{0})\right\\}+\mathbb{P}_{n}\left\\{\Gamma_{0,h,b,a_{0}}\theta_{0}-\gamma_{0,h,b,a_{0}}\right\\}$
$\displaystyle\qquad+\mathbb{P}_{n}\left\\{\int\Gamma_{0,h,b,a_{0}}\left(\mu_{\infty}-\int\mu_{\infty}\,dQ_{0}\right)\,dF_{0}\right\\}.$

In the proof of Theorem 1, we show that the second and third terms in the
above decomposition are $o_{p}(\\{n/h\\}^{-1/2})$ and $O_{p}(n^{-1/2})$,
respectively. Both of these are $o_{p}(\\{nh\\}^{-1/2})$, so under the
conditions of Theorem 1, the simpler representation
$\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})=\mathbb{P}_{n}\left\\{\Gamma_{0,h,b,a_{0}}(\xi_{\infty}-\theta_{0})\right\\}+o_{p}(\\{nh\\}^{-1/2})$
holds, and therefore only the first component of $\phi_{\infty,h,b,a_{0}}^{*}$
contributes to the limit distribution of the estimator. This suggests
$h\mathbb{P}_{n}\left\\{\Gamma_{n,h,b,a_{0}}(\xi_{n}-\theta_{n})\right\\}^{2}$
as an alternative variance estimator. While this variance estimator would
still yield asymptotically valid confidence intervals, including the
additional asymptotically negligible terms in the variance estimator better
captures the finite-sample behavior of the estimator.

We now discuss the conditions used in Theorem 1. Condition (A1) is a standard
condition for kernel smoothing and is satisfied for many common kernel
functions. The bounded support condition is technically convenient but may be
avoidable. The subgraph requirement of the condition is only used for uniform
inference but is relatively mild. We impose this condition so that the class
of functions $\left\\{a\mapsto
K\left(\frac{a-a_{0}}{h}\right):a_{0}\in\mathbb{R}\right\\}$ is of VC-type
(Giné and Guillou,, 2002). The condition is satisfied in particular if $K$ is
of the form $\phi\circ p$, where $p$ is a polynomial and $\phi$ is a bounded
real function of bounded variation, which is the case for many standard
kernels including the triangular, Epanechnikov, and truncated Gaussian
kernels.

The requirements on $h$ in condition (A2) are standard in kernel smoothing.
They require that the bandwidth goes to zero, so that the estimator properly
localizes around $a_{0}$, but that it goes to zero slower than $n^{-1}$ so
that the estimator does not localize too much. In order to ensure that
$o(h^{2})=o(\\{nh\\}^{-1/2})$, the convergence in distribution part of the
result also requires that $nh^{5}=O(1)$, which means that $h$ goes to zero at
least as fast as $n^{-1/5}$. This permits but does not require undersmoothing.
The second part of condition (A2) requires that the bandwidth $b$ used for
estimating the second derivative in the bias correction goes to zero, but that
it not go to zero faster than $h$.

Condition (A3) requires that the nuisance estimators be contained in uniformly
bounded function classes, and in the case of $g_{n}$, that the function class
be uniformly bounded away from zero as well. Furthermore, (A3)(b) restricts
the uniform entropy of these function classes. The uniform entropy condition
for $\mathcal{F}_{g}$ is standard in empirical process theory since it
guarantees that the uniform entropy integral is finite. The uniform entropy
condition for $\mathcal{F}_{\mu}$ is slightly stronger in order to control
empirical U-processes associated with the integrated term
$\int\mu_{n}(a,w)\,dQ_{n}(w)$ in the estimated pseudo-outcomes.

Cross-fitting could be used to avoid the entropy condition (A3)(b) (van der
Laan et al.,, 2018; Westling et al.,, 2020; Semenova and Chernozhukov,, 2021;
Colangelo and Lee,, 2020). However, cross-fitting carries a higher
computational cost if the nuisance estimators are estimated for each training
set, and the nuisance estimators use a smaller training set. Hence, it is of
interest to determine whether the results can be obtained under entropy
conditions. Furthermore, the conditions for both classes notably do not
restrict the nuisance estimator to VC classes; hence, our conditions permit
large function classes typically associated with data-adaptive estimators.

Condition (A4) is a doubly-robust condition similar to, but slightly more
flexible than, that required by Kennedy et al., (2016). It requires that at
least one of $\mu_{n}$ and $g_{n}$ be consistent for $\mu_{0}$ or $g_{0}$,
respectively, in a neighborhood of $a_{0}$ and for almost all covariate
values. It also requires that the product of the rates of convergence of
$\mu_{n}-\mu_{0}$ and $g_{n}-g_{0}$ be faster than $(nh)^{-1/2}$ in order to
ensure negligibility of a second-order remainder term. For points at which
only one of the nuisance estimators is consistent, that estimator must achieve
this rate alone. If the covariates are low-dimensional, these rates can be
guaranteed by many nonparametric estimators. For moderate or high-dimensional
covariates, the nuisance estimators need to take advantage of additional
smoothness or structure of the true nuisance parameters to ensure these rates
of convergence are attainable (Bonvini and Kennedy,, 2022). In practice, we
recommend leveraging multiple candidate estimators in an ensemble estimator
such as SuperLearner (van der Laan et al.,, 2007).

Finally, condition (A5) imposes smoothness conditions on features of the true
distribution. Most importantly, (A5)(a) requires that $\theta_{0}$ be twice
continuously differentiable in a neighborhood of $a_{0}$. Hence, as in
Calonico et al., (2018), our estimator does not require additional smoothness
to yield asymptotically valid inference. Intuition for how this is possible,
despite using a third-order local polynomial estimator for the second
derivative estimator, was provided in Section 2. The bandwidth conditions
permit the estimator to obtain the optimal rate of convergence relative to its
assumed smoothness. Conditions (A5)(c) and (A5)(d) require slightly more than
two bounded moments for $Y$, but do not require that $Y$ be uniformly bounded.

###  3.2 Inference on causal effects

The pointwise results of Theorem 1 allow us to construct asymptotically valid
confidence intervals for $\theta_{0}(a_{0})$ for any $a_{0}$ for which the
conditions hold. However, in many cases, it is also of interest to draw
simultaneous inference on a finite collection of values
$\\{\theta_{0}(a_{1}),\dotsc,\theta_{0}(a_{m})\\}$. This can be used, for
instance, to construct confidence intervals for causal effects of the form
$\theta_{0}(a_{2})-\theta_{0}(a_{1})$. For this, joint convergence of the
estimator at several points is necessary. The next result demonstrates joint
convergence in distribution of our estimator at a finite collection of points.

######  Theorem 2.

If (A1)–(A5) hold for each $a_{0}$ in the finite and fixed collection
$\\{a_{1},\dotsc,a_{m}\\}$ and $nh^{5}=O(1)$, then

$(nh)^{1/2}\begin{pmatrix}\theta_{n,h,b}(a_{1})-\theta_{0}(a_{1})\\\ \vdots\\\
\theta_{n,h,b}(a_{m})-\theta_{0}(a_{m})\end{pmatrix}$

converges in distribution to a mean-zero multivariate normal distribution with
diagonal covariance matrix and variances as defined in Theorem 1.

Theorem 2 demonstrates that the estimator is asymptotically independent at any
two distinct points. This is because the covariance
$hP_{0}(\phi_{\infty,h,b,a_{1}}^{*}\phi_{\infty,h,b,a_{2}}^{*})$ converges to
zero as $h\longrightarrow 0$ for any $a_{1}\neq a_{2}$. Intuitively, this is
due to the fact that the estimator localizes around each point as the sample
size increases. Although this justifies estimating the covariance matrix using
a diagonal matrix with variance estimators $\sigma_{n,h,b}^{2}(a_{j})$ on the
diagonal, we recommend instead estimating the $(j,k)$ element of the
covariance matrix using the estimator
$h\mathbb{P}_{n}(\phi_{n,h,b,a_{j}}^{*}\phi_{n,h,b,a_{k}}^{*})$. This is
reminiscent of the influence-function based estimator of the covariance
between two or more asymptotically linear estimators. Though the true
covariance is going to zero, it is not necessarily zero in finite samples, and
we expect this estimator to better capture the finite-sample covariance. In
numerical studies, we demonstrate that this estimator can yield substantial
finite-sample improvements over an estimator that utilizes asymptotic
independence. To obtain inference on
$\nu(\theta_{0}(a_{1}),\dotsc,\theta_{0}(a_{m}))$ for a differentiable
function $\nu:\mathbb{R}^{m}\to\mathbb{R}$, we can combine this covariance
estimator with the delta method. Hence, Theorem 2 enables us to perform
asymptotically valid inference on effects of the form
$\theta_{0}(a_{1})-\theta_{0}(a_{2})$ without undersmoothing.

###  3.3 Uniform inference

We now turn to the uniform behavior of the estimator over a compact set
$\mathcal{A}_{0}\subset\mathcal{A}$. Our goal is to construct an
asymptotically valid uniform confidence band for $\theta_{0}$ over
$\mathcal{A}_{0}$. A standard approach to this problem would be to demonstrate
that
$\\{(nh)^{1/2}[\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})]:a_{0}\in\mathcal{A}_{0}\\}$
converges weakly as a process to a tight limit in the space
$\ell^{\infty}(\mathcal{A}_{0})$ of uniformly bounded functions on
$\mathcal{A}_{0}$ equipped with the supremum norm. However, by the asymptotic
independence of $(nh)^{1/2}[\theta_{n,h,b}(a_{1})-\theta_{0}(a_{1})]$ and
$(nh)^{1/2}[\theta_{n,h,b}(a_{2})-\theta_{0}(a_{2})]$ demonstrated in Theorem
2, the only possible limit is a white noise process, which is not tight in
$\ell^{\infty}(\mathcal{A}_{0})$. Stupfler et al., (2016) explored this
phenomenon in more depth for kernel density estimators. Instead, we will use
the finite-sample approximation theory developed in Chernozhukov et al.,
(2014). Using this theory and the first-order representation
$\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})=\mathbb{P}_{n}\phi_{0,h,b,a_{0}}^{*}+o_{p}(\\{nh\\}^{-1/2})$,
we approximate the distribution of
$\sup_{\mathcal{A}_{0}}(nh)^{1/2}|\theta_{n,h,b}-\theta_{0}|/\sigma_{n,h,b}$
with that of $\max_{\mathcal{A}_{n}}|Z_{n,h,b}|$, where $\mathcal{A}_{n}$ is a
finite subset of $\mathcal{A}_{0}$ and conditional on the data, $Z_{n,h,b}$ is
a multivariate Gaussian random vector on $\mathcal{A}_{n}$ with covariance
given by

$\mathrm{Cov}\left(Z_{n,h,b}(a_{1}),Z_{n,h,b}(a_{2})\right)=h\mathbb{P}_{n}\left(\phi_{n,h,b,a_{1}}^{*}\phi_{n,h,b,a_{2}}^{*}\right)/\left[\sigma_{n,h,b}(a_{1})\sigma_{n,h,b}(a_{2})\right].$

Letting $t_{1-\alpha,n}$ represent the $1-\alpha$ quantile of the conditional
distribution of $\max_{\mathcal{A}_{n}}|Z_{n,h,b}|$ given the data, the lower
and upper limits of our asymptotic $(1-\alpha)$ confidence band for
$\theta_{0}$ over $a_{0}\in\mathcal{A}_{0}$ are then given by
$\theta_{n,h,b}(a_{0})\pm(nh)^{-1/2}t_{1-\alpha,n}\sigma_{n,h,b}(a_{0})$.
Notably, the limits of this confidence band are proportional to the limits of
the confidence interval defined in (7).

We now introduce additional conditions we will use. For $\delta>0$, we define
$\mathcal{A}_{\delta}$ as the $\delta$-enlargement of $\mathcal{A}_{0}$, that
is, the set of $a\in\mathbb{R}$ such that there exists
$a_{0}\in\mathcal{A}_{0}$ with $|a-a_{0}|\leq\delta$.

  1. (A6)

The constant $V:=\max\\{V_{\mu},V_{g}\\}$, for $V_{\mu}$ and $V_{g}$ defined
in (A3) and the bandwidth $h$ satisfy $n\left[h/(\log
n)\right]^{\frac{2+V}{2-V}}\longrightarrow\infty$ and
$nh^{3}\longrightarrow\infty$.

  2. (A7)

There exist $\mu_{\infty}\in\mathcal{F}_{\mu}$,
$g_{\infty}\in\mathcal{F}_{g}$, $\delta_{3}>0$ and subsets
$\mathcal{S}_{1}^{\prime},\mathcal{S}_{2}^{\prime}$, and
$\mathcal{S}_{3}^{\prime}$ of $\mathcal{A}_{\delta_{3}}\times\mathcal{W}$ such
that
$\mathcal{S}_{1}^{\prime}\cup\mathcal{S}_{2}^{\prime}\cup\mathcal{S}_{3}^{\prime}=\mathcal{A}_{\delta_{3}}\times\mathcal{W}$
and:

    1. (a)

$\mu_{\infty}(a,w)=\mu_{0}(a,w)$ for all
$(a,w)\in\mathcal{S}_{1}^{\prime}\cup\mathcal{S}_{3}^{\prime}$ and
$g_{\infty}(a,w)=g_{0}(a,w)$ for all
$(a,w)\in\mathcal{S}_{2}^{\prime}\cup\mathcal{S}_{3}^{\prime}$;

    2. (b)

$d(\mu_{n},\mu_{\infty};\mathcal{A}_{\delta_{3}},\mathcal{S}_{1}^{\prime})=o_{p}\left(\\{nh\log
n\\}^{-1/2}\right)$;

    3. (c)

$d(g_{n},g_{\infty};\mathcal{A}_{\delta_{3}},\mathcal{S}_{2}^{\prime})\\}=o_{p}\left(\\{nh\log
n\\}^{-1/2}\right)$;

    4. (d)

$d(\mu_{n},\mu_{\infty};\mathcal{A}_{\delta_{3}},\mathcal{S}_{3})d(g_{n},g_{\infty};\mathcal{A}_{\delta_{3}},\mathcal{S}_{3})=o_{p}\left(\\{nh\log
n\\}^{-1/2}\right)$; and

    5. (e)

$d(\mu_{n},\mu_{\infty};\mathcal{A}_{\delta_{3}},\mathcal{A}\times\mathcal{W})$
and
$d(g_{n},g_{\infty};\mathcal{A}_{\delta_{3}},\mathcal{A}\times\mathcal{W})$
are both $o_{p}\left(h^{\frac{V}{2(2-V)}}\\{\log
n\\}^{-\frac{1}{2-V}}\right)$.

  3. (A8)

It holds that:

    1. (a)

$\theta_{0}$ is twice continuously differentiable with Hölder-continuous
second derivative on $\mathcal{A}_{\delta_{3}}$;

    2. (b)

$f_{0}$ is Lipschitz continuous and bounded away from 0 and $\infty$ on
$\mathcal{A}_{\delta_{3}}$;

    3. (c)

$|Y|$ is $P_{0}$-almost surely bounded; and

    4. (d)

$a\mapsto E_{0}\left[\\{\xi_{\infty}(Y,A,W)-\theta_{0}(A)\\}^{2}\mid
A=a\right]$ is continuous on $\mathcal{A}_{\delta_{3}}$.

We now present a result demonstrating asymptotic validity of our proposed
uniform confidence band. We define $\mathbf{O}_{n}:=(O_{1},\dotsc,O_{n})$ and
$\omega_{n}:=\sup_{a_{0}\in\mathcal{A}_{0}}\min_{a\in\mathcal{A}_{n}}|a-a_{0}|$
as the mesh of $\mathcal{A}_{n}$ in $\mathcal{A}_{0}$.

######  Theorem 3.

If (A1)–(A3) and (A6)–(A8) hold, $nh^{5}=O(1)$,
$\omega_{n}=o\left(h^{p}\right)$ for some $p>1$, and $mn^{d}=O(1)$ for some
$d\in(0,\infty)$, then

$\displaystyle\sup_{t\in\mathbb{R}}\left|P_{0}\left(\sup_{a_{0}\in\mathcal{A}_{0}}\\{nh\\}^{1/2}\left|\frac{\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})}{\sigma_{n,h,b}(a_{0})}\right|\leq
t\right)-P_{0}\left(\max_{a_{0}\in\mathcal{A}_{n}}\left|Z_{n,h,b}(a_{0})\right|\leq
t\mid\mathbf{O}_{n}\right)\right|=o_{p}(1).$

Theorem 3 implies that our proposed uniform confidence band has asymptotically
valid coverage. This is because by definition,
$P_{0}\left(\max_{a_{0}\in\mathcal{A}_{n}}\left|Z_{n,h,b}(a_{0})\right|\leq
t_{1-\alpha,n}\mid\mathbf{O}_{n}\right)=1-\alpha$, and combined with Theorem 3
this implies that

$P_{0}\left(\sup_{a_{0}\in\mathcal{A}_{0}}\\{nh\\}^{1/2}\left|\frac{\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})}{\sigma_{n,h,b}(a_{0})}\right|\leq
t_{1-\alpha,n}\right)\longrightarrow 1-\alpha.$

The expression inside the probability is equivalent to the statement that
$\theta_{0}(a_{0})$ is contained in the proposed confidence bands for all
$a_{0}\in\mathcal{A}_{0}$.

To prove Theorem 3, we first use the results of Chernozhukov et al., (2014) to
demonstrate that the distribution of
$\sup_{a_{0}\in\mathcal{A}_{0}}\left|\mathbb{G}_{n}h^{1/2}\phi^{*}_{\infty,h,b,a_{0}}/\sigma_{\infty,h,b}(a_{0})\right|$
can be approximated by that of $\sup_{\mathcal{A}_{0}}|Z_{\infty,h,b}|$, where
$Z_{\infty,h,b}$ is a mean-zero Gaussian process on $\mathcal{A}_{0}$ with
covariance

$\mathrm{Cov}(Z_{\infty,h,b}(a_{1}),Z_{\infty,h,b}(a_{2})):=hP_{0}(\phi^{*}_{\infty,h,b,a_{1}}\phi^{*}_{\infty,h,b,a_{2}})/[\sigma_{\infty,h,b}(a_{1})\sigma_{\infty,h,b}(a_{2})],$

where
$\sigma_{\infty,h,b}^{2}(a_{0}):=hP_{0}(\phi^{*}_{\infty,h,b,a_{0}})^{2}$ and
$\mathbb{G}_{n}:=n^{1/2}(\mathbb{P}_{n}-P_{0})$. Notably, the results of
Chernozhukov et al., (2014) cover situations where the empirical process is
not converging weakly to a tight limit, which is the case for our process.
Since we also establish that
$\sup_{a_{0}\in\mathcal{A}_{0}}|\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})-\mathbb{P}_{n}\phi_{\infty,h,b,a_{0}}^{*}|=o_{p}\left(\\{nh\log
n\\}^{-1/2}\right)$, this means that the distribution of
$(nh)^{1/2}\sup_{\mathcal{A}_{0}}|\theta_{n,h,b}-\theta_{0}|$ can also be
approximated by that of $\sup_{\mathcal{A}_{0}}|Z_{\infty,h,b}|$. Finally, we
use the results of Chernozhukov et al., (2015) to demonstrate that the
distribution of $\sup_{\mathcal{A}_{0}}|Z_{\infty,h,b}|$ can be approximated
by that of $\max_{\mathcal{A}_{n}}|Z_{n,h,b}|$.

In order to show that
$\sup_{a_{0}\in\mathcal{A}_{0}}|\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})-\mathbb{P}_{n}\phi_{\infty,h,b,a_{0}}^{*}|=o_{p}\left(\\{nh\log
n\\}^{-1/2}\right)$, we need faster rates of convergence for the remainder
terms in the first-order expansion of our estimator. For some remainder terms,
this is straightforward. However, for the empirical process remainder term
$(\mathbb{P}_{n}-P_{0})(\phi_{n,h,b,a_{0}}^{*}-\phi_{\infty,h,b,a_{0}}^{*})$,
it is a challenging task because the usual method of demonstrating
negligibility of this remainder, namely uniform asymptotic equicontinuity,
only yields that it is $o_{p}\left(\\{nh\\}^{-1/2}\right)$. In order to
achieve the extra $\\{\log n\\}^{-1/2}$ term in the rate, we use the local
maximal inequalities of van der Vaart and Wellner, (2011). By assuming rates
of convergence for $\mu_{n}-\mu_{\infty}$ and $g_{n}-g_{\infty}$ in condition
(A7)(e), we are able to establish a rate of convergence for
$P_{0}(\phi_{n,h,b,a_{0}}^{*}-\phi_{\infty,h,b,a_{0}}^{*})^{2}$, which then
permits the use of the results of van der Vaart and Wellner, (2011).

We now discuss the additional conditions used in Theorem 3 beyond those used
in Theorem 1. Condition (A6)(a) further restricts the exponent of the uniform
entropy of the nuisance classes, and is used to control the empirical process
remainder using local maximal inequalities as discussed above. If $h$ is
chosen at the optimal rate $n^{-1/5}$, then the condition is satisfied if
$V_{\mu},V_{g}\in(0,4/3)$. If undersmoothing is employed, then the requirement
is stricter. The requirement that $nh^{3}\longrightarrow\infty$ can be relaxed
somewhat, especially if $V_{\mu}$ is much less than 1, but not beyond
$nh^{2}\longrightarrow\infty$. Hence, severe undersmoothing in addition to
debiasing is possible without sample-splitting if the nuisance estimators fall
in smaller classes, though we expect this is of less interest because a main
point of debiasing is to avoid undersmoothing.

Condition (A7) is a uniform version of the doubly-robust condition (A4). The
rates of convergence of the nuisance estimators on the sets on which they are
consistent in parts (b), (c), and (d) are slightly faster in order to
establish the needed rate of convergence of the second-order remainder term.
Part (e) also requires a rate of convergence of the nuisance estimators
towards their limiting objects. As discussed above, this is also to enable the
use of local maximal inequalities. If $h$ is chosen at the optimal rate
$n^{-1/5}$ and $V_{g}<4/3$ as required by (A6)(a), then part (e) is satisfied
if $d(\mu_{n},\mu_{\infty};A_{\delta_{3}},\mathcal{A}\times\mathcal{W})$ and
$d(g_{n},g_{\infty};A_{\delta_{3}},\mathcal{A}\times\mathcal{W})$ are each
$o_{p}(n^{-1/5})$.

Condition (A8) places conditions on the true distribution that are stronger
than those required in condition (A5) for pointwise convergence. The assumed
smoothness of $\theta_{0}$ and $f_{0}$ is standard in the literature (see,
e.g., Assumption 3 in Calonico et al.,, 2018 and R1 in Cheng and Chen,, 2019),
and we notably still do not require more than two derivatives for
$\theta_{0}$. Part (c) requires $Y$ to be uniformly bounded, which is again
used to apply the local maximal inequality. However, this could be relaxed to
$|Y|$ having an almost surely finite conditional $q$th moment for $q>4$ at the
expense of stronger and more complicated restrictions on the complexity and
rate of convergence of the nuisance estimators.

Theorem 3 also requires that the mesh $\omega_{n}$ of $\mathcal{A}_{n}$
decrease faster than $h^{p}$ for some $p>1$ but that the number of points $m$
in $\mathcal{A}_{n}$ increase at most at a polynomial rate. The former is to
ensure that the distribution of $\sup_{\mathcal{A}_{n}}|Z_{\infty,h,b}|$ is a
good enough approximation to that of $\sup_{\mathcal{A}_{0}}|Z_{\infty,h,b}|$,
while the latter is to ensure that the distribution of
$\max_{\mathcal{A}_{n}}|Z_{n,h,b}|$ is a good enough approximation of
$\max_{\mathcal{A}_{n}}|Z_{\infty,h,b}|$. These conditions can be
simultaneously achieved, for instance, with a uniform grid of $n$ points on
$\mathcal{A}_{0}$.

##  4 Numerical studies

###  4.1 Study design

In this section, we conduct numerical studies to investigate the finite-sample
behavior of the proposed methods. We begin by describing our data-generating
process. First, we generate covariates $W\in\mathbb{R}^{4}$ from a standard
multivariate normal distribution. Given $W$, we then generate $A$ from the
distribution whose conditional density function is given by $p_{0}(a\mid
w)\coloneqq I_{[0,1]}(a)[\lambda(w)+2a\\{1-\lambda(w)\\}]$ for
$\lambda(w)\coloneqq 0.1+1.8\,\mathrm{expit}(\beta^{T}w)$, where
$\mathrm{expit}(x):=1/(1+e^{-x})$. This construction guarantees that
$0.1<p_{0}(a\mid w)<1.9$ for all $a\in[0,1]$ and $w\in\mathbb{R}^{4}$ and that
$A$ is marginally $\mathrm{Uniform}(0,1)$. Finally, we generate $Y$ given
$A=a$ and $W=w$ as a Bernoulli random variable with mean
$\mu_{0}(a,w)\coloneqq\text{expit}\left(\gamma_{1}^{T}\bar{w}+\gamma_{2}^{T}\bar{w}a+\gamma_{3}a^{2}+\gamma_{4}T(a)\right),$
where $T(a)\coloneqq\sin(3\pi\\{2a-1\\}/2)/(a^{2}+1)$ and $\bar{w}:=(1,w)$. We
set $\beta=(-1,-1,1,1)^{T}$, $\gamma_{1}=(-1,-1,-1,1,1)^{T}$,
$\gamma_{2}=(3,-1,-1,1,1)^{T}$, $\gamma_{3}=1$ and $\gamma_{4}=3$. Figure 1
displays $\theta_{0}$ and its second derivative $\theta_{0}^{(2)}$ implied by
these settings. The resulting curve is non-monotonic and has large second
derivatives for some regions of $[0,1]$.

Figure 1: The true function $\theta_{0}$ and its second derivative
$\theta_{0}^{(2)}$ used in the numerical experiements.

We simulated 1000 datasets using the above process for each
$n\in\\{500,1000,2500\\}$. For each dataset, we constructed nuisance
estimators $\mu_{n}$ and $g_{n}$ in the following ways. First, we estimated
$\mu_{0}$ using a correctly specified logistic regression model to obtain
estimators of $\gamma_{1}$, $\gamma_{2}$, $\gamma_{3}$ and $\gamma_{4}$. We
estimated $g_{0}$ using maximum likelihood estimation with a correctly
specified parametric model to obtain an estimator of $\beta$. To investigate
the double-robustness of our estimator, we also estimated $\mu_{0}$ and
$g_{0}$ using incorrectly specified parametric models in which the covariates
$W_{3}$ and $W_{4}$ and the corresponding interaction terms were omitted. We
then considered three settings: (i) both $\mu_{n}$ and $g_{n}$ are based on
correctly specified parametric models, (ii) $\mu_{n}$ is based on a correctly
specified model, while $g_{n}$ is based on an incorrectly specified model, and
(iii) $g_{n}$ is based on a correctly specified model, while $\mu_{n}$ is
based on an incorrectly specified model.

We estimated $\theta_{0}$ using both the local linear estimator of Kennedy et
al., (2016) and our debiased estimator. For the local linear estimator, we
considered two bandwidth selection mechanisms. First, we used the leave-one-
out cross-validation bandwidth selection proposed in Section 3.5 of Kennedy et
al., (2016), which we refer to as “Local linear (LOOCV)”. Second, we used the
plug-in methodology described in Section 2.5, which we refer to as “Local
linear (Plug-in)”. For the debiased estimator, we used all three selection
procedures described in Section 2.5. We refer to the estimator with LOOCV
bandwidths $(h_{\mathrm{cv}},b_{\mathrm{cv}})$ as “Debiased (LOOCV)”, the
estimator with LOOCV bandwidths $(h_{\mathrm{cv},1},h_{\mathrm{cv},1})$
selected by minimizing IMSEcv over $h$ alone with $b=h$ as “Debiased (LOOCV,
h=b)”, and the estimator with bandwidths $(h_{\mathrm{plug-
in}},h_{\mathrm{plug-in}})$ based on the plug-in methodology as “Debiased
(Plug-in)”. We constructed 95% confidence intervals for each
$a_{0}\in\\{0.1,0.15,0.2,\dots,0.9\\}$. We also constructed 95% confidence
intervals for $\theta_{0}(a_{0})-\theta_{0}(0.5)$ for
$a_{0}\in\\{0.51,0.52,\dotsc,0.9\\}$ using our estimator and both variance
estimators described following Theorem 2. Finally, we constructed a 95%
confidence band over $\mathcal{A}_{0}=[0.1,0.9]$ using the method described in
Section 3.3.

###  4.2 Results of numerical studies

For ease of viewing, we focus here on the results for $n=1000$. Results for
$n=500$ and $n=2500$ are provided in Supplementary Material.

Figure 2 displays the pointwise empirical bias and variance of the two
methods. As expected, the magnitude of the bias of both local linear
estimators is generally larger than that of all three debiased estimators. In
particular, the local linear estimator has a large bias where the second
derivative of $\theta_{0}$ is large. The bias of the estimators using the
plug-in bandwidth selection is generally larger than the other two methods,
especially for the local linear estimator, because the plug-in bandwidths are
on average larger than those based on LOOCV (see Figure 11 in Supplementary
Material).

The variance of the local linear estimator is smaller than that of the
debiased estimator when using the same bandwidth selection procedure. In
particular, the variance of the local linear estimator is about one-half that
of the debiased estimator with the cross-validated bandwidth and $\tau=1$
fixed, which agrees with the constants computed after Theorem 1. The variance
of the debiased estimator with the cross-validated bandwidth optimized over
both $h$ and $b$ is larger than the same optimized over just $h$. This is
somewhat surprising since we would expect optimizing over $b$ as well to
reduce variance. This is due to the large variance of the bandwidths selected
by this method (see Figure 11 in Supplementary Material). The variance of both
estimators using plug-in bandwidth selection is smaller than the same using
LOOCV bandwidth selection, which again is due to the plug-in bandwidths being
larger on average. The variance of all estimators is larger when the outcome
regression model is misspecified.

Figure 2: Top: squared empirical bias of the estimators on the log scale. The
local linear estimator (shown in red) has large bias at points where the
magnitude of the second derivative is large. Bottom: empirical variance of the
estimators on the log scale.

Figure 3 displays the empirical coverage of pointwise 95% confidence intervals
for $\theta_{0}(a_{0})$. The asymptotic bias of the local linear estimator
results in poor coverage at points where the second derivative is large in
magnitude. The debiased estimator has good coverage for all bandwidth
selection procedures considered. The plug-in method has the best coverage at
interior points, but over-covers near the boundaries.

Figure 3: Empirical coverage of 95% pointwise confidence intervals based on
the debiased local linear estimator and the local linear estimator.

Figure 4 displays the empirical coverage of confidence intervals for the
causal effect $\theta_{0}(a)-\theta_{0}(0.5)$ based on the debiased estimator
and the two variance estimators described following Theorem 2. The confidence
intervals based on the variance estimator using asymptotic independence (top
row) are conservative when $a$ is close to 0.5 because the finite-sample
covariance between the estimators is positive when the distance between the
evaluation points is small. When $a$ is farther from $0.5$, the confidence
intervals have better coverage. The confidence intervals using the variance
estimator based on the influence functions (bottom row) have much better
coverage for all values of $a$, especially at large sample sizes, because the
variance estimator captures some of the finite-sample covariance between the
estimators.

Figure 4: Empirical coverage of 95% pointwise confidence intervals for
$\theta_{0}(a)-\theta_{0}(0.5)$ based on the debiased estimator. The intervals
in the top row use the sum of two variance estimators; those in the bottom row
use the influence function-based variance estimator.

Finally, Figure 5 displays the empirical coverage of the uniform confidence
bands. For this simulation, we considered the augmented set of sample sizes
$n\in\\{500,750,1000,1250,1500,2000\\}$. The plug-in bandwidth selection
method exhibits slight undercoverage for sample sizes less than 1000, but
generally performs very well. However, both bandwidth selection methods based
on cross-validation have serious undercoverage at sample sizes less than 1500.
The coverage improves with sample size, but both methods still exhibit
undercoverage by sample size 2000. We conjecture that this undercoverage is a
result of the bandwidth $h$ selected using LOOCV being smaller on average than
that selected by the plug-in methodology. Smaller bandwidths result in a
process with smaller correlation between points and whose supremum is
stochastically larger. This results in a slower rate of convergence for the
approximation in Theorem 3. This further illustrates the benefit of permitting
that bandwidths be selected at the optimal rate.

Figure 5: Empirical coverage of $95\%$ uniform confidence bands based on the
debiased estimator.

##  5 Effect of fine particulate matter on cardiovascular mortality

Fine particulate matter is a common air pollutant, and its concentration in
the United States is regulated by the Environmental Protection Agency under
the Clean Air Act. Numerous scientific studies have reported an association
between exposure to air pollution and adverse health outcomes. The reduction
in the concentration of particulate matter in the atmosphere over the past
several decades is considered one of the contributors to the declining
cardiovascular mortality rate (Pope et al.,, 2002, 2009; Correia et al.,,
2013; Roth et al.,, 2017; Corrigan et al.,, 2018).

Wyatt et al., 2020a  recently conducted an observational study to investigate
the association between particulate matter less than 2.5 microns in diameter
(PM2.5) and cardiovascular mortality rate. Socioeconomic factors are a
potential confounding of this relationship because they impact both exposure
to PM2.5 and cardiovascular mortality. Using data recorded at the county level
in the United States between 1990 and 2010, the authors found a positive
association between PM2.5 and cardiovascular mortality rate after adjusting
for socioeconomic characteristics using regression models.

We used the publicly available data compiled by Wyatt et al., 2020b  to
estimate the covariate-adjusted relationship between PM2.5 and cardiovascular
mortality rate using the methods presented here. The data contains information
about $n=2132$ counties. Our exposure $A$ was the county-level annual PM2.5
(in $\mu$g/m3) averaged over twenty observations from 1990 to 2010, as
measured by US Environmental Protection Agency’s Community Multi-scale Air
Quality modeling system (Gan et al.,, 2015). Our outcome $Y$ was county-level
cardiovascular mortality rate (CMR, deaths/100,000 people) in 2010, as
measured by the National Center for Health Statistics. Our covariates $W$
consisted of county-level socioeconomic factors based on 1990 and 2000 census
data: population, households below the poverty level, proportion of female-
headed households with dependent children, vacant housing units, owner-
occupied housing units, median household income, high education rate,
unemployment rate, and density of medical facilities. We note that there are
very likely to be violations of the stable unit treatment value assumption in
this data since air pollution in one county could impact the health of
residents in neighboring counties. In addition, the data are unlikely to be
independent. It would be of interest in future research to extend the methods
presented here to deal with treatment spillover and dependent data.

We estimated the outcome regression using SuperLearner (van der Laan et al.,,
2007) with a library consisting of generalized linear models, multivariate
adaptive regression splines, generalized additive models, and regression
trees. We estimated the conditional density using the version of SuperLearner
developed by Díaz Muñoz and van der Laan, (2011) with the same library. We
used the Epanechnikov kernel and selected the bandwidth $h$ using the plug-in
method discussed in Section 2.5 and set $b=h$. We focus on values of PM2.5
between $3.5\mu g/m^{3}$ and $10\mu g/m^{3}$, which approximately corresponds
to the 0.05 and 0.95 quantiles of the marginal empirical distribution of
PM2.5, respectively.

Figure 6: Estimated covariate-adjusted relationship between the concentration
of fine particulate matter (PM2.5) in the atmosphere and cardiovascular
mortality rate at the county level. The regression function is adjusted for
socioeconomic factors. The figure shows 95% pointwise confidence intervals as
vertical lines and a 95% uniform confidence band as dashed lines.

Figure 6 displays the estimated covariate-adjusted regression function using
our debiased method. Pointwise 95% confidence intervals are displayed using
vertical lines, and a 95% uniform confidence band is displayed using dashed
lines. The adjusted expected CMR appears to be monotone increasing as a
function of PM2.5, and its rate of increase is fastest for PM2.5 values
between roughly 6 and 8 $\mu g/m^{3}$. Our method shows larger uncertainty for
PM2.5 levels below 4 $\mu g/m^{3}$ and above 9 $\mu g/m^{3}$ due to the low
density of PM2.5. It is unclear if average CMR increases or plateaus over
these regions. We estimate that increasing PM2.5 from 6 to 8 $\mu g/m^{3}$
corresponds to an increase of 11.62 in expected CMR after adjusting for
socioeconomic factors (95% CI: 8.39–14.85), and increasing PM2.5 from 5 to 9
$\mu g/m^{3}$ corresponds to an increase of 17.48 in expected CMR (95% CI:
13.55–21.41). Our conclusions generally agree with those reported in Wyatt et
al., 2020a , but our flexible approach to estimation and principled approach
to inference ensure that the conclusions are not the result of statistical
bias.

##  6 Concluding remarks

In this article, we studied nonparametric inference for a covariate-adjusted
regression function. This problem has wide applications in observational
studies when the exposure of interest is continuous. In particular, under
appropriate causal assumptions, the covariate-adjusted regression function
corresponds to the average counterfactual outcome had all units been assigned
to the same fixed exposure level. We presented conditions under which our
proposed procedures yield valid pointwise and uniform inference. Our
conditions do not require undersmoothing and permit the use of data-adaptive
estimators for nuisance functions, and our results do not require more
smoothness conditions than the original local linear estimator.

There are several natural extensions to our work. Cross-fitting the nuisance
estimators would remove the empirical process conditions. Debiasing higher-
order local polynomial estimators would yield faster rates of convergence
under stronger smoothness assumptions. Debiased estimators of derivatives of
the parameter could be obtained using similar methodology. Debiasing a higher-
order corrected estimator would yield weaker assumptions for the rates of
convergence of the nuisance estimators (Bonvini and Kennedy,, 2022).
Alternative procedures for targeting the smoothed parameter, such as targeted
minimum loss-based estimation (van der Laan and Rose,, 2011) could yield
improved finite-sample or asymptotic properties of the methods. Finally,
twicing kernels (Newey et al.,, 2004; Zhang and Xia,, 2012) may be an
alternative approach to bias correction.

## Acknowledgements

The authors gratefully acknowledge support from the University of
Massachusetts Amherst Department of Mathematics and Statistics startup fund
and NSF Award 2113171.

## References

  * Bickel, (1982) Bickel, P. J. (1982).  On adaptive estimation.  The Annals of Statistics, pages 647–671.
  * Bickel and Klaassen, (1993) Bickel, P. J. and Klaassen, C. A. (1993).  Efficient and adaptive estimation for semiparametric models, volume 4.  Springer.
  * Bonvini and Kennedy, (2022) Bonvini, M. and Kennedy, E. H. (2022).  Fast convergence rates for dose-response estimation.  arXiv preprint arXiv:2207.11825.
  * Calonico et al., (2018) Calonico, S., Cattaneo, M. D., and Farrell, M. H. (2018).  On the Effect of Bias Estimation on Coverage Accuracy in Nonparametric Inference.  Journal of the American Statistical Association, 113(522):767–779.
  * Cheng and Chen, (2019) Cheng, G. and Chen, Y.-C. (2019).  Nonparametric inference via bootstrapping the debiased estimator.  Electronic Journal of Statistics, 13(1):2194–2256.
  * Chernozhukov et al., (2014) Chernozhukov, V., Chetverikov, D., and Kato, K. (2014).  Gaussian approximation of suprema of empirical processes.  The Annals of Statistics, 42(4):1564–1597.
  * Chernozhukov et al., (2015) Chernozhukov, V., Chetverikov, D., and Kato, K. (2015).  Comparison and anti-concentration bounds for maxima of gaussian random vectors.  Probability Theory and Related Fields, 162(1):47–70.
  * Colangelo and Lee, (2020) Colangelo, K. and Lee, Y.-Y. (2020).  Double Debiased Machine Learning Nonparametric Inference with Continuous Treatments.  arXiv preprint arXiv:2004.03036.
  * Correia et al., (2013) Correia, A. W., Pope III, C. A., Dockery, D. W., Wang, Y., Ezzati, M., and Dominici, F. (2013).  The effect of air pollution control on life expectancy in the united states: an analysis of 545 us counties for the period 2000 to 2007.  Epidemiology, 24(1):23–31.
  * Corrigan et al., (2018) Corrigan, A. E., Becker, M. M., Neas, L. M., Cascio, W. E., and Rappold, A. G. (2018).  Fine Particulate Matters: The Impact of Air Quality Standards on Cardiovascular Mortality.  Environmental Research, 161:364–369.
  * Díaz and van der Laan, (2013) Díaz, I. and van der Laan, M. J. (2013).  Targeted Data Adaptive Estimation of the Causal Dose–Response Curve.  Journal of Causal Inference, 1(2):171–192.
  * Díaz Muñoz and van der Laan, (2011) Díaz Muñoz, I. and van der Laan, M. J. (2011).  Super learner based conditional density estimation with application to marginal structural models.  The International Journal of Biostatistics, 7(1).
  * Fan and Gijbels, (1996) Fan, J. and Gijbels, I. (1996).  Local Polynomial Modelling and Its Applications: Monographs on Statistics and Applied Probability 66 (1st ed.).  Routledge.
  * Gan et al., (2015) Gan, C.-M., Pleim, J., Mathur, R., Hogrefe, C., Long, C., Xing, J., Wong, D., Gilliam, R., and Wei, C. (2015).  Assessment of long-term WRF–CMAQ simulations for understanding direct aerosol effects on radiation” brightening” in the United States.  Atmospheric Chemistry and Physics, 15(21):12193–12209.
  * Gill and Robins, (2001) Gill, R. D. and Robins, J. M. (2001).  Causal Inference for Complex Longitudinal Data: The Continuous Case.  The Annals of Statistics, 29(6):1785–1811.
  * Giné and Guillou, (2002) Giné, E. and Guillou, A. (2002).  Rates of strong uniform consistency for multivariate kernel density estimators.  Annales de l’Institut Henri Poincare (B) Probability and Statistics, 38(6):907–921.
  * Hall, (1992) Hall, P. (1992).  Effect of Bias Estimation on Coverage Accuracy of Bootstrap Confidence Intervals for a Probability Density.  The Annals of Statistics, 20(2):675–694.
  * Hirano and Imbens, (2004) Hirano, K. and Imbens, G. W. (2004).  The Propensity Score with Continuous Treatments.  In Applied Bayesian Modeling and Causal Inference from Incomplete‐Data Perspectives, chapter 7, pages 73–84. John Wiley & Sons.
  * Imai and Van Dyk, (2004) Imai, K. and Van Dyk, D. A. (2004).  Causal inference with general treatment regimes: Generalizing the propensity score.  Journal of the American Statistical Association, 99(467):854–866.
  * Kennedy et al., (2016) Kennedy, E. H., Ma, Z., McHugh, M. D., and Small, D. S. (2016).  Non-parametric methods for doubly robust estimation of continuous treatment effects.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79(4):1229–1245.
  * Mize et al., (2019) Mize, T. D., Doan, L., and Long, J. S. (2019).  A General Framework for Comparing Predictions and Marginal Effects across Models.  Sociological Methodology, 49(1):152–189.
  * Neugebauer and van der Laan, (2007) Neugebauer, R. and van der Laan, M. (2007).  Nonparametric causal effects based on marginal structural models.  Journal of Statistical Planning and Inference, 137(2):419 – 434\.
  * Newey et al., (2004) Newey, W. K., Hsieh, F., and Robins, J. M. (2004).  Twicing Kernels and a Small Bias Property of Semiparametric Estimators.  Econometrica, 72(3):947–962.
  * Pfanzagl and Wefelmeyer, (1985) Pfanzagl, J. and Wefelmeyer, W. (1985).  Contributions to a general asymptotic statistical theory.  Statistics & Risk Modeling, 3(3-4):379–388.
  * Pollard, (1984) Pollard, D. (1984).  Convergence of Stochastic Processes. Springer.
  * Pope et al., (2002) Pope, C. A., Burnett, R. T., Thun, M. J., Calle, E. E., Krewski, D., Ito, K., and Thurston, G. D. (2002).  Lung Cancer, Cardiopulmonary Mortality, and Long-term Exposure to Fine Particulate Air Pollution.  JAMA, 287(9):1132–1141.
  * Pope et al., (2009) Pope, C. A., Ezzati, M., and Dockery, D. W. (2009).  Fine-Particulate Air Pollution and Life Expectancy in the United States.  New England Journal of Medicine, 360(4):376–386.
  * Robins, (1986) Robins, J. (1986).  A new approach to causal inference in mortality studies with a sustained exposure period—application to control of the healthy worker survivor effect.  Mathematical Modelling, 7(9-12):1393–1512.
  * Robins, (2000) Robins, J. M. (2000).  Marginal Structural Models versus Structural nested Models as Tools for Causal inference.  In Halloran, M. E. and Berry, D., editors, Statistical Models in Epidemiology, the Environment, and Clinical Trials, pages 95–133. Springer New York.
  * Roth et al., (2017) Roth, G. A., Dwyer-Lindgren, L., Bertozzi-Villa, A., Stubbs, R. W., Morozoff, C., Naghavi, M., Mokdad, A. H., and Murray, C. J. L. (2017).  Trends and Patterns of Geographic Variation in Cardiovascular Mortality Among US Counties, 1980-2014.  JAMA, 317(19):1976–1992.
  * Rubin and van der Laan, (2006) Rubin, D. and van der Laan, M. J. (2006).  Extending Marginal Structural Models through Local, Penalized, and Additive Learning.  Working Paper 212, Division of Biostatistics, University of California at Berkeley, Berkeley, California.
  * Semenova and Chernozhukov, (2021) Semenova, V. and Chernozhukov, V. (2021).  Debiased machine learning of conditional average treatment effects and other causal functions.  The Econometrics Journal, 24(2):264–289.
  * Stupfler et al., (2016) Stupfler, G. et al. (2016).  On the weak convergence of the kernel density estimator in the uniform topology.  Electronic Communications in Probability, 21.
  * van der Laan and Rose, (2011) van der Laan, M. and Rose, S. (2011).  Targeted Learning: Causal Inference for Observational and Experimental Data.  Springer New York.
  * van der Laan et al., (2018) van der Laan, M. J., Bibaut, A., and Luedtke, A. R. (2018).  CV-TMLE for Nonpathwise Differentiable Target Parameters.  In van der Laan, M. J. and Rose, S., editors, Targeted Learning in Data Science: Causal Inference for Complex Longitudinal Studies, pages 455–481. Springer International Publishing.
  * van der Laan et al., (2007) van der Laan, M. J., Polley, E. C., and Hubbard, A. E. (2007).  Super Learner.  Statistical Applications in Genetics and Molecular Biology, 6(1).
  * van der Vaart and van der Laan, (2006) van der Vaart, A. and van der Laan, M. J. (2006).  Estimating a survival distribution with current status data and high-dimensional covariates.  The International Journal of Biostatistics, 2(1).
  * van der Vaart and Wellner, (2011) van der Vaart, A. and Wellner, J. A. (2011).  A local maximal inequality under uniform entropy.  Electron. J. Statist., 5:192–203.
  * van der Vaart and Wellner, (1996) van der Vaart, A. W. and Wellner, J. A. (1996).  Weak Convergence and Empirical Processes.
  * Wasserman, (2006) Wasserman, L. (2006).  All of Nonparametric Statistics.  Springer Science & Business Media.
  * Weng et al., (2022) Weng, G., Doss, C. R., Wang, L., Moscovice, I., and Chantarat, T. (2022).  A nonparametric doubly robust test for a continuous treatment effect.  arXiv preprint arXiv:2202.03369.
  * Westling, (2022) Westling, T. (2022).  Nonparametric Tests of the Causal Null With Nondiscrete Exposures.  Journal of the American Statistical Association, 117(539):1551–1562.
  * Westling and Carone, (2020) Westling, T. and Carone, M. (2020).  A unified study of nonparametric inference for monotone functions.  Annals of Statistics, 48(2):1001–1024.
  * Westling et al., (2020) Westling, T., Gilbert, P., and Carone, M. (2020).  Causal isotonic regression.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 82(3):719–747.
  * (45) Wyatt, L. H., Peterson, G. C., Wade, T. J., Neas, L. M., and Rappold, A. G. (2020a).  The contribution of improved air quality to reduced cardiovascular mortality: Declines in socioeconomic differences over time.  Environment International, 136:105430.
  * (46) Wyatt, L. H., Peterson, G. C. L., Wade, T. J., Neas, L. M., and Rappold, A. G. (2020b).  Annual PM2.5 and cardiovascular mortality rate data: Trends modified by county socioeconomic status in 2,132 US counties.  Data in Brief, 30:105318.
  * Zhang and Xia, (2012) Zhang, W. and Xia, Y. (2012).  Twicing local linear kernel regression smoothers.  Journal of Nonparametric Statistics, 24(2):399–417.
  * Zhang et al., (2016) Zhang, Z., Zhou, J., Cao, W., and Zhang, J. (2016).  Causal inference with a quantitative exposure.  Statistical Methods in Medical Research, 25(1):315–335.

## Supplementary Material

##  Supplement A Guide to notation

We recall the data unit is $O:=(Y,A,W)$, and takes value in the sample space
$\mathcal{Y}\times\mathcal{A}\times\mathcal{W}$, where
$\mathcal{Y}\subseteq\mathbb{R}$, $\mathcal{A}\subseteq\mathbb{R}$, and
$\mathcal{W}\subseteq\mathbb{R}^{d}$. We observe $n$ IID replicates of this
data unit, which we call $O_{1},\dotsc,O_{n}$, and we define the unknown
distribution of $O_{i}$ as $P_{0}$. We index objects and functions by $P$ when
they depend on a distribution $P$ generating $O$. We index by $0$ when
$P=P_{0}$. We index by $n$ when the object is random depending on the sample
$O_{1},\dotsc,O_{n}$. Finally, we index by $\infty$ when the object is a limit
of an object indexed by $n$ that may or may not be equal to its true
counterpart depending on $P_{0}$.

  * •

$F_{P}(a):=P(A\leq a)$ and $f_{P}(a):=\frac{\partial}{\partial a}F_{P}(a)$ is
the marginal CDF and PDF of $A$ under $P$

  * •

$Q_{P}(\mathcal{X}):=P(W\in\mathcal{X})$ is the marginal distribution of $W$
under $P$

  * •

$\mu_{P}(a,w):=E_{P}(Y\mid A=a,W=w)$ as the outcome regression function under
$P$

  * •

$\theta_{P}(a):=E_{P}[\mu_{P}(a,W)]=\int\mu_{P}(a,w)\,dQ_{P}(w)$ is the
$G$-computed regression function

  * •

$g_{P}(a,w):=\left[\frac{\partial}{\partial a}P(A\leq a\mid
W=w)\right]/f_{P}(a)$ is the standardized conditional density of $A$ given
$W=w$

  * •

$Pf:=\int f(o)\,dP(o)$

  * •

$\|\beta\|_{P,q}:=(P|\beta|^{q})^{1/q}$

  * •

$\|\beta\|_{\mathcal{Z}}:=\sup_{z\in\mathcal{Z}}|\beta(z)|$

  * •

$\mathbb{P}_{n}$ is the empirical probability measure corresponding to
$O_{1},\dotsc,O_{n}$

  * •

$\mathbb{G}_{n}:=n^{1/2}(\mathbb{P}_{n}-P_{0})$

  * •

$Q_{n}$ is the marginal empirical distribution of $W_{1},\dotsc,W_{n}$

  * •

$F_{n}$ is the marginal empirical distribution of $A_{1},\dotsc,A_{n}$

  * •

$\mu_{n}$ is an estimator of $\mu_{0}$ with limit $\mu_{\infty}$

  * •

$g_{n}$ is an estimator of $g_{0}$ with limit $g_{\infty}$

  * •

$\psi_{P}(y,a,w):=\\{y-\mu_{P}(a,w)\\}/g_{P}(a,w)$

  * •

$\psi_{n}(y,a,w):=\\{y-\mu_{n}(a,w)\\}/g_{n}(a,w)$

  * •

$\psi_{\infty}(y,a,w):=\\{y-\mu_{\infty}(a,w)\\}/g_{\infty}(a,w)$

  * •

$\xi_{P}(y,a,w):=\psi_{P}(y,a,w)+\int\mu_{P}(a,\bar{w})\,dQ_{P}(\bar{w})$

  * •

$\xi_{n}(y,a,w):=\psi_{n}(y,a,w)+\int\mu_{n}(a,\bar{w})\,dQ_{n}(\bar{w})$

  * •

$\xi_{\infty}(y,a,w):=\psi_{\infty}(y,a,w)+\int\mu_{\infty}(a,\bar{w})\,dQ_{0}(\bar{w})$

  * •

$K_{h,a_{0}}(a):=K((a-a_{0})/h)/h$ for $K$ the kernel function

  * •

$c_{j}:=\int u^{j}K(u)\,du$.

  * •

$c_{j}^{*}:=\int u^{j}K^{2}(u)\,du$.

  * •

$c_{j,\tau}^{*}:=\int u^{j}K(u)K(\tau u)\,du$.

  * •

$\mathbf{S}_{k}:=(c_{i+j-2})_{1\leq i,j\leq k}$; i.e. the $k\times k$ matrix
with $(i,j)$ element $c_{i+j-2}$

  * •

$\mathbf{S}_{\tau}^{*}:=\int(1,u)^{T}(1,\tau u,[\tau u]^{2},[\tau
u]^{3})K(u)K(\tau u)\,du$

  * •

$e_{j}$ is the ($j+1$)-vector of zero’s with one at the $j$th position for any
integer $j\geq 1$.

  * •

$w_{h,a_{0},j}(a):=\left(1,[a-a_{0}]/h,\dotsc,[a-a_{0}]^{j}/h^{j}\right)^{T}$
for an integer $j\geq 1$;

  * •

$\mathbf{D}_{P,h,a_{0},j}:=P\left(w_{h,a_{0},j}K_{h,a_{0}}w^{T}_{h,a_{0},j}\right)$
for an integer $j\geq 1$;

  * •

$\mathbf{D}_{n,h,a_{0},j}:=\mathbb{P}_{n}\left(w_{h,a_{0},j}K_{h,a_{0}}w^{T}_{h,a_{0},j}\right)$
for an integer $j\geq 1$;

  * •

$\Gamma_{P,h,b,a_{0}}(a):=e_{1}^{T}\mathbf{D}_{P,h,a_{0},1}^{-1}w_{h,a_{0},1}(a)K_{h,a_{0}}(a)-e_{3}^{T}c_{2}(h/b)^{2}\mathbf{D}_{P,b,a_{0},3}^{-1}w_{b,a_{0},3}(a)K_{b,a_{0}}(a)$

  * •

$\Gamma_{n,h,b,a_{0}}(a):=e_{1}^{T}\mathbf{D}_{n,h,a_{0},1}^{-1}w_{h,a_{0},1}(a)K_{h,a_{0}}(a)-e_{3}^{T}c_{2}(h/b)^{2}\mathbf{D}_{n,b,a_{0},3}^{-1}w_{b,a_{0},3}(a)K_{b,a_{0}}(a)$

  * •

$\theta_{P,h,b}(a_{0}):=P(\Gamma_{P,h,b,a_{0}}\theta_{P})$

  * •

$\theta_{n,h,b}(a_{0}):=\mathbb{P}_{n}(\Gamma_{n,h,b,a_{0}}\xi_{n})$

  * •

$\begin{aligned}
\gamma_{P,h,b,a_{0}}(a)&:=e_{1}^{T}\mathbf{D}^{-1}_{P,h,a_{0},1}w_{h,a_{0},1}(a)K_{h,a_{0}}(a)w^{T}_{h,a_{0},1}(a)\mathbf{D}^{-1}_{P,h,a_{0},1}P\left(w_{h,a_{0},1}K_{h,a_{0}}\theta_{P}\right)\\\
&\qquad-
e_{3}^{T}c_{2}(h/b)^{2}\mathbf{D}_{P,b,a_{0},3}^{-1}w_{b,a_{0},3}(a)K_{b,a_{0}}(a)w^{T}_{b,a_{0},3}(a)\mathbf{D}^{-1}_{P,b,a_{0},3}P\left(w_{b,a_{0},3}K_{b,a_{0}}\theta_{P}\right)\end{aligned}$

  * •

$\begin{aligned}
\phi^{*}_{P,h,b,a_{0}}(y,a,w)&:=\Gamma_{P,h,b,a_{0}}(a)\xi_{P}(y,a,w)-\gamma_{P,h,b,a_{0}}(a)\\\
&\qquad+\int\Gamma_{P,h,b,a_{0}}(\bar{a})\left\\{\mu_{P}(\bar{a},w)-\int\mu_{P}(\bar{a},\bar{w})\,dQ_{P}(\bar{w})\right\\}\,dF_{P}(\bar{a})\end{aligned}$

  * •

$\begin{aligned}
\phi_{\infty,h,b,a_{0}}^{*}(y,a,w)&:=\Gamma_{0,h,b,a_{0}}(a)\xi_{\infty}(y,a,w)-\gamma_{0,h,b,a_{0}}(a)\\\
&\qquad+\int\Gamma_{0,h,b,a_{0}}(\bar{a})\left\\{\mu_{\infty}(\bar{a},w)-\int\mu_{\infty}(\bar{a},\bar{w})\,dQ_{0}(\bar{w})\right\\}\,dF_{0}(\bar{a})\end{aligned}$

  * •

$\begin{aligned}
\phi_{n,h,b,a_{0}}^{*}(y,a,w)&:=\Gamma_{n,h,b,a_{0}}(a)\xi_{n}(y,a,w)-\gamma_{n,h,b,a_{0}}(a)\\\
&\qquad+\int\Gamma_{n,h,b,a_{0}}(\bar{a})\left\\{\mu_{n}(\bar{a},w)-\int\mu_{n}(\bar{a},\bar{w})\,dQ_{n}(\bar{w})\right\\}\,dF_{n}(\bar{a})\end{aligned}$

  * •

$\tau_{n}:=h_{n}/b_{n}$

  * •

$\sigma_{0}^{2}(a):=E_{0}\left[\left\\{\xi_{\infty}(Y,A,W)-\theta_{0}(A)\right\\}^{2}\mid
A=a\right]$

  * •

$\sigma_{\infty,h,b}^{2}(a_{0}):=hP_{0}\left(\phi_{\infty,h,b,a_{0}}^{*}\right)^{2}$.

  * •

$\sigma_{n,h,b}^{2}(a_{0}):=h\mathbb{P}_{n}\left(\phi_{n,h,b,a_{0}}^{*}\right)^{2}$.

  * •

$\omega_{n}:=\sup_{a_{0}\in\mathcal{A}_{0}}\inf_{a\in\mathcal{A}_{n}}|a_{0}-a_{n}|$.

Throughout, we use $\lesssim$ to mean “less than up to a constant not
depending on $n$, $h$, or $b$.

##  Supplement B Conditions

Here, we re-state the conditions we use in our asymptotic results for
convenience. These are the same conditions as stated in the main text.

  1. (A1)

The kernel $K$ is a mean-zero, symmetric, nonnegative, and Lipschitz
continuous density function with support contained in $[-1,1]$. Additionally,
$K$ belongs to the linear span of the functions whose subgraph can be
represented as a finite number of Boolean operations among sets of the form
$\\{(s,u)\in\mathbb{R}\times\mathbb{R}:p(s,u)\leq\varphi(u)\\}$ where $p$ is a
polynomial and $\varphi$ is an arbitrary real function.

  2. (A2)

As $n\longrightarrow\infty$, the bandwidths $h=h_{n}$ and $b=b_{n}$ satisfy
$h_{n}\longrightarrow 0$, $nh_{n}\longrightarrow\infty$, $b_{n}\longrightarrow
0$, and $\tau_{n}:=h_{n}/b_{n}\longrightarrow\tau\in[0,\infty)$.

We define the following pseudo-distance for any $P_{0}$-square integrable
functions $\gamma_{1},\gamma_{2}:\mathcal{A}\times\mathcal{W}\to\mathbb{R}$,
$\mathcal{A}_{0}\subset\mathcal{A}$, and
$\mathcal{S}\subseteq\mathcal{A}\times\mathcal{W}$:

$\displaystyle
d(\gamma_{1},\gamma_{2};\mathcal{A}_{0},\mathcal{S}):=\sup_{a\in\mathcal{A}_{0}}\left\\{E_{0}[I_{\mathcal{S}}(a,W)\\{\gamma_{1}(a,W)-\gamma_{2}(a,W)\\}^{2}]\right\\}^{1/2}$

We also define $B_{\varepsilon}(a_{0})$ as the closed ball of radius
$\varepsilon$ centered at $a_{0}$.

  1. (A3)

There exist classes of functions $\mathcal{F}_{\mu}$ and $\mathcal{F}_{g}$
such that almost surely for all $n$ large enough,
$\mu_{0},\mu_{n}\in\mathcal{F}_{\mu}$, $g_{0},g_{n}\in\mathcal{F}_{g}$, and
for some constants $C_{j}\in(0,\infty)$, $V_{\mu}\in(0,1)$, and
$V_{g}\in(0,2)$:

    1. (a)

$|\mu|\leq C_{1}$ for all $\mu\in\mathcal{F}_{\mu}$, and $C_{2}\leq g\leq
C_{3}$ for all $g\in\mathcal{F}_{g}$;

    2. (b)

$\sup_{Q}\log N(\varepsilon,\mathcal{F}_{\mu},L_{2}(Q))\leq
C_{4}\varepsilon^{-V_{\mu}}$ and $\sup_{Q}\log
N(\varepsilon,\mathcal{F}_{g},L_{2}(Q))\leq C_{5}\varepsilon^{-V_{g}}$ for all
$\varepsilon>0$.

  2. (A4)

There exist $\mu_{\infty}\in\mathcal{F}_{\mu}$ and
$g_{\infty}\in\mathcal{F}_{g}$, $\delta_{1}>0$, and subsets
$\mathcal{S}_{1},\mathcal{S}_{2}$ and $\mathcal{S}_{3}$ of
$B_{\delta_{1}}(a_{0})\times\mathcal{W}$ such that
$\mathcal{S}_{1}\cup\mathcal{S}_{2}\cup\mathcal{S}_{3}=B_{\delta_{1}}(a_{0})\times\mathcal{W}$
and:

    1. (a)

$\mu_{\infty}(a,w)=\mu_{0}(a,w)$ for all
$(a,w)\in\mathcal{S}_{1}\cup\mathcal{S}_{3}$ and $g_{\infty}(a,w)=g_{0}(a,w)$
for all $(a,w)\in\mathcal{S}_{2}\cup\mathcal{S}_{3}$;

    2. (b)

$d(\mu_{n},\mu_{\infty};B_{\delta_{1}}(a_{0}),\mathcal{S}_{1})=o_{p}\left(\\{nh\\}^{-1/2}\right)$,
and $d(g_{n},g_{\infty};B_{\delta_{1}}(a_{0}),\mathcal{S}_{1})=o_{p}(1)$;

    3. (c)

$d(g_{n},g_{\infty};B_{\delta_{1}}(a_{0}),\mathcal{S}_{2})\\}=o_{p}\left(\\{nh\\}^{-1/2}\right)$,
and $d(\mu_{n},\mu_{\infty};B_{\delta_{1}}(a_{0}),\mathcal{S}_{1})=o_{p}(1)$;

    4. (d)

$d(\mu_{n},\mu_{\infty};B_{\delta_{1}}(a_{0}),\mathcal{S}_{3})d(g_{n},g_{\infty};B_{\delta_{1}}(a_{0}),\mathcal{S}_{3})=o_{p}\left(\\{nh\\}^{-1/2}\right)$.

  3. (A5)

It holds that:

    1. (a)

$\theta_{0}$ is twice continuously differentiable on $B_{\delta_{1}}(a_{0})$;

    2. (b)

$f_{0}$ is positive and Lipschitz continuous on $B_{\delta_{1}}(a_{0})$;

    3. (c)

there exist $\delta_{2}>0$ and $C_{6}<\infty$ such that
$E_{0}[|Y|^{2+\delta_{2}}\mid A=a,W=w]\leq C_{6}$ for all $a\in
B_{\delta_{1}}(a_{0})$ and $P_{0}$-almost every $w$ and
$E_{0}[|Y|^{4}]<\infty$; and

    4. (d)

$a\mapsto\sigma_{0}^{2}(a):=E_{0}\left[\\{\xi_{\infty}(Y,A,W)-\theta_{0}(A)\\}^{2}\mid
A=a\right]$ is bounded and continuous on $B_{\delta_{1}}(a_{0})$, where
$\xi_{\infty}:=\xi_{\mu_{\infty},g_{\infty},Q_{0}}$ is the limiting pseudo-
outcome.

Let $\mathcal{A}_{0}$ be a compact subset of $\mathcal{A}$ over which we wish
to perform uniform inference. For $\delta>0$ define $\mathcal{A}_{\delta}$ as
the $\delta$-enlargement of $\mathcal{A}_{0}$; that is, the set of
$a\in\mathbb{R}$ such that there exists $a_{0}\in\mathcal{A}_{0}$ with
$|a-a_{0}|\leq\delta$. We then assume there exists $\delta_{3}>0$ such that
the following hold.

  1. (A6)

The constant $V:=\max\\{V_{\mu},V_{g}\\}$, for $V_{\mu}$ and $V_{g}$ defined
in (A3) and the bandwidth $h$ satisfies $n\left[h/(\log
n)\right]^{\frac{2+V}{2-V}}\longrightarrow\infty$ and
$nh^{3}\longrightarrow\infty$.

  2. (A7)

There exist $\mu_{\infty}\in\mathcal{F}_{\mu}$,
$g_{\infty}\in\mathcal{F}_{g}$, $\delta_{3}>0$ and subsets
$\mathcal{S}_{1}^{\prime},\mathcal{S}_{2}^{\prime}$, and
$\mathcal{S}_{3}^{\prime}$ of $\mathcal{A}_{\delta_{3}}\times\mathcal{W}$ such
that
$\mathcal{S}_{1}^{\prime}\cup\mathcal{S}_{2}^{\prime}\cup\mathcal{S}_{3}^{\prime}=\mathcal{A}_{\delta_{3}}\times\mathcal{W}$
and:

    1. (a)

$\mu_{\infty}(a,w)=\mu_{0}(a,w)$ for all
$(a,w)\in\mathcal{S}_{1}^{\prime}\cup\mathcal{S}_{3}^{\prime}$ and
$g_{\infty}(a,w)=g_{0}(a,w)$ for all
$(a,w)\in\mathcal{S}_{2}^{\prime}\cup\mathcal{S}_{3}^{\prime}$;

    2. (b)

$d(\mu_{n},\mu_{\infty};\mathcal{A}_{\delta_{3}},\mathcal{S}_{1}^{\prime})=o_{p}\left(\\{nh\log
n\\}^{-1/2}\right)$;

    3. (c)

$d(g_{n},g_{\infty};\mathcal{A}_{\delta_{3}},\mathcal{S}_{2}^{\prime})\\}=o_{p}\left(\\{nh\log
n\\}^{-1/2}\right)$;

    4. (d)

$d(\mu_{n},\mu_{\infty};\mathcal{A}_{\delta_{3}},\mathcal{S}_{3})d(g_{n},g_{\infty};\mathcal{A}_{\delta_{3}},\mathcal{S}_{3})=o_{p}\left(\\{nh\log
n\\}^{-1/2}\right)$.

    5. (e)

$d(\mu_{n},\mu_{\infty};\mathcal{A}_{\delta_{3}},\mathcal{A}\times\mathcal{W})$
and
$d(g_{n},g_{\infty};\mathcal{A}_{\delta_{3}},\mathcal{A}\times\mathcal{W})$
are both $o_{p}\left(h^{\frac{V}{2(2-V)}}\\{\log n\\}^{-\frac{1}{2-V}}\right)$

  3. (A8)

It holds that:

    1. (a)

$\theta_{0}$ is twice continuously differentiable with Hölder-continuous
second derivative on $\mathcal{A}_{\delta_{3}}$;

    2. (b)

$f_{0}$ is Lipschitz continuous and bounded away from 0 and $\infty$ on
$\mathcal{A}_{\delta_{3}}$;

    3. (c)

$|Y|$ is $P_{0}$-almost surely bounded; and

    4. (d)

$a\mapsto E_{0}\left[\\{\xi_{\infty}(Y,A,W)-\theta_{0}(A)\\}^{2}\mid
A=a\right]$ is continuous on $\mathcal{A}_{\delta_{3}}$.

##  Supplement C Efficient influence function of the local parameter

###### Proof of Lemma 1.

We let $\\{P_{\varepsilon}:|\varepsilon|\leq\delta\\}$ be a one-dimensional
quadratic mean differentiable path contained in the model $\mathcal{M}$ such
that $P_{\epsilon=0}=P_{0}$. We let $(y,a,w)\mapsto\dot{\ell}_{0}(y,a,w)$ be
the score function of the path at $\varepsilon=0$, which necessarily satisfies
$P_{0}\dot{\ell}_{0}=0$ and $P_{0}\dot{\ell}^{2}_{0}<\infty$. We define the
marginal score functions $\dot{\ell}_{0}(a,w):=E_{0}[\dot{\ell}_{0}\mid
A=a,W=w]$, $\dot{\ell}_{0}(a):=E_{0}[\dot{\ell}_{0}\mid A=a]$, and
$\dot{\ell}_{0}(w):=E[\dot{\ell}_{0}\mid W=w]$, which all have mean zero under
$P_{0}$ by the tower property. We also define the conditional score function
$\dot{\ell}_{0}(y\mid a,w):=\dot{\ell}_{0}(y,a,w)-\dot{\ell}(a,w)$, which
satisfies $E_{0}[\dot{\ell}_{0}(Y\mid A,W)\mid A=a,W=w]=0$ for $P_{0}$-almost
every $(a,w)$. For convenience, we denote objects depending on
$P_{\varepsilon}$ with a subscript $\varepsilon$.

To demonstrate that $\phi^{*}_{0,h,b,a_{0}}$ is the nonparametric efficient
influence function of $\theta_{0,h,b}(a_{0})$, we show that

$\left.\frac{\partial}{\partial\varepsilon}\theta_{\varepsilon,h,b}(a_{0})\right|_{\varepsilon=0}=P_{0}(\phi^{*}_{0,h,b,a_{0}}\dot{\ell}_{0})$

for each fixed $h,b>0$ and $a_{0}\in\mathcal{A}$. We define

$\displaystyle\theta_{P,h}(a_{0})$
$\displaystyle:=e_{1}^{T}\mathbf{D}^{-1}_{P,h,a_{0}}P_{0}\left(w_{h,a_{0},1}K_{h,a_{0}}\theta_{P}\right)\text{
and}$ $\displaystyle\theta_{P,b}^{(2)}(a_{0})$
$\displaystyle:=2b^{-2}e_{3}^{T}\mathbf{D}^{-1}_{P,b,a_{0}}P\left(w_{b,a_{0},3}K_{b,a_{0}}\theta_{P}\right),$

so that $\theta_{P,h,b}(a_{0})$ can be expressed as

$\displaystyle\theta_{P,h,b}(a_{0})=\theta_{P,h}(a_{0})-\tfrac{1}{2}c_{2}h^{2}\theta^{(2)}_{P,b}(a_{0}).$

Therefore,

$\left.\frac{\partial}{\partial\varepsilon}\theta_{\varepsilon,h,b}(a_{0})\right|_{\varepsilon=0}=\left.\frac{\partial}{\partial\varepsilon}\theta_{\varepsilon,h}(a_{0})\right|_{P_{\varepsilon}=0}-\tfrac{1}{2}c_{2}h^{2}\left.\frac{\partial}{\partial\varepsilon}\theta^{(2)}_{\varepsilon,b}(a_{0})\right|_{\varepsilon=0}.$

We now provide the derivation of
$\left.\frac{\partial}{\partial\varepsilon}\theta_{\varepsilon,h}(a_{0})\right|_{\varepsilon=0}$.
We first have by the product rule that

$\displaystyle\left.\frac{\partial}{\partial\varepsilon}\theta_{\varepsilon,h}(a_{0})\right|_{\varepsilon=0}$
$\displaystyle=\left.\frac{\partial}{\partial\varepsilon}e_{1}^{T}\mathbf{D}^{-1}_{\varepsilon,h,a_{0}}P_{\varepsilon}\left(w_{h,a_{0},1}K_{h,a_{0}}\theta_{\varepsilon}\right)\right|_{\varepsilon=0}$
$\displaystyle=\left.\frac{\partial}{\partial\varepsilon}e_{1}^{T}\mathbf{D}^{-1}_{\varepsilon,h,a_{0}}\right|_{\varepsilon=0}P_{0}\left(w_{h,a_{0},1}K_{h,a_{0}}\theta_{0}\right)+e_{1}^{T}\mathbf{D}^{-1}_{0,h,a_{0},1}\left.\frac{\partial}{\partial\varepsilon}P_{\varepsilon}\left(w_{h,a_{0},1}K_{h,a_{0}}\theta_{0}\right)\right|_{\varepsilon=0}$
$\displaystyle\qquad+e_{1}^{T}\mathbf{D}^{-1}_{0,h,a_{0},1}\left.\frac{\partial}{\partial\varepsilon}P_{0}\left(w_{h,a_{0},1}K_{h,a_{0}}\theta_{\varepsilon}\right)\right|_{\varepsilon=0}.$

For the first term, by the definition of $\mathbf{D}_{P,h,a_{0}}$ as the mean
of a fixed and uniformly bounded function and the chain rule, we have

$\displaystyle\left.\frac{\partial}{\partial\varepsilon}e_{1}^{T}\mathbf{D}^{-1}_{\varepsilon,h,a_{0}}\right|_{\varepsilon=0}$
$\displaystyle=-e_{1}^{T}\mathbf{D}^{-1}_{0,h,a_{0}}\left.\frac{\partial}{\partial\varepsilon}\mathbf{D}_{\varepsilon,h,a_{0}}\right|_{\varepsilon=0}\mathbf{D}^{-1}_{0,h,a_{0}}$
$\displaystyle=-e_{1}^{T}\mathbf{D}^{-1}_{0,h,a_{0}}P_{0}\left(w_{h,a_{0},1}K_{h,a_{0}}w^{T}_{h,a_{0},1}\dot{\ell}_{0}\right)\mathbf{D}^{-1}_{0,h,a_{0}}.$

Similarly, for the second term, we have

$\displaystyle\left.\frac{\partial}{\partial\varepsilon}P_{\varepsilon}\left(w_{h,a_{0},1}K_{h,a_{0}}\theta_{0}\right)\right|_{\varepsilon=0}=P_{0}\left(w_{h,a_{0},1}K_{h,a_{0}}\theta_{0}\dot{\ell}_{0}\right).$

For the third term, we have by definition of $\theta_{P}$

$\displaystyle\left.\frac{\partial}{\partial\varepsilon}P_{0}\left(w_{h,a_{0},1}K_{h,a_{0}}\theta_{\varepsilon}\right)\right|_{\varepsilon=0}$
$\displaystyle=\left.\frac{\partial}{\partial\varepsilon}\int
w_{h,a_{0},1}(a)K_{h,a_{0}}(a)\iint y\,dP_{\varepsilon}(y\mid
a,w)\,dQ_{\varepsilon}(w)\,dF_{0}(a)\right|_{\varepsilon=0}$
$\displaystyle=\int w_{h,a_{0},1}(a)K_{h,a_{0}}(a)\iint y\dot{\ell}_{0}(y\mid
a,w)\,dP_{0}(y\mid a,w)\,dQ_{0}(w)\,dF_{0}(a)$ $\displaystyle\qquad+\int
w_{h,a_{0},1}(a)K_{h,a_{0}}(a)\iint y\,dP_{0}(y\mid
a,w)\dot{\ell}_{0}(w)\,dQ_{0}(w)\,dF_{0}(a).$

By properties of score functions and since $g_{0}(a,w)>0$ for $|a-a_{0}|\leq
h$ and $Q_{0}$-a.e. $w$, the first term equals

$\displaystyle\iiint
w_{h,a_{0},1}(a)K_{h,a_{0}}(a)\frac{y}{g_{0}(a,w)}\dot{\ell}_{0}(y\mid
a,w)\,dP_{0}(y,a,w)$ $\displaystyle\qquad=\iiint
w_{h,a_{0},1}(a)K_{h,a_{0}}(a)\frac{y-\mu_{0}(a,w)}{g_{0}(a,w)}\dot{\ell}_{0}(y,a,w)\,dP_{0}(y,a,w).$

The second term equals

$\displaystyle\iint
w_{h,a_{0},1}(a)K_{h,a_{0}}(a)\mu_{0}(a,w)\,dF_{0}(a)\dot{\ell}_{0}(w)\,dQ_{0}(w)$
$\displaystyle=\int\eta_{0}(w)\dot{\ell}_{0}(w)\,dQ_{0}(w)$
$\displaystyle=\iiint\eta_{0}(w)\dot{\ell}_{0}(y,a,w)\,dP_{0}(y,a,w),$

where $\eta_{0}(w):=\int
w_{h,a_{0},1}(a)K_{h,a_{0}}(a)\mu_{0}(a,w)\,dF_{0}(a)$. Putting it together,
we have that
$\left.\frac{\partial}{\partial\varepsilon}\theta_{\varepsilon,h}(a_{0})\right|_{\varepsilon=0}$
equals

$\displaystyle-
e_{1}^{T}\mathbf{D}^{-1}_{0,h,a_{0},1}P_{0}\left(w_{h,a_{0},1}K_{h,a_{0}}w^{T}_{h,a_{0},1}\dot{\ell}_{0}\right)\mathbf{D}^{-1}_{0,h,a_{0},1}P_{0}\left(w_{h,a_{0},1}K_{h,a_{0}}\theta_{0}\right)$
$\displaystyle\qquad+e_{1}^{T}\mathbf{D}^{-1}_{0,h,a_{0},1}P_{0}\left(w_{h,a_{0},1}K_{h,a_{0}}\theta_{0}\dot{\ell}_{0}\right)$
$\displaystyle\qquad+e_{1}^{T}\mathbf{D}^{-1}_{0,h,a_{0},1}\iiint\left\\{\left[w_{h,a_{0},1}(a)K_{h,a_{0}}(a)\left\\{\frac{y-\mu_{0}(a,w)}{g_{0}(a,w)}\right\\}+\eta_{0}(w)\right]\dot{\ell}_{0}(y,a,w)\right\\}\,dP_{0}(y,a,w)$
$\displaystyle=e_{1}^{T}\mathbf{D}^{-1}_{0,h,a_{0},1}\iiint\left\\{\left[w_{h,a_{0},1}(a)K_{h,a_{0}}(a)\left\\{\frac{y-\mu_{0}(a,w)}{g_{0}(a,w)}+\theta_{0}(a)-\gamma_{0,h,a_{0}}(a)\right\\}+\eta_{0}(w)\right]\dot{\ell}_{0}(y,a,w)\right\\}\,dP_{0}(y,a,w)$

where
$\gamma_{0,h,a_{0}}(a):=w^{T}_{h,a_{0},1}(a)\mathbf{D}^{-1}_{0,h,a_{0},1}P_{0}\left(w_{h,a_{0},1}K_{h,a_{0}}\theta_{0}\right)$.
Therefore, the uncentered influence function of $\theta_{P,h}(a_{0})$ at
$P=P_{0}$ is

$\phi_{0,h,a_{0}}(y,a,w):=e_{1}^{T}\mathbf{D}^{-1}_{0,h,a_{0},1}\left[w_{h,a_{0},1}K_{h,a_{0}}\left\\{\frac{y-\mu_{0}(a,w)}{g_{0}(a,w)}+\theta_{0}(a)-\gamma_{0,h,a_{0}}(a)\right\\}+\eta_{0}(w)\right]$

Noting that

$P_{0}\left(w_{h,a_{0},1}K_{h,a_{0}}\gamma_{0,h,a_{0}}\right)=P_{0}\eta_{0}=P_{0}\left(w_{h,a_{0},1}K_{h,a_{0}}\theta_{0}\right),$

the mean of $\phi_{0,h,a_{0}}$ is

$\displaystyle P_{0}\phi_{0,h,a_{0}}$
$\displaystyle=e_{1}^{T}\mathbf{D}^{-1}_{0,h,a_{0},1}P_{0}\left(w_{h,a_{0},1}K_{h,a_{0}}\theta_{0}\right).$

Thus, the centered influence function of $\theta_{P,h}$ at $P_{0}$ is

$\displaystyle\phi^{*}_{0,h,a_{0}}(y,a,w)$
$\displaystyle=\phi_{0,h,a_{0}}(y,a,w)-P_{0}\phi_{0,h,a_{0}}$
$\displaystyle=e_{1}^{T}\mathbf{D}^{-1}_{0,h,a_{0},1}w_{h,a_{0},1}(a)K_{h,a_{0}}(a)\left\\{\frac{y-\mu_{0}(a,w)}{g_{0}(a,w)}+\theta_{0}(a)-\gamma_{0,h,a_{0}}(a)\right\\}$
$\displaystyle\qquad+e_{1}^{T}\mathbf{D}^{-1}_{0,h,a_{0},1}\int
w_{h,a_{0},1}(a)K_{h,a_{0}}(a)\\{\mu_{0}(a,w)-\theta_{0}(a)\\}\,dF_{0}(a)$

By a similar argument, the influence function of $\theta_{P,b}^{(2)}(a_{0})$
at $P_{0}$ is

$\displaystyle\phi^{*}_{0,b,a_{0}}(y,a,w)$
$\displaystyle=2b^{-2}e_{3}^{T}\mathbf{D}^{-1}_{0,b,a_{0},3}w_{b,a_{0},3}(a)K_{b,a_{0}}(a)\left\\{\frac{y-\mu_{0}(a,w)}{g_{0}(a,w)}+\theta_{0}(a)-\gamma_{0,b,a_{0}}(a)\right\\}$
$\displaystyle\qquad+2b^{-2}e_{3}^{T}\mathbf{D}^{-1}_{0,b,a_{0},3}\int
w_{b,a_{0},3}(a)K_{b,a_{0}}(a)\\{\mu_{0}(a,w)-\theta_{0}(a)\\}\,dF_{0}(a)$

where
$\gamma_{0,b,a_{0}}(a):=w^{T}_{b,a_{0},3}(a)\mathbf{D}^{-1}_{0,b,a_{0},3}P_{0}\left(w_{b,a_{0},3}K_{b,a_{0}}\theta_{0}\right)$.
Hence, the influence function of $\theta_{P,h,b}$ at $P_{0}$ is

$\displaystyle\phi^{*}_{0,h,b,a_{0}}(y,a,w)=\phi^{*}_{0,h,a_{0}}(y,a,w)-\tfrac{1}{2}c_{2}h^{2}\phi^{*}_{0,b,a_{0}}(y,a,w)$
$\displaystyle\qquad=\Gamma_{0,h,b,a_{0}}(a)\left\\{\frac{y-\mu_{0}(a,w)}{g_{0}(a,w)}+\theta_{0}(a)\right\\}-\gamma_{0,h,b,a_{0}}(a)+\int\Gamma_{0,h,b,a_{0}}(a)\\{\mu_{0}(a,w)-\theta_{0}(a)\\}\,dF_{0}(a).$

∎

##  Supplement D Decomposition of the estimator

In this section, we present a decomposition of the estimator into six
remainder terms. Before we present the decomposition, we provide a supporting
Lemma showing that the pseudo-outcome function possesses a double robustness
property under (A4)(a).

######  Lemma 2.

If (A4)(a) holds, then

$\displaystyle\theta_{0}(a)=E_{0}\left\\{\frac{Y-\mu_{\infty}(A,W)}{g_{\infty}(A,W)}+\int\mu_{\infty}(A,w)\,dQ_{0}(w)\,\Big{|}\,A=a\right\\}$

for $F_{0}$-almost every $a\in B_{\delta_{1}}(a_{0})$.

###### Proof of Lemma 2.

By the tower property and since $dP_{0}(w\mid A=a)=g_{0}(a,w)\,dQ_{0}(w)$, we
have for any $a$ in the support of $A$ that

$\displaystyle E_{0}\left\\{\xi_{\infty}(Y,A,W)\mid A=a\right\\}$
$\displaystyle=E_{0}\left\\{\frac{Y-\mu_{\infty}(A,W)}{g_{\infty}(A,W)}+\int\mu_{\infty}(A,w)\,dQ_{0}(w)\,\Big{|}\,A=a\right\\}$
$\displaystyle=\int\frac{\mu_{0}(a,w)-\mu_{\infty}(a,w)}{g_{\infty}(a_{0},w)}g_{0}(a,w)\,dQ_{0}(w)+\int\mu_{\infty}(a,w)\,dQ_{0}(w)$
$\displaystyle=\int\\{\mu_{0}(a,w)-\mu_{\infty}(a,w)\\}\left\\{\frac{g_{0}(a,w)}{g_{\infty}(a,w)}-1\right\\}\,dQ_{0}(w)+\int\mu_{0}(a,w)\,dQ_{0}(w)$
$\displaystyle=\int\\{\mu_{0}(a,w)-\mu_{\infty}(a,w)\\}\left\\{\frac{g_{0}(a,w)}{g_{\infty}(a,w)}-1\right\\}\,dQ_{0}(w)+\theta_{0}(a).$

Hence,

$\displaystyle\int_{B_{\delta_{1}}(a_{0})}\left|E_{0}\left\\{\xi_{\infty}(Y,A,W)\mid
A=a\right\\}-\theta_{0}(a)\right|\,dF_{0}(a)$
$\displaystyle\qquad=\int_{B_{\delta_{1}}(a_{0})}\left|\int\\{\mu_{0}(a,w)-\mu_{\infty}(a,w)\\}\left\\{\frac{g_{0}(a,w)}{g_{\infty}(a,w)}-1\right\\}\,dQ_{0}(w)\right|\,dF_{0}(a)$
$\displaystyle\qquad\leq\iint_{B_{\delta_{1}}(a_{0})\times\mathcal{W}}\left|\\{\mu_{0}(a,w)-\mu_{\infty}(a,w)\\}\left\\{\frac{g_{0}(a,w)}{g_{\infty}(a,w)}-1\right\\}\right|\,d(F_{0}\times
Q_{0})(a,w)$
$\displaystyle\qquad=\sum_{j=1}^{3}\iint_{\mathcal{S}_{j}}\left|\\{\mu_{0}(a,w)-\mu_{\infty}(a,w)\\}\left\\{\frac{g_{0}(a,w)}{g_{\infty}(a,w)}-1\right\\}\right|\,d(F_{0}\times
Q_{0})(a,w).$

By (A4)(a), at least one of the two terms inside the integral is zero for each
of the three integrals over $\mathcal{S}_{1}$, $\mathcal{S}_{2}$, and
$\mathcal{S}_{3}$. Hence, the above expression is zero, and the claim follows.
∎

We now provide the decomposition of the estimator that we will use throughout
our results. The decomposition involves a leading empirical mean,
$\mathbb{P}_{n}\phi^{*}_{\infty,h,b,a_{0}}$, which drives the first-order
asymptotic behavior of the estimator, and six remainder terms. We separate the
remainder into six terms because the methods of controlling these remainders
are conceptually distinct, as we will see in the ensuing results.

######  Lemma 3.

If (A1) and (A4)(a) hold, then the following expansion of the estimator holds
for all $h\in(0,\delta_{1})$:

$\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})=\mathbb{P}_{n}\phi^{*}_{\infty,h,b,a_{0}}+\sum_{j=1}^{6}\mathrel{R}_{n,h,b,a_{0},j},$

where

$\displaystyle\mathrel{R}_{n,h,b,a_{0},1}$
$\displaystyle:=\theta_{0,h,b}(a_{0})-\theta_{0}(a_{0}),$
$\displaystyle\mathrel{R}_{n,h,b,a_{0},2}$
$\displaystyle:=(\mathbb{P}_{n}-P_{0})\left[\left\\{\Gamma_{0,h,b,a_{0}}\left(\psi_{n}+\int\mu_{n}\,dQ_{0}\right)+\int\Gamma_{0,h,b,a_{0}}\mu_{n}\,dF_{0}\right\\}\right.$
$\displaystyle\qquad\qquad\qquad\qquad\left.-\left\\{\Gamma_{0,h,b,a_{0}}\left(\psi_{\infty}+\int\mu_{\infty}\,dQ_{0}\right)+\int\Gamma_{0,h,b,a_{0}}\mu_{\infty}\,dF_{0}\right\\}\right],$
$\displaystyle\mathrel{R}_{n,h,b,a_{0},3}$
$\displaystyle:=(\mathbb{P}_{n}-P_{0})\left\\{\left(\Gamma_{n,h,b,a_{0}}-\Gamma_{0,h,b,a_{0}}\right)\left(\psi_{n}+\int\mu_{n}\,dQ_{0}\right)+\int\left(\Gamma_{n,h,b,a_{0}}-\Gamma_{0,h,b,a_{0}}\right)\mu_{n}\,dF_{0}\right\\},$
$\displaystyle\mathrel{R}_{n,h,b,a_{0},4}$
$\displaystyle:=\iint\Gamma_{n,h,b,a_{0}}\\{\mu_{n}-\mu_{0}\\}\left\\{1-\frac{g_{0}}{g_{n}}\right\\}\,dF_{0}\,dQ_{0},$
$\displaystyle\mathrel{R}_{n,h,b,a_{0},5}$
$\displaystyle:=\iint\Gamma_{n,h,b,a_{0}}(a)\mu_{n}(a,w)\,d(Q_{n}-Q_{0})(w)\,d(F_{0}-F_{n})(a),\text{
and}$ $\displaystyle\mathrel{R}_{n,h,b,a_{0},6}$
$\displaystyle:=e_{1}^{T}\mathbf{D}^{-1}_{0,h,a_{0},1}\left(\mathbf{D}_{0,h,a_{0},1}-\mathbf{D}_{n,h,a_{0},1}\right)\left(\mathbf{D}^{-1}_{n,h,a_{0},1}-\mathbf{D}^{-1}_{0,h,a_{0},1}\right)P_{0}\left(w_{h,a_{0},1}K_{1,a_{0}}\theta_{0}\right)$
$\displaystyle\qquad-
e_{3}^{T}c_{2}\tau_{n}^{2}\mathbf{D}^{-1}_{0,b,a_{0},3}\left(\mathbf{D}_{0,b,a_{0},3}-\mathbf{D}_{n,b,a_{0},3}\right)\left(\mathbf{D}^{-1}_{n,b,a_{0},3}-\mathbf{D}^{-1}_{0,b,a_{0},3}\right)P_{0}\left(w_{b,a_{0},3}K_{b,a_{0}}\theta_{0}\right).$

###### Proof of Lemma 3.

It is straightforward to see that
$P_{0}\gamma_{0,h,b,a_{0}}=P_{0}\left(\Gamma_{0,h,b,a_{0}}\theta_{0}\right)$.
Since (A4)(a) holds by assumption, Lemma 2 implies that
$E_{0}\left\\{\xi_{\infty}(Y,A,W)\mid A=a\right\\}=\theta_{0}(a)$ for
$F_{0}$-almost every $a\in B_{\delta_{1}}(a_{0})$, and since (A1) holds, for
all $h\leq\delta_{1}$, $\Gamma_{0,h,b,a_{0}}(a)=0$ for $a\notin
B_{\delta_{1}}(a_{0})$. Thus,

$P_{0}\left(\Gamma_{0,h,b,a_{0}}\xi_{\infty}\right)=E_{0}\left[\Gamma_{0,h,b,a_{0}}(A)E_{0}\left\\{\xi_{\infty}(Y,A,W)\mid
A\right\\}\right]=E_{0}\left[\Gamma_{0,h,b,a_{0}}(A)\theta_{0}(A)\right]=P_{0}\left(\Gamma_{0,h,b,a_{0}}\theta_{0}\right).$

Therefore,

$\displaystyle P_{0}\phi_{\infty,h,b,a_{0}}^{*}$
$\displaystyle=P_{0}\left(\Gamma_{0,h,b,a_{0}}\xi_{\infty}\right)-P_{0}\gamma_{0,h,b,a_{0}}+\iint\Gamma_{0,h,b,a_{0}}\mu_{\infty}\,dF_{0}\,dQ_{0}-\iint\Gamma_{0,h,b,a_{0}}\mu_{\infty}\,dF_{0}\,dQ_{0}$
$\displaystyle=P_{0}\left(\Gamma_{0,h,b,a_{0}}\theta_{0}\right)-P_{0}\left(\Gamma_{0,h,b,a_{0}}\theta_{0}\right)=0.$

We now define

$\displaystyle\gamma_{n,h,b,a_{0}}(a)$
$\displaystyle:=e_{1}^{T}\mathbf{D}^{-1}_{n,h,a_{0},1}w_{h,a_{0},1}(a)K_{h,a_{0}}(a)w_{h,a_{0},1}^{T}(a)\mathbf{D}^{-1}_{n,h,a_{0},1}\mathbb{P}_{n}\left(w_{h,a_{0},1}K_{h,a_{0}}\xi_{n}\right)$
$\displaystyle\qquad-
e_{3}^{T}c_{2}\tau_{n}^{2}\mathbf{D}_{n,b,a_{0},3}^{-1}w_{b,a_{0},3}(a)K_{b,a_{0}}(a)w_{b,a_{0},3}^{T}(a)\mathbf{D}_{n,b,a_{0},3}^{-1}\mathbb{P}_{n}\left(w_{b,a_{0},3}K_{b,a_{0}}\xi_{n}\right),$

and we note that

$\displaystyle\mathbb{P}_{n}\gamma_{n,h,b,a_{0}}$
$\displaystyle=e_{1}^{T}\mathbf{D}^{-1}_{n,h,a_{0},1}\mathbb{P}_{n}\left(w_{h,a_{0},1}K_{h,a_{0}}w_{h,a_{0},1}^{T}\right)\mathbf{D}^{-1}_{n,h,a_{0},1}\mathbb{P}_{n}\left(w_{h,a_{0},1}K_{h,a_{0}}\xi_{n}\right)$
$\displaystyle\qquad-
e_{3}^{T}c_{2}\tau_{n}^{2}\mathbf{D}_{n,b,a_{0},3}^{-1}\mathbb{P}_{n}\left(w_{b,a_{0},3}K_{b,a_{0}}w_{b,a_{0},3}^{T}\right)\mathbf{D}_{n,b,a_{0},3}^{-1}\mathbb{P}_{n}\left(w_{b,a_{0},3}K_{b,a_{0}}\xi_{n}\right)$
$\displaystyle=e_{1}^{T}\mathbf{D}^{-1}_{n,h,a_{0},1}\mathbf{D}_{n,h,a_{0},1}\mathbf{D}^{-1}_{n,h,a_{0},1}\mathbb{P}_{n}\left(w_{h,a_{0},1}K_{h,a_{0}}\xi_{n}\right)$
$\displaystyle\qquad-
e_{3}^{T}c_{2}\tau_{n}^{2}\mathbf{D}_{n,b,a_{0},3}^{-1}\mathbf{D}_{n,b,a_{0},3}\mathbf{D}_{n,b,a_{0},3}^{-1}\mathbb{P}_{n}\left(w_{b,a_{0},3}K_{b,a_{0}}\xi_{n}\right)$
$\displaystyle=e_{1}^{T}\mathbf{D}^{-1}_{n,h,a_{0},1}\mathbb{P}_{n}\left(w_{h,a_{0},1}K_{h,a_{0}}\xi_{n}\right)-e_{3}^{T}c_{2}\tau_{n}^{2}\mathbf{D}_{n,b,a_{0},3}^{-1}\mathbb{P}_{n}\left(w_{b,a_{0},3}K_{b,a_{0}}\xi_{n}\right)$
$\displaystyle=\mathbb{P}_{n}\left(e_{1}^{T}\mathbf{D}^{-1}_{n,h,a_{0},1}w_{h,a_{0},1}K_{h,a_{0}}\xi_{n}-e_{3}^{T}c_{2}\tau_{n}^{2}\mathbf{D}_{n,b,a_{0},3}^{-1}w_{b,a_{0},3}K_{b,a_{0}}\xi_{n}\right)$
$\displaystyle=\mathbb{P}_{n}\left(\Gamma_{n,h,b,a_{0}}\xi_{n}\right).$

We then define

$\displaystyle\phi_{n,h,b,a_{0}}^{*}(y,a,w)$
$\displaystyle:=\Gamma_{n,h,b,a_{0}}(a)\xi_{n}(y,a,w)-\gamma_{n,h,b,a_{0}}(a)+\int\Gamma_{n,h,b,a_{0}}(\bar{a})\left\\{\mu_{n}(\bar{a},w)-\int\mu_{n}(\bar{a},\bar{w})\,dQ_{n}(\bar{w})\right\\}\,dF_{n}(\bar{a}),$

and we note that

$\displaystyle\mathbb{P}_{n}\phi_{n,h,b,a_{0}}^{*}$
$\displaystyle=\mathbb{P}_{n}\left(\Gamma_{n,h,b,a_{0}}\xi_{n}\right)-\mathbb{P}_{n}\gamma_{n,h,b,a_{0}}+\iint\Gamma_{n,h,b,a_{0}}\mu_{n}\,dF_{n}\,dQ_{n}-\iint\Gamma_{n,h,b,a_{0}}\mu_{n}\,dQ_{n}\,dF_{n}$
$\displaystyle=\mathbb{P}_{n}\gamma_{n,h,b,a_{0}}-\mathbb{P}_{n}\gamma_{n,h,b,a_{0}}=0.$

By adding and subtracting terms and using the fact that
$P_{0}\phi_{\infty,h,b,a_{0}}^{*}=\mathbb{P}_{n}\phi_{n,h,b,a_{0}}^{*}=0$, we
can now write

$\displaystyle\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})$
$\displaystyle=\mathbb{P}_{n}\phi_{\infty,h,b,a_{0}}^{*}+(\mathbb{P}_{n}-P_{0})\left(\phi_{n,h,b,a_{0}}^{*}-\phi_{\infty,h,b,a_{0}}^{*}\right)+\left\\{\theta_{n,h,b}(a_{0})-\theta_{0,h,b}(a_{0})+P_{0}\phi_{n,h,b,a_{0}}^{*}\right\\}$
$\displaystyle\qquad+\left\\{\theta_{0,h,b}(a_{0})-\theta_{0}(a_{0})\right\\}.$

We note that the last summand in braces equals $R_{n,h,b,a_{0},1}$. Up to this
last term, this is a standard first-order expansion of the estimator
$\theta_{n,h,b}(a_{0})$ of the pathwise differentiable parameter
$\theta_{0,h,b}(a_{0})$ for fixed $h$ and $b$, though our analysis in
subsequent results will consider the case where $h$ and $b$ go to zero as $n$
grows.

We now enter a calculation showing that the second and third summands above
equal $\sum_{j=2}^{6}R_{n,h,b,a_{0},j}$. For convenience, in this derivation
we omit the subscripts $h$, $b$, and $a_{0}$ when it is clear. We have

$\displaystyle(\mathbb{P}_{n}-P_{0})\left(\phi_{n,h,b,a_{0}}^{*}-\phi_{\infty,h,b,a_{0}}^{*}\right)+\left\\{\theta_{n,h,b}(a_{0})-\theta_{0,h,b}(a_{0})+P_{0}\phi_{n,h,b,a_{0}}^{*}\right\\}$
$\displaystyle\qquad=(\mathbb{P}_{n}-P_{0})\left(\left\\{\Gamma_{n}\psi_{n}+\Gamma_{n}\int\mu_{n}\,dQ_{n}-\gamma_{n}+\int\Gamma_{n}\mu_{n}\,dF_{n}\right\\}\right.$
$\displaystyle\qquad\qquad\qquad\qquad\left.-\left\\{\Gamma_{0}\psi_{\infty}+\Gamma_{0}\int\mu_{\infty}\,dQ_{0}-\gamma_{0}+\int\Gamma_{0}\mu_{\infty}\,dF_{0}\right\\}\right)$
$\displaystyle\qquad\qquad+\left\\{\mathbb{P}_{n}(\Gamma_{n}\xi_{n})-P_{0}(\Gamma_{0}\theta_{0})+P_{0}\left(\Gamma_{n}\psi_{n}+\Gamma_{n}\int\mu_{n}\,dQ_{n}-\gamma_{n}+\int\Gamma_{n}\mu_{n}\,dF_{n}-\iint\Gamma_{n}\mu_{n}\,dF_{n}\,dQ_{n}\right)\right\\}$
$\displaystyle\qquad=(\mathbb{P}_{n}-P_{0})(\Gamma_{n}\psi_{n}-\Gamma_{0}\psi_{\infty})+\left\\{-(\mathbb{P}_{n}-P_{0})\gamma_{n}+\mathbb{P}_{n}(\Gamma_{n}\xi_{n})-P_{0}\gamma_{n}\right\\}+\left\\{(\mathbb{P}_{n}-P_{0})\gamma_{0}+P_{0}(\Gamma_{n}\theta_{0})-P_{0}(\Gamma_{0}\theta_{0})\right\\}$
$\displaystyle\qquad\qquad+\iint\Gamma_{n}\mu_{n}\left\\{dQ_{n}\,d(F_{n}-F_{0})+dF_{n}\,d(Q_{n}-Q_{0})+dQ_{n}\,dF_{0}+dF_{n}\,dQ_{0}-dF_{n}\,dQ_{n}-dQ_{0}\,dF_{0}\right\\}$
$\displaystyle\qquad\qquad-\iint\Gamma_{0}\mu_{\infty}\left\\{dQ_{0}\,d(F_{n}-F_{0})+dF_{0}\,d(Q_{n}-Q_{0})\right\\}$
$\displaystyle\qquad\qquad+\left\\{P_{0}(\Gamma_{n}\psi_{n})+\iint\Gamma_{n}\mu_{n}\,dQ_{0}\,dF_{0}-P_{0}(\Gamma_{n}\theta_{0})\right\\}$
(8)

We note that we added and subtracted the terms
$\iint\Gamma_{n}\mu_{n}\,dQ_{0}\,dF_{0}$ and $P_{0}\Gamma_{n}\theta_{0}$
above. We address each of the summands in equation (8) in turn. First, we
write

$(\mathbb{P}_{n}-P_{0})(\Gamma_{n}\psi_{n}-\Gamma_{0}\psi_{\infty})=(\mathbb{P}_{n}-P_{0})\left\\{(\Gamma_{n}-\Gamma_{0})\psi_{n}\right\\}+(\mathbb{P}_{n}-P_{0})\left\\{\Gamma_{0}(\psi_{n}-\psi_{\infty})\right\\}.$

Next, we note that since
$\mathbb{P}_{n}\gamma_{n}=\mathbb{P}_{n}(\Gamma_{n}\xi_{n})$,

$-(\mathbb{P}_{n}-P_{0})\gamma_{n}+\mathbb{P}_{n}(\Gamma_{n}\xi_{n})-P_{0}\gamma_{n}=0.$

Next, we have

$\displaystyle(\mathbb{P}_{n}-P_{0})$
$\displaystyle\gamma_{0}+P_{0}(\Gamma_{n}\theta_{0})-P_{0}(\Gamma_{0}\theta_{0})$
$\displaystyle=P_{0}\left\\{e_{1}^{T}\mathbf{D}^{-1}_{0,h}(\mathbf{D}_{n,h}-\mathbf{D}_{0,h})\mathbf{D}^{-1}_{0,h}w_{h,a_{0},1}K_{h,a_{0}}\theta_{0}\right.$
$\displaystyle\qquad\qquad\left.-e_{3}^{T}c_{2}\tau_{n}^{2}\mathbf{D}^{-1}_{0,b}(\mathbf{D}_{n,b}-\mathbf{D}_{0,b})\mathbf{D}^{-1}_{0,b}w_{b,a_{0},3}K_{b,a_{0}}\theta_{0}\right\\}$
$\displaystyle\qquad+P_{0}\left\\{e_{1}^{T}(\mathbf{D}^{-1}_{n,h}-\mathbf{D}^{-1}_{0,h})w_{h,a_{0},1}K_{h,a_{0}}\theta_{0}-e_{3}^{T}c_{2}\tau_{n}^{2}(\mathbf{D}^{-1}_{n,b}-\mathbf{D}^{-1}_{0,b})w_{b,a_{0},3}K_{b,a_{0}}\theta_{0}\right\\}.$

Now we note that
$\mathbf{D}_{0}^{-1}(\mathbf{D}_{n}-\mathbf{D}_{0})\mathbf{D}_{0}^{-1}=\mathbf{D}_{0}^{-1}\mathbf{D}_{n}(\mathbf{D}_{0}^{-1}-\mathbf{D}_{n}^{-1}),$
so the above equals

$\displaystyle P_{0}$
$\displaystyle\left\\{e_{1}^{T}\mathbf{D}^{-1}_{0,h}\mathbf{D}_{n,h}(\mathbf{D}^{-1}_{0,h}-\mathbf{D}^{-1}_{n,h})w_{h,a_{0},1}K_{h,a_{0}}\theta_{0}\right.$
$\displaystyle\qquad\qquad\left.-e_{3}^{T}c_{2}\tau_{n}^{2}\mathbf{D}^{-1}_{0,b}\mathbf{D}_{n,b}(\mathbf{D}^{-1}_{0,b}-\mathbf{D}^{-1}_{n,b})w_{b,a_{0},3}K_{b,a_{0}}\theta_{0}\right\\}$
$\displaystyle\qquad+P_{0}\left\\{e_{1}^{T}(\mathbf{D}^{-1}_{n,h}-\mathbf{D}^{-1}_{0,h})w_{h,a_{0},1}K_{h,a_{0}}\theta_{0}-e_{3}^{T}c_{2}\tau_{n}^{2}(\mathbf{D}^{-1}_{n,b}-\mathbf{D}^{-1}_{0,b})w_{b,a_{0},3}K_{b,a_{0}}\theta_{0}\right\\}$
$\displaystyle=P_{0}\left\\{e_{1}^{T}(I-\mathbf{D}^{-1}_{0,h}\mathbf{D}_{n,h})(\mathbf{D}^{-1}_{n,h}-\mathbf{D}^{-1}_{0,h})w_{h,a_{0},1}K_{h,a_{0}}\theta_{0}\right.$
$\displaystyle\qquad\qquad\left.-e_{3}^{T}c_{2}\tau_{n}^{2}(I-\mathbf{D}^{-1}_{0,b}\mathbf{D}_{n,b})(\mathbf{D}^{-1}_{n,b}-\mathbf{D}^{-1}_{0,b})w_{b,a_{0},3}K_{b,a_{0}}\theta_{0}\right\\}$
$\displaystyle=R_{n,h,b,a_{0},6}.$

We address the next two terms in equation (8) together:

$\displaystyle\iint\Gamma_{n}\mu_{n}$
$\displaystyle\left\\{dQ_{n}\,d(F_{n}-F_{0})+dF_{n}\,d(Q_{n}-Q_{0})+dQ_{n}\,dF_{0}+dF_{n}\,dQ_{0}-dF_{n}\,dQ_{n}-dQ_{0}\,dF_{0}\right\\}$
$\displaystyle\qquad-\iint\Gamma_{0}\mu_{\infty}\left\\{dQ_{0}\,d(F_{n}-F_{0})+dF_{0}\,d(Q_{n}-Q_{0})\right\\}$
$\displaystyle=\iint\Gamma_{n}\mu_{n}\left\\{dQ_{n}\,dF_{n}-dQ_{0}\,dF_{0}\right\\}-\iint\Gamma_{0}\mu_{n}\left\\{dQ_{0}\,d(F_{n}-F_{0})+dF_{0}\,d(Q_{n}-Q_{0})\right\\}$
$\displaystyle\qquad+\iint\Gamma_{0}(\mu_{n}-\mu_{\infty})\left\\{dQ_{0}\,d(F_{n}-F_{0})+dF_{0}\,d(Q_{n}-Q_{0})\right\\}$
$\displaystyle=\iint\Gamma_{n}\mu_{n}\left\\{dQ_{n}\,dF_{n}-dQ_{0}\,dF_{0}-dQ_{0}\,d(F_{n}-F_{0})-dF_{0}\,d(Q_{n}-Q_{0})\right\\}$
$\displaystyle\qquad+\iint\Gamma_{0}(\mu_{n}-\mu_{\infty})\left\\{dQ_{0}\,d(F_{n}-F_{0})+dF_{0}\,d(Q_{n}-Q_{0})\right\\}$
$\displaystyle\qquad+\iint(\Gamma_{n}-\Gamma_{0})\mu_{n}\left\\{dQ_{0}\,d(F_{n}-F_{0})+dF_{0}\,d(Q_{n}-Q_{0})\right\\}.$

We then note that

$\displaystyle\iint\Gamma_{n}\mu_{n}$
$\displaystyle\left\\{dQ_{n}\,dF_{n}-dQ_{0}\,dF_{0}-dQ_{0}\,d(F_{n}-F_{0})-dF_{0}\,d(Q_{n}-Q_{0})\right\\}$
$\displaystyle=\iint\Gamma_{n}\mu_{n}\left\\{dQ_{n}\,dF_{n}-dQ_{0}\,dF_{n}-dF_{0}\,dQ_{n}+dQ_{0}\,dF_{0}\right\\}$
$\displaystyle=\iint\Gamma_{n}\mu_{n}\,d(Q_{n}-Q_{0})\,d(F_{n}-F_{0})$
$\displaystyle=R_{n,h,b,a_{0},5}.$

Finally, we have

$\displaystyle
P_{0}(\Gamma_{n}\psi_{n})+\iint\Gamma_{n}\mu_{n}\,dQ_{0}\,dF_{0}-P_{0}(\Gamma_{n}\theta_{0})$
$\displaystyle=P_{0}\left(\Gamma_{n}\frac{\mu_{0}-\mu_{n}}{g_{n}}\right)+\iint\Gamma_{n}(\mu_{n}-\mu_{0})\,dQ_{0}\,dF_{0}$
$\displaystyle=\iint\Gamma_{n}\left\\{\frac{g_{0}}{g_{n}}(\mu_{0}-\mu_{n})+(\mu_{n}-\mu_{0})\right\\}\,dQ_{0}\,dF_{0}$
$\displaystyle=\iint\Gamma_{n}\left(1-\frac{g_{0}}{g_{n}}\right)(\mu_{n}-\mu_{0})\,dQ_{0}\,dF_{0}$
$\displaystyle=R_{n,h,b,a_{0},4}.$

Putting these derivations back into equation (8), we have

$\displaystyle(\mathbb{P}_{n}-P_{0})$
$\displaystyle\left(\phi_{n,h,b,a_{0}}^{*}-\phi_{\infty,h,b,a_{0}}^{*}\right)+\left\\{\theta_{n,h,b}(a_{0})-\theta_{0,h,b}(a_{0})+P_{0}\phi_{n,h,b,a_{0}}^{*}\right\\}$
$\displaystyle=(\mathbb{P}_{n}-P_{0})\left\\{(\Gamma_{n}-\Gamma_{0})\psi_{n}\right\\}+(\mathbb{P}_{n}-P_{0})\left\\{\Gamma_{0}(\psi_{n}-\psi_{\infty})\right\\}$
$\displaystyle\qquad+\iint\Gamma_{0}(\mu_{n}-\mu_{\infty})\left\\{dQ_{0}\,d(F_{n}-F_{0})+dF_{0}\,d(Q_{n}-Q_{0})\right\\}$
$\displaystyle\qquad+\iint(\Gamma_{n}-\Gamma_{0})\mu_{n}\left\\{dQ_{0}\,d(F_{n}-F_{0})+dF_{0}\,d(Q_{n}-Q_{0})\right\\}$
$\displaystyle\qquad+R_{n,h,b,a_{0},4}+R_{n,h,b,a_{0},5}+R_{n,h,b,a_{0},6}$
$\displaystyle=(\mathbb{P}_{n}-P_{0})\left[\left\\{\Gamma_{0}\left(\psi_{n}+\int\mu_{n}\,dQ_{0}\right)+\int\Gamma_{0}\mu_{n}\,dF_{0}\right\\}-\left\\{\Gamma_{0}\left(\psi_{\infty}+\int\mu_{\infty}\,dQ_{0}\right)+\int\Gamma_{0}\mu_{\infty}\,dF_{0}\right\\}\right]$
$\displaystyle\qquad+(\mathbb{P}_{n}-P_{0})\left\\{(\Gamma_{n}-\Gamma_{0})\left(\psi_{n}+\int\mu_{n}\,dQ_{0}\right)+\int(\Gamma_{n}-\Gamma_{0})\mu_{n}\,dF_{0}\right]$
$\displaystyle\qquad+R_{n,h,b,a_{0},4}+R_{n,h,b,a_{0},5}+R_{n,h,b,a_{0},6}$
$\displaystyle=\sum_{j=2}^{6}R_{n,h,b,a_{0},j}.$

∎

##  Supplement E Proof of Theorems

###### Proof of Theorem 2.

By Lemma 3, we have

$\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})=\mathbb{P}_{n}\phi_{\infty,h,b,a_{0}}^{*}+\sum_{j=1}^{6}R_{n,h,b,a_{0},j}.$

By Lemma 7, $R_{n,h,b,a_{0},1}=o\left(h^{2}\right)$. By Corollary 4 and Lemmas
17, 19, and 20, $R_{n,h,b,a_{0},j}=o_{p}\left(\\{nh\\}^{-1/2}\right)$ for each
$j\in\\{1,\dotsc,6\\}$. This establishes the first claim. By Lemma 5,
$(nh)^{1/2}\mathbb{P}_{n}\phi_{\infty,h,b,a_{0}}^{*}$ converges in
distribution to the claimed limit distribution. If $nh^{5}=O(1)$, then
$(nh)^{1/2}R_{n,h,b,a_{0},1}=o\left(\\{nh^{5}\\}^{1/2}\right)=o(1)$, so

$(nh)^{1/2}\left[\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})\right]=(nh)^{1/2}\mathbb{P}_{n}\phi_{\infty,h,b,a_{0}}^{*}+o_{p}(1),$

so $(nh)^{1/2}\left[\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})\right]$ converges
to this same limit. For the final statement, by Lemma 21,
$\sigma_{n,h,b}^{2}(a_{0})=h\mathbb{P}_{n}(\phi_{n,h,b,a_{0}}^{*})^{2}$
satisfies
$\sigma_{n,h,b}^{2}(a_{0})-hP_{0}(\phi_{\infty,h,b,a_{0}}^{*})^{2}=o_{p}(1)$.
By Lemma 5, $hP_{0}(\phi_{\infty,h,b,a_{0}}^{*})^{2}$ converges to the
variance of the limit distribution,
$V_{K,\tau}\sigma_{0}^{2}(a_{0})/f_{0}(a_{0})$. Hence, $\sigma_{n,h,b}(a_{0})$
converges in probability to the standard deviation of the limit distribution,
so the final statement follows by Slutsky’s theorem. ∎

###### Proof of Theorem 3.

By Theorem 1 and since $nh^{5}=O(1)$, we have

$(nh)^{1/2}\begin{pmatrix}\theta_{n,h,b}(a_{1})-\theta_{0}(a_{1})\\\ \vdots\\\
\theta_{n,h,b}(a_{m})-\theta_{0}(a_{m})\end{pmatrix}=(nh)^{1/2}\mathbb{P}_{n}\begin{pmatrix}\phi_{\infty,h,b,a_{1}}^{*}\\\
\vdots\\\ \phi_{\infty,h,b,a_{m}}^{*}\end{pmatrix}+o_{p}(1).$

The result follows by Lemma 6 and Slutsky’s theorem. ∎

###### Proof of Theorem 4.

By the triangle inequality, we first have

$\displaystyle\sup_{t\in\mathbb{R}}\left|P_{0}\left(\sup_{a_{0}\in\mathcal{A}_{0}}(nh)^{1/2}\left|\frac{\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})}{\sigma_{n,h,b}(a_{0})}\right|\leq
t\right)-P_{0}\left(\max_{a_{0}\in\mathcal{A}_{n}}\left|Z_{n,h,b}(a_{0})\right|\leq
t\mid\mathbf{O}_{n}\right)\right|$
$\displaystyle\qquad\leq\sup_{t\in\mathbb{R}}\left|P_{0}\left(\sup_{a_{0}\in\mathcal{A}_{0}}(nh)^{1/2}\left|\frac{\theta_{n,h,b}(a_{0})-\theta_{0}(a_{0})}{\sigma_{n,h,b}(a_{0})}\right|\leq
t\right)-P_{0}\left(\sup_{a_{0}\in\mathcal{A}_{0}}\left|Z_{\infty,h,b}(a_{0})\right|\leq
t\right)\right|$
$\displaystyle\qquad\qquad+\sup_{t\in\mathbb{R}}\left|P_{0}\left(\sup_{a_{0}\in\mathcal{A}_{0}}\left|Z_{\infty,h,b}(a_{0})\right|\leq
t\right)-P_{0}\left(\max_{a_{0}\in\mathcal{A}_{n}}\left|Z_{\infty,h,b}(a_{0})\right|\leq
t\right)\right|$
$\displaystyle\qquad\qquad+\sup_{t\in\mathbb{R}}\left|P_{0}\left(\max_{a_{0}\in\mathcal{A}_{n}}\left|Z_{\infty,h,b}(a_{0})\right|\leq
t\right)-P_{0}\left(\max_{a_{0}\in\mathcal{A}_{n}}\left|Z_{n,h,b}(a_{0})\right|\leq
t\mid\mathbf{O}_{n}\right)\right|.$

The first term on the right hand side is $o(1)$ by Lemma 27. The second term
on the right hand side is $o(1)$ by Lemma 28. The last term on the right hand
side is $o_{p}(1)$ by Lemma 29. ∎

##  Supplement F Analysis of the leading term

We first establish the following result on the behavior of
$\mathbf{D}^{-1}_{0,h,a_{0},1}$ and $\mathbf{D}^{-1}_{0,b,a_{0},3}$ as
$h,b\longrightarrow 0$.

######  Lemma 4.

If conditions (A1) and (A5) hold, then
$\mathbf{D}^{-1}_{0,h,a_{0},1}=f_{0}(a_{0})^{-1}\mathbf{S}_{2}^{-1}+O(h)$ as
$h\longrightarrow 0$ and
$\mathbf{D}^{-1}_{0,b,a_{0},3}=f_{0}(a_{0})^{-1}\mathbf{S}_{4}^{-1}+O(b)$ as
$b\longrightarrow 0$. If conditions (A1) and (A8) hold, then
$\sup_{a_{0}\in\mathcal{A}_{0}}\left\|\mathbf{D}^{-1}_{0,h,a_{0},1}-f_{0}(a_{0})^{-1}\mathbf{S}_{2}^{-1}\right\|_{\infty}=O(h)$
and
$\sup_{a_{0}\in\mathcal{A}_{0}}\left\|\mathbf{D}^{-1}_{0,h,a_{0},3}-f_{0}(a_{0})^{-1}\mathbf{S}_{4}^{-1}\right\|_{\infty}=O(h)$
as $h\longrightarrow 0$.

###### Proof of Lemma 4.

The proof is analogous to that of Section 3.7 of Fan and Gijbels, (1996) and
many others. By definition of $\mathbf{D}_{0,h,a_{0},1}$ and the change of
variables $u=(t-a_{0})/h$,

$\displaystyle\mathbf{D}_{0,h,a_{0},1}[j,k]$
$\displaystyle=\int\left(\frac{t-a_{0}}{h}\right)^{j+k-2}h^{-1}K\left(\frac{t-a_{0}}{h}\right)\,dF_{0}(t)=\int