# Nonparametric Inference on Dose-Response Curves Without the Positivity
Condition

Yikun Zhanglabel=e1 [<EMAIL_ADDRESS>Yen-Chi Chenlabel=e2 [
<EMAIL_ADDRESS>Alexander Giessinglabel=e3 [<EMAIL_ADDRESS>Department of Statistics, University of Washington

# Supplement to “Nonparametric Inference on Dose-Response Curves Without the
Positivity Condition”

Yikun Zhanglabel=e1 [<EMAIL_ADDRESS>Yen-Chi Chenlabel=e2 [
<EMAIL_ADDRESS>Alexander Giessinglabel=e3 [<EMAIL_ADDRESS>Department of Statistics, University of Washington

###### Abstract

Existing statistical methods in causal inference often rely on the assumption
that every individual has some chance of receiving any treatment level
regardless of its associated covariates, which is known as the positivity
condition. This assumption could be violated in observational studies with
continuous treatments. In this paper, we present a novel integral estimator of
the causal effects with continuous treatments (_i.e._ , dose-response curves)
without requiring the positivity condition. Our approach involves estimating
the derivative function of the treatment effect on each observed data sample
and integrating it to the treatment level of interest so as to address the
bias resulting from the lack of positivity condition. The validity of our
approach relies on an alternative weaker assumption that can be satisfied by
additive confounding models. We provide a fast and reliable numerical recipe
for computing our estimator in practice and derive its related asymptotic
theory. To conduct valid inference on the dose-response curve and its
derivative, we propose using the nonparametric bootstrap and establish its
consistency. The practical performances of our proposed estimators are
validated through simulation studies and an analysis of the effect of air
pollution exposure (PM2.5) on cardiovascular mortality rates.

62G05,

62D20, 62G20,

Causal inference,

dose-response curve,

positivity,

kernel smoothing,

###### keywords:

[class=MSC2020]

###### keywords:

, and

##  1 Introduction

In observational studies, the causal effect of interest does not always result
from a standard binary intervention but rather comes as a consequence of a
continuous treatment or exposure. Such a causal effect on the outcome variable
$Y\in\mathcal{Y}\subset\mathbb{R}$ from a continuous treatment variable
$T\in\mathcal{T}\subset\mathbb{R}$ is known as the (causal) dose-response
curve or relationship due to its major application in studying biological
effects (Waud, 1975). More precisely, a dose-response curve characterizes the
average outcome if all units would have been assigned to a certain treatment
level $T=t$. In practice, it is fairly common that an additional set of
covariates $\bm{S}\in\mathcal{S}\subset\mathbb{R}^{d}$ is collected, which, to
some extent, contains all the possible confounding variables that influence
both the treatment $T$ and the outcome $Y$. Under regularity conditions (see
Section 2.2), the dose-response curve coincides with the so-called covariate-
adjusted regression function $t\mapsto
m(t)=\mathbb{E}\left[\mu(t,\bm{S})\right]$ with
$\mu(t,\bm{s})=\mathbb{E}\left(Y|T=t,\bm{S}=\bm{s}\right)$ and is thus
identifiable from the observed data
$\left\\{(Y_{i},T_{i},\bm{S}_{i})\right\\}_{i=1}^{n}$ for any
$t\in\mathcal{T}$ (Robins et al., 2000; Neugebauer and van der Laan, 2007;
Díaz and van der Laan, 2013; Kennedy et al., 2017). Alternatively, the
covariate-adjusted regression function $m(t)$ can be viewed as a univariate
summary function of the outcome against a continuous covariate $T$ when we
average the regression function $\mu(t,\bm{s})$ over all other covariates
$\bm{S}$ (Takatsu and Westling, 2022).

However, the regularity conditions for identifying and estimating the dose-
response curve might not be verifiable or could even be violated in
observational studies. In particular, it is commonly assumed that there is a
sufficient amount of variability in the treatment assignment within each
strata of the covariates, which is captured by the following positivity or
overlapping condition.

######  Assumption A0 (Positivity).

The conditional density $p(t|\bm{s})$ is bounded above and away from zero
almost surely for all $t\in\mathcal{T}$ and $\bm{s}\in\mathcal{S}$.

The positivity condition (A0) may be violated in observational studies for
either of the following two reasons: (i) theoretically, it is impossible for
some individuals with certain covariate values to receive some levels of the
treatment, and (ii) practically, individuals at some levels of the treatment
may not be collected in the finite data sample; see Cole and Hernán (2008);
Westreich and Cole (2010); Petersen et al. (2012) for related discussions.
These problems are particularly pervasive under the context of continuous
treatments. For instance, air pollution levels can vary greatly across larger
regions but remain relatively consistent within smaller and nearby areas.
Therefore, individuals at the same location are typically only get exposed to
one level of air pollution, _i.e._ , spatial confounding variables may change
at a finer scale than the variation of exposure, thereby violating the
positivity condition (Paciorek, 2010; Schnell and Papadogeorgou, 2020; Keller
and Szpiro, 2020). As a more concrete example, consider a single confounder
model

$Y=T^{2}+T+1+10S+\epsilon,\quad T=\sin(\pi S)+E,\quad\text{ and }\quad
S\sim\text{Uniform}[-1,1]\subset\mathbb{R},$ (1)

where $E\sim\text{Uniform}[-0.3,0.3]$ is an independent treatment variation
and $\epsilon\sim\mathcal{N}(0,1)$ is an exogenous normal noise. The marginal
supports of $T$ and $S$ are $\mathcal{T}=[-1.3,1.3]$ and $\mathcal{S}=[-1,1]$
respectively, while the joint support of $(T,S)$ only covers a thin band
region of the product space $\mathcal{T}\times\mathcal{S}$; see the left panel
of Fig 1. The conditional density $p(t|s)$ for any $s\in\mathcal{S}$ is 0
within the gray regions, and the positivity condition (A0) clearly fails.
Without (A0), the existing approaches for estimating the dose-response curve
$m(t)$ and its derivative $\theta(t)=m^{\prime}(t)$ can be very unstable at
some specific treatment levels; see the usual regression adjustment estimators
(defined in (6) and Remark 1 below) in the middle and right panels of Fig 1
for illustrations.

Figure 1: Simulation results under the single confounder model (1). Left: The
support of the joint distribution of $(T,S)$. Middle: The estimated dose-
response curves by usual regression adjustment and our proposed estimators
overlain with the true one $m(t)$. Right: The estimated derivatives of the
dose-response curve by usual regression adjustment and our proposed estimators
overlain with the true one $\theta(t)$. The middle and right panels also
present the 95% confidence intervals and/or uniform confidence bands from our
proposed estimators as shaded regions and dashed lines, respectively.

In this paper, we propose a novel integral estimator that can consistently
recover the entire dose-response curve $m(t)$ and a localized estimator of its
derivative $\theta(t)=m^{\prime}(t)$ even when the positivity condition (A0)
fails to hold in some regions of the product space
$\mathcal{T}\times\mathcal{S}$; see the middle and right panels of Fig 1. Our
main contributions are summarized as follows:

1\. Methodology: After discussing the identification conditions of $m(t)$ as a
causal dose-response curve in Section 2, we introduce our integral estimator
of $m(t)$ in Section 3.1, which is constructed from a localized estimator of
$\theta(t)$ around the observed data and extrapolate to any treatment level
$t\in\mathcal{T}$ through the fundamental theorem of calculus. Such an
integral estimator can also be efficiently computed in practice via the
Riemann sum approximation in Section 3.2 and be reliably inferred through
nonparametric bootstrap in Section 3.3.

2\. Asymptotic Theory: We establish uniform consistencies of our proposed
integral estimator and localized derivative estimator under the context of
kernel smoothing methods in Section 4.2. We also prove the validity of
nonparametric bootstrap inference in Section 4.3.

3\. Experiments: We demonstrate the finite-sample performances of our proposed
estimators through simulation studies and a case study of the effect of fine
particulate matter (PM2.5) on cardiovascular mortality rates in Section 5. All
the code for our experiments is available at
https://github.com/zhangyk8/NPDoseResponse/tree/main/Paper_Code.

###  1.1 Other Related Works

Estimating the dose-response curve $m(t)$ is a technically difficult problem
in causal inference due to the fact that $m(t)$ is not pathwise differentiable
and cannot be consistently estimated in a $\sqrt{n}$ rate (Chamberlain, 1986;
van der Vaart, 1991). Parametrically, Robins et al. (2000) pioneered a
marginal structural model to estimate $m(t)$, whose dependence on the correct
specification of a parametric model was later relaxed by van der Laan and
Robins (2003); Neugebauer and van der Laan (2007) through a projection of
$m(t)$ to the parametric model space. Nonparametrically, a regression
adjustment approach to estimating $m(t)=\mathbb{E}\left[\mu(t,\bm{S})\right]$
under a two-stage kernel smoothing estimator was first studied by Newey (1994)
and later adapted to the estimation of dose-response curves by Flores (2007).
However, the consistency and asymptotic normality of their kernel-based
estimators rely on an undersmoothing bandwidth, which is asymptotically
smaller than the one that attains the optimal asymptotic bias-variance trade-
off and is difficult to select in practice; see Section 5.7 in Wasserman
(2006). On the contrary, our proposed integral estimator of $m(t)$ under the
context of kernel smoothing is based on the derivative estimation and require
neither explicit undersmoothing nor bias correction (Calonico et al., 2018;
Takatsu and Westling, 2022), permitting the use of any bandwidth selection
method in nonparametric regression.

There are many existing works about the estimation of average derivative
effects in the literature; see Härdle and Stoker (1989); Powell et al. (1989);
Newey and Stoker (1993); Cattaneo et al. (2010); Hirshberg and Wager (2020)
and references therein. Under some regularity conditions, average derivative
effects are identical to the so-called incremental treatment effects, which
are closely related to the derivatives $\theta(t)=m^{\prime}(t)$ of dose-
response curves; see Proposition 1 in Rothenhäusler and Yu (2019) and Section
6.1 in Hines et al. (2023).

All the aforementioned methods for estimating $m(t)$ and $\theta(t)$ assume
the positivity condition (A0). Under the context of binary treatments, there
are previous researches studying the dependence of common estimators on (A0)
from the perspectives of convergence rates (Khan and Tamer, 2010; D’Amour et
al., 2021) and empirical performances (Busso et al., 2014; Léger et al.,
2022). To address the violation of positivity, some approaches considered
trimming the extreme values of (estimated) propensity scores (_i.e._ ,
discrete versions of the conditional density $p(t|\bm{s})$; Dehejia and Wahba
1999; Crump et al. 2009; Yang and Ding 2018), while others proposed robust
estimation methods against the violation (Rothe, 2017; Ma and Wang, 2020). In
addition, Kennedy (2019) considered switching the causal estimand to the
incremental treatment effect for discrete treatments, which can be
nonparametrically estimated without positivity. Rothenhäusler and Yu (2019)
later studied the same estimand for continuous treatments. To the best of our
knowledge, none of existing works address the nonparametric inference on
$m(t)$ and $\theta(t)$ without assuming (A0).

###  1.2 Notations

Throughout the paper, we consider a real-valued outcome variable
$Y\in\mathcal{Y}\subset\mathbb{R}$, univariate continuous treatment
$T\in\mathcal{T}\subset\mathbb{R}$, and a vector of covariates
$\bm{S}=(S_{1},...,S_{d})\in\mathcal{S}\subset\mathbb{R}^{d}$ with a fixed
dimension $d$. We write $Y{\perp\\!\\!\\!\perp}X$ when the random variables
$Y,X$ are independent. The data sample consists of independent and identically
distributed (i.i.d.) observations
$\bm{U}_{i}=(Y_{i},T_{i},\bm{S}_{i}),i=1,...,n$ with the common distribution
$\mathrm{P}$ and Lebesgue density $p(y,t,\bm{s})=p(y|t,\bm{s})\cdot
p(t|\bm{s})\cdot p_{S}(\bm{s})$. Here, $p_{T}(t)$ and $p_{S}(\bm{s})$ are the
marginal densities of $T$ and $\bm{S}$, respectively, and
$p(t|\bm{s})=\frac{\partial}{\partial t}\mbox{$\mathrm{P}$}\left(T\leq
t|\bm{S}=\bm{s}\right)$ is the conditional density of $T$ given covariates
$\bm{S}=\bm{s}$. We also denote the joint density of $(T,\bm{S})$ by
$p(t,\bm{s})=p(t|\bm{s})\cdot p_{S}(\bm{s})=p(\bm{s}|t)\cdot p_{T}(t)$. Let
$\mathbb{P}_{n}$ be the empirical distribution of the observed data so that
$\mathbb{P}_{n}g=\int
g(\bm{u})\,d\mathbb{P}_{n}(\bm{u})=\frac{1}{n}\sum_{i=1}^{n}g(\bm{U}_{i})$ for
some measurable function $g$. At the population level, we denote
$\mbox{$\mathrm{P}$}g=\mathbb{E}\left[g(\bm{U})\right]=\int
g(\bm{u})\,d\mbox{$\mathrm{P}$}(\bm{u})$ for $\mathrm{P}$-integrable function
$g$. In addition, we define the empirical process evaluated at a
$\mathrm{P}$-integrable function $g$ as
$\mathbb{G}_{n}(g)=\sqrt{n}\left(\mathbb{P}_{n}-\mbox{$\mathrm{P}$}\right)g=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left[g(\bm{U}_{i})-\mathbb{E}\left(g(\bm{U}_{i})\right)\right]$.
Finally, we use $\mathbbm{1}_{A}$ to denote the indicator function of a set
$A$.

We use the big-$O$ notation $h_{n}=O(g_{n})$ if $|h_{n}|$ is upper bounded by
a positive constant multiple of $g_{n}>0$ when $n$ is sufficiently large. In
contrast, $h_{n}=o(g_{n})$ when $\lim_{n\to\infty}\frac{|h_{n}|}{g_{n}}=0$.
For random variables, the notation $o_{P}(1)$ is short for a sequence of
random variables converging to zero in probability, while the expression
$O_{P}(1)$ denotes the sequence that is bounded in probability. We also use
the notation $a_{n}\lesssim b_{n}$ or $b_{n}\gtrsim a_{n}$ when there exists
an absolute constant $A>0$ such that $a_{n}\leq Ab_{n}$ when $n$ is large. If
$a_{n}\gtrsim b_{n}$ and $a_{n}\lesssim b_{n}$, then $a_{n},b_{n}$ are
asymptotically equal and it is denoted by $a_{n}\asymp b_{n}$.

##  2 Model Setups and Identification Conditions

We assume that $\bm{U}_{i}=(Y_{i},T_{i},\bm{S}_{i}),i=1,...,n$ are i.i.d.
observations generated from the following model:

$\displaystyle\begin{split}Y&=\mu(T,\bm{S})+\epsilon\quad\text{ and }\quad
T=f(\bm{S})+E,\end{split}$ (2)

where $E\in\mathbb{R}$ is the treatment variation with $\mathbb{E}(E)=0$,
$E{\perp\\!\\!\\!\perp}\bm{S}$, and $\epsilon\in\mathbb{R}$ is an exogenous
noise variable with
$\epsilon{\perp\\!\\!\\!\perp}E,\epsilon{\perp\\!\\!\\!\perp}\bm{S}$ and
$\mathbb{E}(\epsilon)=0,\mathrm{Var}(\epsilon)=\sigma^{2}>0,\mathbb{E}(\epsilon^{4})<\infty$.
The function $f$ determines how the covariates (or confounding variables)
$\bm{S}$ influence the treatment $T$. Under model (2), the covariate-adjusted
regression function (or the dose-response curve under identification
conditions in Section 2.2) is given by

$t\mapsto
m(t)=\mathbb{E}\left[\mathbb{E}(Y|T=t,\bm{S})\right]=\mathbb{E}\left[\mu(t,\bm{S})\right],$

and its derivative is written as
$\theta(t)=m^{\prime}(t)=\frac{d}{dt}\mathbb{E}\left[\mu(t,\bm{S})\right]$.
Notice that the derivative function $\theta(t)$ can be viewed as the
continuous version of the average treatment effect in causal inference. In the
presence of confounding variables $\bm{S}$, $m(t)$ cannot be identified by
regressing $Y$ only with respect to $T$, because
$m(t)=\mathbb{E}\left[\mu(t,\bm{S})\right]\neq\mathbb{E}\left(Y|T=t\right)=\mathbb{E}\left[\mu(T,\bm{S})|T=t\right]$.

###  2.1 Motivating Example: Additive Confounding Model

One important exemplification of model (2) that we will frequently refer to in
this paper is the following additive confounding model:

$\displaystyle\begin{split}Y&=\bar{m}(T)+\eta(\bm{S})+\epsilon\quad\text{ and
}\quad T=f(\bm{S})+E,\end{split}$ (3)

where $\bar{m}(t)$ is the primary treatment effect of interest and
$\eta(\bm{s})$ is the random effect with
$\mbox{$\mathbb{E}$}\left[\eta(\bm{S})\right]=0$. Such an additive form is a
common working model in spatial confounding problems (Paciorek, 2010; Schnell
and Papadogeorgou, 2020) and also known as the geoadditive structural equation
model (Kammann and Wand, 2003; Thaden and Kneib, 2018; Wiecha and Reich,
2024), where $\bm{S}\in\mathbb{R}^{d}$ are the spatial locations (usually with
$d=2$) or other spatially correlated covariates that affect both the treatment
$T$ and the outcome $Y$. We summarize key properties of the additive
confounding model (3) in the following proposition.

######  Proposition 1 (Properties of the additive confounding model (3)).

Let $\theta_{M}(t)=\mathbb{E}\left[\frac{\partial}{\partial
t}\mu(t,\bm{S})\right]$ and
$\theta_{C}(t)=\mathbb{E}\left[\frac{\partial}{\partial
t}\mu(t,\bm{S})\big{|}T=t\right]$, where
$\mu(t,\bm{s})=\mathbb{E}\left[Y|T=t,\bm{S}=\bm{s}\right]$. Under the additive
confounding model (3) with $\mbox{$\mathbb{E}$}\left[\eta(\bm{S})\right]=0$,
the following results hold for all $t\in\mathcal{T}$:

  1. (a)

$m(t)=\bar{m}(t)$.

  2. (b)

$m(t)\neq\mbox{$\mathbb{E}$}(Y|T=t)=\bar{m}(t)+\mbox{$\mathbb{E}$}\left[\eta(\bm{S})|T=t\right]$
when $\mbox{$\mathbb{E}$}\left[\eta(\bm{S})|T=t\right]\neq 0$.

  3. (c)

$\theta(t)=\theta_{M}(t)=\theta_{C}(t)$.

  4. (d)

$\theta(t)\neq\frac{d}{dt}\mbox{$\mathbb{E}$}\left[\mu(t,\bm{S})|T=t\right]=\theta(t)+\frac{d}{dt}\mbox{$\mathbb{E}$}\left[\eta(\bm{S})|T=t\right]$
when $\frac{d}{dt}\mbox{$\mathbb{E}$}\left[\eta(\bm{S})|T=t\right]\neq 0$.

  5. (e)

$\mathbb{E}\left[\mu(T,\bm{S})\right]=\mathbb{E}\left[m(T)\right]$ even when
$\mbox{$\mathbb{E}$}\left[\eta(\bm{S})\right]\neq 0$.

The above results hold even if the treatment variation $E=0$ almost surely.

The proof of Proposition 1 is in Section B.1. Note that in Proposition 1(e),
the first expectation is taken with respect to the joint distribution of
$(T,\bm{S})$ while the second expectation is taken with respect to the
marginal distribution of $T$.

###  2.2 Identification Conditions

While $m(t)$ can be defined through model (2), we need some additional
assumptions in order to express $m(t)$ as a causal effect of $T$ on $Y$ and
identify it from the observed data
$\left\\{(Y_{i},T_{i},\bm{S}_{i})\right\\}_{i=1}^{n}$. Following the potential
outcome framework (Rubin, 1974), we let $Y(t,\bm{s})$ denote the potential
outcome that would have been observed under treatment level $T=t$ and
covariate vector $\bm{S}=\bm{s}$. We summarize the required assumptions as
follows.

######  Assumption A1 (Identification conditions for $m(t)$).

  1. (a)

(Treatment effect and consistency) $Y(t,\bm{s})=Y(t)=Y$ for any
$t\in\mathcal{T}$ and $\bm{s}\in\mathcal{S}$.

  2. (b)

(Ignorability or unconfoundingness)
$Y(t){\perp\\!\\!\\!\perp}T\,\big{|}\,\bm{S}$ for all $t\in\mathcal{T}$.

  3. (c)

(Treatment variation) The treatment variation $E$ has nonzero variance, _i.e._
, $\mathrm{Var}(E)>0$.

Assumption A1(a) is a continuous version of the stable unit treatment value
assumption (SUTVA; Page 19 of Cox 1958 and Rubin 1980). The first equality
ensures that the potential outcome $Y(t)$ under treatment level $T=t$ is well-
defined whatever value the covariate vector $\bm{S}$ takes, while the second
equality guarantees that the treatment level of any subject does not affect
the potential outcomes of others (_i.e._ , no interference) and there are no
different versions of treatments. Given that both the treatment $T$ and the
outcome $Y$ are continuous, the essential consistency condition that we need
is the identity of conditional distributions
$\mbox{$\mathrm{P}$}\left(Y(T)\leq
y|T=t\right)=\mbox{$\mathrm{P}$}\left(Y(t)\leq|T=t\right)$ for any
$t\in\mathcal{T}$ and $y\in\mathcal{Y}$; see Gill and Robins (2001) for
details. The ignorability condition (Assumption A1(b)) was first generalized
to continuous treatments by Hirano and Imbens (2004), stating that the
potential outcome variable is independent of the treatment level within any
specific strata of covariates. In the context of spatial confounding, the
ignorability condition holds when the spatial locations are taken into account
(Gilbert et al., 2023). It also implies that the mean potential outcome under
$T=t$ remains the same across all treatment levels when we condition on
$\bm{S}=\bm{s}$. Finally, the treatment variation condition (Assumption A1(c))
is crucial for identifying the conditional mean outcome (or regression)
function $\mu(t,\bm{s})$ on a non-degenerate region of
$\mathcal{T}\times\mathcal{S}$. When $\mathrm{Var}(E)=0$, $\mu(t,\bm{s})$ can
only be identified on the lower dimensional surface
$\left\\{(t,\bm{s})\in\mathcal{T}\times\mathcal{S}:t=f(\bm{s})\right\\}$,
which is also the support of $(T,\bm{S})$. We derive nonparametric bounds on
$m(t)=\mbox{$\mathbb{E}$}\left[\mu(t,\bm{S})\right]$ and its derivative
$\theta(t)=\frac{d}{dt}\mbox{$\mathbb{E}$}\left[\mu(t,\bm{S})\right]$ in
Section A. We also demonstrate how $\mathrm{Var}(E)=0$ can lead to ambiguous
definitions of the associated dose-response curves in Example 1 below.

######  Example 1 (Necessity of $\mathrm{Var}(E)>0$).

Suppose that $\mathrm{Var}(E)=0$ and $T=f(\bm{S})+E=S_{1}$ almost surely,
where $S_{1}$ is the first component of
$\bm{S}\in\mathcal{S}\subset\mathbb{R}^{d}$. We further assume that
$\mathbb{E}(S_{1})=0$. Now, consider two equivalent conditional mean outcome
functions

$\displaystyle\mu_{1}(T,\bm{S})\equiv T+2S_{1}\quad\text{ and
}\quad\mu_{2}(T,\bm{S})\equiv 2T+S_{1},$

both of which are equal to $3S_{1}$ and agree on the support
$\left\\{(t,\bm{s})\in\mathcal{T}\times\mathcal{S}:t=f(\bm{s})=s_{1}\right\\}$.
However, these two conditional mean outcome functions lead to two distinct
treatment effects:

$m_{1}(t)=\mbox{$\mathbb{E}$}\left[\mu_{1}(t,\bm{S})\right]=t\quad\text{ and
}\quad m_{2}(t)=\mbox{$\mathbb{E}$}\left[\mu_{2}(t,\bm{S})\right]=2t,$

whose derivatives are different as well.

Under Assumption A1, the dose-response curve
$t\mapsto\mbox{$\mathbb{E}$}\left[Y(t)\right]$ is identical to the standard
covariate-adjusted regression function
$m(t)=\mbox{$\mathbb{E}$}\left[\mu(t,\bm{S})\right]=\int_{\mathcal{S}}\mu(t,\bm{s})\cdot
p_{S}(\bm{s})\,d\bm{s}$ and can thus be estimated from the observed data
$\left\\{(Y_{i},T_{i},\bm{S}_{i})\right\\}_{i=1}^{n}$. To estimate the
derivative
$\theta(t)=m^{\prime}(t)=\frac{d}{dt}\mbox{$\mathbb{E}$}\left[\mu(t,\bm{S})\right]$
of the dose-response curve from the observed data sample, we impose an
additional assumption on the conditional mean outcome function
$\mu(t,\bm{s})=\mbox{$\mathbb{E}$}(Y|T=t,\bm{S}=\bm{s})$ that enables us to
bypass the positivity condition (A0).

######  Assumption A2 (Identification condition for $\theta(t)$).

The function $\mu(t,\bm{s})$ is continuously differentiable with respect to
$t$ for any $(t,\bm{s})\in\mathcal{T}\times\mathcal{S}$ and the following two
equalities hold true:

$\theta(t)=\mbox{$\mathbb{E}$}\left[\frac{\partial}{\partial
t}\mu(t,\bm{S})\right]=\mbox{$\mathbb{E}$}\left[\frac{\partial}{\partial
t}\mu(t,\bm{S})\Big{|}T=t\right]$ (4)

and

$\mathbb{E}\left[\mu(T,\bm{S})\right]=\mathbb{E}\left[m(T)\right].$ (5)

The first equality in (4) is a mild condition and can be satisfied under
various settings. In particular, it holds when $\left|\frac{\partial}{\partial
t}\mu(t,\bm{S})\right|$ is bounded by an integrable function
$\bar{\mu}(\bm{S})$ with respect to the distribution of $\bm{S}$; see Theorem
1.1 and Example 1.8 in Shao (2003). The second equality in (4) as well as the
equation (5) of Assumption A2 are stricter but still valid under the additive
confounding model in Section 2.1; recall Proposition 1(c,e). Under Assumption
A2, we can express $\theta(t)$ in three different but also equivalent ways

$\theta(t)=\underbrace{\frac{d}{dt}\mathbb{E}\left[\mu(t,\bm{S})\right]}_{m^{\prime}(t)}=\underbrace{\mbox{$\mathbb{E}$}\left[\frac{\partial}{\partial
t}\mu(t,\bm{S})\right]}_{\theta_{M}(t)}=\underbrace{\mbox{$\mathbb{E}$}\left[\frac{\partial}{\partial
t}\mu(t,\bm{S})\Big{|}T=t\right]}_{\theta_{C}(t)}.$

Estimating the derivative effect $\theta(t)$ through the form of
$\theta_{C}(t)$ is our key technique to bypass the particularly strong
positivity condition (A0) under the continuous treatment setting. After
$\theta(t)$ is consistently recovered, we refer the estimation back to the
dose-response curve $m(t)$ by (5) and the fundamental theorem of calculus; see
Section 3.1 for details.

##  3 Nonparametric Inference Without the Positivity Condition

Under the positivity condition (A0), together with Assumption A1, the dose-
response curve $m(t)=\mbox{$\mathbb{E}$}\left[Y(t)\right]$ can be identified
through the covariate-adjusted regression function
$\mbox{$\mathbb{E}$}\left[\mu(t,\bm{S})\right]$ from the observed data
$\left\\{(Y_{i},T_{i},\bm{S}_{i})\right\\}_{i=1}^{n}$, suggesting the
following regression adjustment (RA) or G-computation (Robins, 1986; Gill and
Robins, 2001) estimator as:

$\widehat{m}_{RA}(t)=\frac{1}{n}\sum_{i=1}^{n}\widehat{\mu}(t,\bm{S}_{i}),$
(6)

where $\widehat{\mu}(t,\bm{s})$ is any consistent estimator of the conditional
mean outcome function $\mu(t,\bm{s})$. However, when the positivity condition
(A0) does not hold for some region in $\mathcal{T}\times\mathcal{S}$, the
above estimator (6) will be unstable and even inconsistent. This is because
without (A0), the joint density $p(t,\bm{S}_{i})=p(t|\bm{S}_{i})\cdot
p_{S}(\bm{S}_{i})$ can be close to 0 for some $i=1,...,n$, and $\mu(t,\bm{s})$
cannot be consistently estimated for those query points
$(t,\bm{S}_{i}),i=1,...,n$. Other existing methods for estimating $m(t)$, such
as the inverse probability of treatment weighted (IPTW) or its augmented
variants, also relies on the validity of (A0) for their consistency and
empirical behaviors (Díaz and van der Laan, 2013; Kennedy et al., 2017; Huber
et al., 2020; Colangelo and Lee, 2020). The same issue incurred by the failure
of (A0) also applies to the estimation of $\theta(t)=m^{\prime}(t)$.

In this section, we introduce a novel integral estimator to resolve the
inconsistency issues of existing estimators in recovering the dose-response
curve $m(t)$ and its derivative $\theta(t)$ without the positivity condition
(A0). On one hand, the estimation of $\theta(t)$ is based on the conditional
expectation
$\theta_{C}(t)\equiv\mbox{$\mathbb{E}$}\left[\frac{\partial}{\partial
t}\mu(t,\bm{S})\Big{|}T=t\right]$ in (4), which only relies on the consistent
estimation of $\frac{\partial}{\partial t}\mu(t,\bm{s})$ in the high density
region of $p(\bm{s}|t)$. On the other hand, our proposed estimator of $m(t)$
generalizes the regression adjustment estimator (6) but remains consistent and
numerically stable even when the conditional density $p(t|\bm{s})$ is zero for
some values of $\bm{s}\in\mathcal{S}$ thanks to its integral formulation. We
also provide a fast algorithm for computing our proposed estimator of $m(t)$
in practice and delineate the bootstrap inference procedures for both
estimators of $\theta(t)$ and $m(t)$.

###  3.1 Proposed Integral Estimator of $m(t)$

Given the observed data $\left\\{(Y_{i},T_{i},\bm{S}_{i})\right\\}_{i=1}^{n}$,
our proposal for removing the reliance of the positivity condition (A0) from
(6) is based on three critical insights.

  * •

Insight 1: Consistent estimation of $\mu(t,\bm{s})$ and
$\frac{\partial}{\partial t}\mu(t,\bm{s})$ at each $(T_{i},\bm{S}_{i})$. Given
that the observed data $(T_{i},\bm{S}_{i}),i=1,...,n$ generally appear in a
high density region of $p(t,\bm{s})$, the conditional mean outcome function
$\mu(t,\bm{s})$ can be well-estimated at each observation
$(T_{i},\bm{S}_{i})$. As a result, one can expect that the partial derivative
$\frac{\partial}{\partial t}\mu(t,\bm{s})$ can be consistently estimated at
each observation $(T_{i},\bm{S}_{i})$ as well.

  * •

Insight 2: Consistent estimation of $\theta(t)$ from a localized form
$\theta_{C}(t)$. The equations (4) in Assumption A2 lead to two possible
approaches to estimating the derivative effect
$\theta(t)=\frac{d}{dt}\mathbb{E}\left[\mu(t,\bm{S})\right]$. One approach is
via $\theta_{M}(t)=\mbox{$\mathbb{E}$}\left[\frac{\partial}{\partial
t}\mu(t,\bm{S})\right]$, while the alternative is based on
$\theta_{C}(t)=\mbox{$\mathbb{E}$}\left[\frac{\partial}{\partial
t}\mu(t,\bm{S})\Big{|}T=t\right]$. Estimating $\theta(t)$ by $\theta_{C}(t)$
is preferable because $\theta_{M}(t)$ relies on a consistent estimator of
$\frac{\partial}{\partial t}\mu(t,\bm{s})$ at each pair $(t,\bm{S}_{i})$,
which is not possible when the positivity condition (A0) fails to hold at some
$(t,\bm{S}_{i})$ for $i=1,...,n$. In contrast, the expression $\theta_{C}(t)$
only requires the estimator of $\frac{\partial}{\partial t}\mu(t,\bm{s})$ to
be accurate at the covariate vector $\bm{s}$ with a high conditional density
value $p(\bm{s}|t)$.

  * •

Insight 3: Integral relation between $\theta(t)$ and $m(t)$. For any
$t\in\mathcal{T}$, the fundamental theorem of calculus reveals that

$m(t)=m(T)+\int_{\widetilde{t}=T}^{\widetilde{t}=t}m^{\prime}(\widetilde{t})\,d\widetilde{t}=m(T)+\int_{\widetilde{t}=T}^{\widetilde{t}=t}\theta(\widetilde{t})\,d\widetilde{t}.$

Under Assumption A2, we take the expectation on both sides of the above
equality to obtain that

$\displaystyle\begin{split}m(t)=\mbox{$\mathbb{E}$}\left[m(T)+\int_{\widetilde{t}=T}^{\widetilde{t}=t}\theta(\widetilde{t})\,d\widetilde{t}\right]&=\mbox{$\mathbb{E}$}\left[\mu(T,\bm{S})\right]+\mbox{$\mathbb{E}$}\left[\int_{\widetilde{t}=T}^{\widetilde{t}=t}\theta_{C}(\widetilde{t})\,d\widetilde{t}\right]\\\
&=\mbox{$\mathbb{E}$}(Y)+\mbox{$\mathbb{E}$}\left\\{\int_{\widetilde{t}=T}^{\widetilde{t}=t}\mbox{$\mathbb{E}$}\left[\frac{\partial}{\partial
t}\mu(\widetilde{t},\bm{S})\Big{|}T=\widetilde{t}\right]\,d\widetilde{t}\right\\}.\end{split}$
(7)

This expression suggests that as long as we have consistent estimators of
$\mu(t,\bm{s})$ at $(T,\bm{S})$ and $\theta(t)$ near $T$, we can then use the
integration to extrapolate the estimation to $\mu(t,\bm{S})$ and $m(t)$ for
any $t\in\mathcal{T}$ even when $p(t,\bm{S})=0$.

According to the integral relation (7), we propose an _integral estimator_ of
the dose-response curve $m(t)$ as:

$\widehat{m}_{\theta}(t)=\frac{1}{n}\sum_{i=1}^{n}\left[Y_{i}+\int_{\widetilde{t}=T_{i}}^{\widetilde{t}=t}\widehat{\theta}_{C}(\widetilde{t})\,d\widetilde{t}\right],$
(8)

where $\widehat{\theta}_{C}(t)$ is a consistent estimator of
$\theta_{C}(t)=\mbox{$\mathbb{E}$}\left[\frac{\partial}{\partial
t}\mu(t,\bm{S})\Big{|}T=t\right]=\int\frac{\partial}{\partial
t}\mu(t,\bm{s})\,d\mbox{$\mathrm{P}$}(\bm{s}|t)$. To construct an estimator of
$\theta_{C}(t)$, we need to estimate two nuisance functions: (i) the partial
derivative $\frac{\partial}{\partial t}\mu(t,\bm{s})$ of the conditional mean
outcome function $\mu(t,\bm{s})$ and (ii) the conditional cumulative
distribution function (CDF) $P(\bm{s}|t)$. In this paper, we leverage the
following two kernel smoothing methods for estimating these two nuisance
functions and leave other possibilities as a future direction.

####  3.1.1 Local Polynomial Regression Estimator of $\frac{\partial}{\partial
t}\mu(t,\mathbf{s})$

We consider estimating $\frac{\partial}{\partial t}\mu(t,\bm{s})$ by the local
polynomial regression (Fan and Gijbels, 1996) because of its robustness around
the boundary of support (also known as the automatic kernel carpentry in
Hastie and Loader 1993).

Let $K_{T}:\mathbb{R}\to[0,\infty),K_{S}:\mathbb{R}^{d}\to[0,\infty)$ be two
symmetric kernel functions and $h,b>0$ be their corresponding smoothing
bandwidth parameters. Some commonly used univariate kernel functions include
the Epanechnikov kernel
$K(u)=\frac{3}{4}\left(1-u^{2}\right)\cdot\mathbbm{1}_{\\{|u|\leq 1\\}}$ and
Gaussian kernel $K(u)=\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{u^{2}}{2}\right)$.
For the multivariate kernel function, one often resorts to the product kernel
technique as $K_{S}(\bm{u})=\prod_{i=1}^{d}K(u_{i})$ for
$\bm{u}\in\mathbb{R}^{d}$. To estimate $\frac{\partial}{\partial
t}\mu(t,\bm{s})$ from the observed data
$\\{(Y_{i},T_{i},\bm{S}_{i})\\}_{i=1}^{n}$, we fit a partial local polynomial
regression of order $q$ ($q\geq 1$) with monomials of
$\\{(T_{i}-t)\\}_{i=1}^{n}$ as the polynomial basis in treatment variable $T$
and the local linear function in covariate vector $\bm{S}$ (Ruppert and Wand,
1994). Specifically, we let
$\bm{X}_{i}(t,\bm{s})=\left(1,(T_{i}-t),...,(T_{i}-t)^{q},(S_{i,1}-s_{1}),...,(S_{i,d}-s_{d})\right)^{T}\in\mathbb{R}^{q+1+d}$
for $i=1,...,n$ and consider

$\displaystyle\left(\widehat{\bm{\beta}}(t,\bm{s}),\widehat{\bm{\alpha}}(t,\bm{s})\right)^{T}$
(9)
$\displaystyle=\operatorname*{arg\,min}_{(\bm{\beta},\bm{\alpha})^{T}\in\mathbb{R}^{q+1}\times\mathbb{R}^{d}}\sum_{i=1}^{n}\left[Y_{i}-\sum_{j=0}^{q}\beta_{j}(T_{i}-t)^{q}-\sum_{\ell=1}^{d}\alpha_{\ell}(S_{i,\ell}-s_{\ell})\right]^{2}K_{T}\left(\frac{T_{i}-t}{h}\right)K_{S}\left(\frac{\bm{S}_{i}-\bm{s}}{b}\right)$
$\displaystyle=\operatorname*{arg\,min}_{(\bm{\beta},\bm{\alpha})^{T}\in\mathbb{R}^{q+1}\times\mathbb{R}^{d}}\sum_{i=1}^{n}\left(Y_{i}-\bm{X}_{i}(t,\bm{s})^{T}\begin{pmatrix}\bm{\beta}\\\
\bm{\alpha}\end{pmatrix}\right)^{2}K_{T}\left(\frac{T_{i}-t}{h}\right)K_{S}\left(\frac{\bm{S}_{i}-\bm{s}}{b}\right).$

For simplicity, we use the same bandwidth parameter $b$ for each coordinate in
$K_{S}$ here. One can generalize the above method and related theoretical
results to a general bandwidth matrix in $K_{S}$ with little effort. Let
$\bm{X}(t,\bm{s})\in\mathbb{R}^{n\times(q+1+d)}$ be a matrix with the $j$-th
row as $\bm{X}_{j}(t,\bm{s})^{T}$ and
$\bm{Y}=(Y_{1},...,Y_{n})^{T}\in\mathbb{R}^{n}$. We also define a diagonal
weight matrix as:

$\bm{W}(t,\bm{s})=\mathtt{Diag}\left(K_{T}\left(\frac{T_{1}-t}{h}\right)K_{S}\left(\frac{\bm{S}_{1}-\bm{s}}{b}\right),...,K_{T}\left(\frac{T_{n}-t}{h}\right)K_{S}\left(\frac{\bm{S}_{n}-\bm{s}}{b}\right)\right)\in\mathbb{R}^{n\times
n}.$

Then, (9) has a closed-form solution from a weighted least square problem as:

$\left(\widehat{\bm{\beta}}(t,\bm{s}),\widehat{\bm{\alpha}}(t,\bm{s})\right)^{T}=\left[\bm{X}^{T}(t,\bm{s})\bm{W}(t,\bm{s})\bm{X}(t,\bm{s})\right]^{-1}\bm{X}(t,\bm{s})^{T}\bm{W}(t,\bm{s})\bm{Y}.$
(10)

Finally, the second component $\widehat{\beta}_{2}(t,\bm{s})$ of the fitted
coefficient $\widehat{\bm{\beta}}(t,\bm{s})$ provides a natural estimator of
$\beta_{2}(t,\bm{s})\equiv\frac{\partial}{\partial t}\mu(t,\bm{s})$. In
practice, we recommend to choose $q$ to be an even number when estimating the
first-order (partial) derivative $\frac{\partial}{\partial t}\mu(t,\bm{s})$
via (9), because there is an increment to the asymptotic variance of
$\widehat{\beta}_{2}(t,\bm{s})$ when $q$ changes from an even number to the
consecutive odd number. Additionally, fitting (partial) local polynomial
regressions of higher orders often give rise to a possible reduction of bias
but also a substantial increase of the variability; see Chapter 3.3 in Fan and
Gijbels (1996). Therefore, we mainly focus on the (partial) local quadratic
regression $q=2$ when constructing our derivative estimator
$\widehat{\beta}_{2}(t,\bm{s})$ in the subsequent analysis.

####  3.1.2 Kernel-Based Conditional CDF Estimator of
$\mbox{$\mathrm{P}$}(\mathbf{s}|t)$

We consider estimating $\mbox{$\mathrm{P}$}(\bm{s}|t)$ through a Nadaraya-
Watson conditional CDF estimator (Hall et al., 1999) defined as:

$\widehat{P}_{\hslash}(\bm{s}|t)=\frac{\sum_{i=1}^{n}\mathbbm{1}_{\\{\bm{S}_{i}\leq\bm{s}\\}}\cdot\bar{K}_{T}\left(\frac{T_{i}-t}{\hslash}\right)}{\sum_{j=1}^{n}\bar{K}_{T}\left(\frac{T_{j}-t}{\hslash}\right)},$
(11)

where $\bar{K}_{T}:\mathbb{R}\to[0,\infty)$ is again a kernel function and
$\hslash>0$ is the smoothing bandwidth parameter that needs not be the same as
the bandwidth parameter $h>0$ for estimating
$\beta_{2}(t,\bm{s})=\frac{\partial}{\partial t}\mu(t,\bm{s})$ by the local
polynomial regression. Practically, there are several strategies for choosing
the bandwidth $\hslash>0$ in (11) as described by Bashtannyk and Hyndman
(2001); Holmes et al. (2012).

Combining the partial derivative estimator $\widehat{\beta}_{2}(t,\bm{s})$ in
(10) with the conditional CDF estimator $\widehat{P}_{\hslash}(\bm{s}|t)$ in
(11), we deduce the final localized estimator of
$\theta_{C}(t)=\int\frac{\partial}{\partial
t}\mu(t,\bm{s})\,d\mbox{$\mathrm{P}$}(\bm{s}|t)$ as:

$\displaystyle\widehat{\theta}_{C}(t)$
$\displaystyle=\int\widehat{\beta}_{2}(t,\bm{s})\,d\widehat{P}_{\hslash}(\bm{s}|t)=\frac{\sum_{i=1}^{n}\widehat{\beta}_{2}(t,\bm{S}_{i})\cdot\bar{K}_{T}\left(\frac{T_{i}-t}{\hslash}\right)}{\sum_{j=1}^{n}\bar{K}_{T}\left(\frac{T_{j}-t}{\hslash}\right)}.$
(12)

In essence, $\widehat{\theta}_{C}(t)$ is a regression adjustment estimator
with two nuisance functions: $\beta_{2}(t,\bm{s})=\frac{\partial}{\partial
t}\mu(t,\bm{s})$ and $P(\bm{s}|t)$. Here, we have showcased how to use the
(partial) local polynomial regression (9) to estimate $\beta_{2}(t,\bm{s})$
and Nadaraya-Watson conditional CDF estimator (11) to estimate $P(\bm{s}|t)$.

######  Remark 1 (Regression adjustment estimator of $\theta(t)$).

Under the additive confounding model (3), one can directly estimate
$\theta(t)$ via $\widehat{\beta}_{2}(t,\bm{s})$ or the regression adjustment
estimator
$\widehat{\theta}_{RA}(t)=\frac{1}{n}\sum_{i=1}^{n}\widehat{\beta}_{2}(t,\bm{S}_{i})$,
because by Proposition 1,

$\theta(t)=\frac{d}{dt}\mathbb{E}\left[\mu(t,\bm{S})\right]=\frac{d}{dt}\mathbb{E}\left[\bar{m}(t)+\eta(\bm{S})\right]=\bar{m}^{\prime}(t)=\frac{\partial}{\partial
t}\mu(t,\bm{s}).$

The estimator $\widehat{\beta}_{2}(t,\bm{s})$ is suboptimal, because it only
uses the value from a single location. Furthermore, $\widehat{\theta}_{RA}(t)$
will not be a stable and consistent estimator of $\theta(t)$ when the
positivity condition fails to hold at those
$(t,\bm{S}_{i})\in\mathcal{T}\times\mathcal{S}$; see Fig 1 for an
illustration. In contrast, what we propose in (12) relies on the conditional
CDF estimator and remains valid and consistent even without the positivity
condition (A0).

######  Remark 2 (Linear smoother).

It is worth noting that our integral estimator (8) under the kernel-based
estimator (12) is that it is a linear smoother. Let
$\bm{e}_{2}=(0,1,0,...,0)^{T}\in\mathbb{R}^{q+1+d}$. Then, with (10) and
$\widehat{\beta}_{2}(t,\bm{s})=\bm{e}_{2}^{T}\widehat{\bm{\beta}}(t,\bm{s})$,
we know that $\widehat{\theta}(t)$ in (12) can also be written as:

$\widehat{\theta}_{C}(t)=\left\\{\int\bm{e}_{2}^{T}\left[\bm{X}^{T}(t,\bm{s})\bm{W}(t,\bm{s})\bm{X}(t,\bm{s})\right]^{-1}\bm{X}(t,\bm{s})^{T}\bm{W}(t,\bm{s})\,d\widehat{P}_{\hslash}(\bm{s}|t)\right\\}\bm{Y}.$

It implies that

$\displaystyle\widehat{m}_{\theta}(t)$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}Y_{i}+\left\\{\frac{1}{n}\sum_{i=1}^{n}\int_{\widetilde{t}=T_{i}}^{\widetilde{t}=t}\int\bm{e}_{2}^{T}\left[\bm{X}^{T}(t,\bm{s})\bm{W}(t,\bm{s})\bm{X}(t,\bm{s})\right]^{-1}\bm{X}(t,\bm{s})^{T}\bm{W}(t,\bm{s})\,d\widehat{P}_{\hslash}(\bm{s}|t)\,d\widetilde{t}\right\\}\bm{Y}$
$\displaystyle\equiv\sum_{i=1}^{n}l_{i}(t)Y_{i}.$

As a consequence, we can also utilize the theory of linear smoothers to derive
its effective degrees of freedom and fine-tune the smoothing bandwidth
parameters (Buja et al., 1989; Wasserman, 2006).

###  3.2 A Fast Computing Algorithm for the Proposed Integral Estimator

Our proposed estimator
$\widehat{m}_{\theta}(t)=\frac{1}{n}\sum_{i=1}^{n}\left[Y_{i}+\int_{T_{i}}^{t}\widehat{\theta}_{C}(\widetilde{t})\,d\widetilde{t}\right]$
of the dose-response curve $m(t)$ involves an integral that may be
analytically difficult to compute in practice. Here, we propose a fast
computing algorithm that can numerically approximate $\widehat{m}_{\theta}(t)$
with an error of the order at most $O_{P}\left(\frac{1}{n}\right)$. The key
idea is to approximate the integral via a Riemann sum, evaluate
$\widehat{m}_{\theta}(t)$ only at the data sample $T_{1},...,T_{n}$, and then
use the linear interpolation to obtain the value at any arbitrary
$t\in\mathcal{T}$.

Let $T_{(1)}\leq\cdots\leq T_{(n)}$ be the order statistics of
$T_{1},...,T_{n}$ and $\Delta_{j}=T_{(j+1)}-T_{(j)}$ for $j=1,...,n-1$ be
their consecutive differences. The integral estimator
$\widehat{m}_{\theta}(t)$ in (8) can be rewritten as:

$\widehat{m}_{\theta}(t)=\frac{1}{n}\sum_{i=1}^{n}Y_{i}+\frac{1}{n}\sum_{i=1}^{n}\int_{\widetilde{t}=T_{(i)}}^{\widetilde{t}=t}\widehat{\theta}_{C}(\widetilde{t})\,d\widetilde{t}.$

To compute $\widehat{m}_{\theta}(T_{(j)})$ evaluated at the $j$-th order
statistic, we consider approximating the second integral term
$\frac{1}{n}\sum_{i=1}^{n}\int_{T_{(i)}}^{T_{(j)}}\widehat{\theta}_{C}(\widetilde{t})\,d\widetilde{t}$
as follows.

  * •

When $i<j$, we have the following Riemann sum approximation as:

$\displaystyle\int_{\widetilde{t}=T_{(i)}}^{\widetilde{t}=T_{(j)}}\widehat{\theta}_{C}(\widetilde{t})d\widetilde{t}$
$\displaystyle\approx\sum_{\ell=i}^{j-1}\widehat{\theta}_{C}(T_{(\ell)})\Delta_{\ell}.$

  * •

When $i>j$, we use another Riemann sum approximation as:

$\displaystyle\int_{\widetilde{t}=T_{(i)}}^{\widetilde{t}=T_{(j)}}\widehat{\theta}_{C}(\widetilde{t})d\widetilde{t}$
$\displaystyle\approx-\sum_{\ell=j}^{i-1}\widehat{\theta}_{C}(T_{(\ell+1)})\Delta_{\ell}.$

Plugging the above results into
$\frac{1}{n}\sum_{i=1}^{n}\int_{T_{(i)}}^{T_{(j)}}\widehat{\theta}_{C}(\widetilde{t})\,d\widetilde{t}$,
we obtain that for $1<j<n$,

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\int_{\widetilde{t}=T_{(i)}}^{\widetilde{t}=T_{(j)}}\widehat{\theta}_{C}(\widetilde{t})d\widetilde{t}$
$\displaystyle\approx\frac{1}{n}\left[\sum_{i=1}^{j-1}\sum_{\ell=i}^{j-1}\widehat{\theta}_{C}(T_{(\ell)})\Delta_{\ell}-\sum_{i=j}^{n}\sum_{\ell=j}^{i-1}\widehat{\theta}_{C}(T_{(\ell+1)})\Delta_{\ell}\right]$
$\displaystyle\stackrel{{\scriptstyle\text{(i)}}}{{=}}\frac{1}{n}\left[\sum_{\ell=1}^{j-1}\ell\cdot\widehat{\theta}_{C}(T_{(\ell)})\Delta_{\ell}-\sum_{\ell=j}^{n-1}(n-\ell)\cdot\widehat{\theta}_{C}(T_{(\ell+1)})\Delta_{\ell}\right]$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n-1}\Delta_{i}\left[i\cdot\widehat{\theta}_{C}(T_{(i)})\mathbbm{1}_{\\{i<j\\}}-(n-i)\cdot\widehat{\theta}_{C}(T_{(i+1)})\mathbbm{1}_{\\{i\geq
j\\}}\right],$

where the equality (i) follows from switching the orders of summations.
Similarly, when $j=1$ or $j=n$, we also have that

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\int_{\widetilde{t}=T_{(i)}}^{\widetilde{t}=T_{(j)}}\widehat{\theta}_{C}(\widetilde{t})d\widetilde{t}$
$\displaystyle\approx\begin{cases}-\frac{1}{n}\sum_{i=2}^{n}\sum_{\ell=1}^{i-1}\widehat{\theta}_{C}(T_{(\ell+1)})\Delta_{\ell}&\text{
when }j=1,\\\
\frac{1}{n}\sum_{i=1}^{n-1}\sum_{\ell=i}^{n-1}\widehat{\theta}_{C}(T_{(\ell)})\Delta_{\ell}&\text{
when }j=n,\end{cases}$
$\displaystyle=\begin{cases}-\frac{1}{n}\sum_{\ell=1}^{n-1}(n-\ell)\cdot\widehat{\theta}_{C}(T_{(\ell+1)})\Delta_{\ell}&\text{
when }j=1,\\\
\frac{1}{n}\sum_{\ell=1}^{n-1}\ell\cdot\widehat{\theta}_{C}(T_{(\ell)})\Delta_{\ell}&\text{
when }j=n,\end{cases}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n-1}\Delta_{i}\left[i\cdot\widehat{\theta}_{C}(T_{(i)})\mathbbm{1}_{\\{i<j\\}}-(n-i)\cdot\widehat{\theta}_{C}(T_{(i+1)})\mathbbm{1}_{\\{i\geq
j\\}}\right].$

Therefore, we propose the following approximation for
$\widehat{m}_{\theta}(T_{(j)})$ as:

$\widehat{m}_{\theta}(T_{(j)})\approx\frac{1}{n}\sum_{i=1}^{n}Y_{i}+\frac{1}{n}\sum_{i=1}^{n-1}\Delta_{i}\left[i\cdot\widehat{\theta}_{C}(T_{(i)})\mathbbm{1}_{\\{i<j\\}}-(n-i)\cdot\widehat{\theta}_{C}(T_{(i+1)})\mathbbm{1}_{\\{i\geq
j\\}}\right].$ (13)

Finally, to evaluate $\widehat{m}_{\theta}(t)$ at any arbitrary $t$, we
conduct the linear interpolation between $\widehat{m}_{\theta}(T_{(j)})$ and
$\widehat{m}_{\theta}(T_{(j+1)})$ on the interval
$t\in\left[T_{(j)},T_{(j+1)}\right]$. One biggest advantage of using the
approximation formula (13) is that we only need to compute the derivative
estimator $\widehat{\theta}_{C}(t)$ at the order statistics
$T_{(1)},...,T_{(n)}$. When $\widehat{m}_{\theta}(t)$ and its derivative
$\widehat{\theta}_{C}(t)$ are Lipschitz and the marginal density $p_{T}(t)$ is
uniformly bounded away from 0 on the region of interest, this approximation
formula has at most $O_{P}\left(\frac{1}{n}\right)$ error, which is
asymptotically negligible compared to the dominating estimation error of
$\widehat{m}_{\theta}(t)$; see Theorem 4 for details.

###  3.3 Bootstrap Inference

Since it is complicated to derive consistent estimators of the (asymptotic)
variances of our integral estimator (8) and localized derivative estimator
(12), we consider conducting inference on $m(t)$ and $\theta(t)$ through the
empirical bootstrap method (Efron, 1979) as follows. Other bootstrap methods,
including residual bootstrap (Freedman, 1981) and wild bootstrap (Wu, 1986),
also work under some modified conditions.

  1. 1.

Compute the integral estimator $\widehat{m}_{\theta}(t)$ and localized
derivative estimator $\widehat{\theta}_{C}(t)$ on the original data
$\\{(Y_{i},T_{i},\bm{S}_{i})\\}_{i=1}^{n}$.

  2. 2.

Generate $B$ bootstrap samples
$\left\\{\left(Y_{i}^{*(b)},T_{i}^{*(b)},\bm{S}_{i}^{*(b)}\right)\right\\}_{i=1}^{n},b=1,...,B$
by sampling with replacement from the original data and compute the integral
estimator $\widehat{m}_{\theta}^{*(b)}(t)$ and localized derivative estimator
$\widehat{\theta}_{C}^{*(b)}(t)$ on each bootstrapped sample for $b=1,...,B$.

  3. 3.

Let $\alpha\in(0,1)$ be a pre-specified significance level.

     * •

For a pointwise inference at $t_{0}\in\mathcal{T}$, we calculate the
$1-\alpha$ quantiles $\zeta_{1-\alpha}^{*}(t_{0})$ and
$\bar{\zeta}_{1-\alpha}^{*}(t_{0})$ of $\\{D_{1}(t_{0}),...,D_{B}(t_{0})\\}$
and $\\{\bar{D}_{1}(t_{0}),...,\bar{D}_{B}(t_{0})\\}$ respectively, where
$D_{b}(t_{0})=\left|\widehat{m}_{\theta}^{*(b)}(t_{0})-\widehat{m}_{\theta}(t_{0})\right|$
and
$\bar{D}_{b}(t_{0})=\left|\widehat{\theta}_{C}^{*(b)}(t_{0})-\widehat{\theta}_{C}(t_{0})\right|$
for $b=1,...,B$.

     * •

For an uniform inference on the entire dose-response curve $m(t)$ and its
derivative $\theta(t)$, we compute the $1-\alpha$ quantiles
$\xi_{1-\alpha}^{*}$ and $\bar{\xi}_{1-\alpha}^{*}$ of
$\\{D_{\sup,1},...,D_{\sup,B}\\}$ and
$\\{\bar{D}_{\sup,1},...,\bar{D}_{\sup,B}\\}$ respectively, where
$D_{\sup,b}=\sup_{t\in\mathcal{T}}\left|\widehat{m}_{\theta}^{*(b)}(t)-\widehat{m}_{\theta}(t)\right|$
and
$\bar{D}_{\sup,b}=\sup_{t\in\mathcal{T}}\left|\widehat{\theta}_{C}^{*(b)}(t)-\widehat{\theta}_{C}(t)\right|$
for $b=1,...,B$.

  4. 4.

Define the $1-\alpha$ confidence intervals for $m(t_{0})$ and $\theta(t_{0})$
as:

$\left[\widehat{m}_{\theta}(t_{0})-\zeta_{1-\alpha}^{*}(t_{0}),\,\widehat{m}_{\theta}(t_{0})+\zeta_{1-\alpha}^{*}(t_{0})\right]\;\text{
and
}\;\left[\widehat{\theta}_{C}(t_{0})-\bar{\zeta}_{1-\alpha}^{*}(t_{0}),\,\widehat{\theta}_{C}(t_{0})+\bar{\zeta}_{1-\alpha}^{*}(t_{0})\right]$

respectively, as well as the simultaneous $1-\alpha$ confidence bands as:

$\left[\widehat{m}_{\theta}(t)-\xi_{1-\alpha}^{*},\,\widehat{m}_{\theta}(t)+\xi_{1-\alpha}^{*}\right]\quad\text{
and
}\quad\left[\widehat{\theta}_{C}(t)-\bar{\xi}_{1-\alpha}^{*},\,\widehat{\theta}_{C}(t)+\bar{\xi}_{1-\alpha}^{*}\right]$

for every $t\in\mathcal{T}$.

In Section 4.3, we will establish the consistency of the above bootstrap
inference procedures.

##  4 Asymptotic Theory

In this section, we study the consistency results of our integral estimator
(8) and localized derivative estimator (12) proposed in Section 3.1 and the
validity of bootstrap inference described in Section 3.3.

###  4.1 Notations and Assumptions

We introduce the regularity conditions under the general confounding model (2)
for our subsequent theoretical analysis. Let
$\mathcal{E}\subset\mathcal{T}\times\mathcal{S}$ be the support of the joint
density $p(t,\bm{s})$, $\mathcal{E}^{\circ}$ be the interior of $\mathcal{E}$,
and $\partial\mathcal{E}$ be the boundary of $\mathcal{E}$.

######  Assumption A3 (Differentiability of the conditional mean outcome
function).

For any $(t,\bm{s})\in\mathcal{T}\times\mathcal{S}$, the conditional mean
outcome function $\mu(t,\bm{s})$ is at least $(q+1)$ times continuously
differentiable with respect to $t$ and at least four times continuously
differentiable with respect to $\bm{s}$, where $q$ is the order of (partial)
local polynomial regression in (9). Furthermore, $\mu(t,\bm{s})$ and all of
its partial derivatives are uniformly bounded on
$\mathcal{T}\times\mathcal{S}$.

######  Assumption A4 (Differentiability of the joint density).

The joint density $p(t,\bm{s})$ is bounded and at least twice continuously
differentiable with bounded partial derivatives up to the second order on
$\mathcal{E}^{\circ}$. All these partial derivatives of $p(t,\bm{s})$ are
continuous up to the boundary $\partial\mathcal{E}$. Furthermore,
$\mathcal{E}$ is compact and $p(t,\bm{s})$ is uniformly bounded away from 0 on
$\mathcal{E}$. Finally, the marginal density $p_{T}(t)$ of $T$ is non-
degenerate, _i.e._ , its support $\mathcal{T}$ has a nonempty interior.

Assumptions A3 and A4 are commonly assumed differentiability conditions in the
literature of local polynomial regression (Ruppert and Wand, 1994; Fan and
Gijbels, 1996), which can be slightly relaxed via the Hölder continuity
condition. It ensures that the bias term of $\widehat{\beta}_{2}(t,\bm{s})$
from the local polynomial regression (10) is at least of the standard order
$O(h^{q})+O(b^{2})$; see Lemma 2 below. Notice that the projection
$\mathrm{proj}_{T}(\mathcal{E})$ of the joint density support $\mathcal{E}$
onto the domain of $T$ coincides with the marginal support $\mathcal{T}$.
Hence, $\mathcal{T}$ will be compact as well under Assumption A4.

To control the boundary effects of the local polynomial regression, we impose
the following conditions on $p(t,\bm{s})$, $\mu(t,\bm{s})$, and the geometric
structure of $\mathcal{E}$ near the boundary $\partial\mathcal{E}$.

######  Assumption A5 (Boundary conditions).

  1. (a)

There exists some constants $r_{1},r_{2}\in(0,1)$ such that for any
$(t,\bm{s})\in\mathcal{E}$ and all $\delta\in(0,r_{1}]$, there is a point
$(t^{\prime},\bm{s}^{\prime})\in\mathcal{E}$ satisfying

$\mathcal{B}\left((t^{\prime},\bm{s}^{\prime}),\,r_{2}\delta\right)\subset\mathcal{B}\left((t,\bm{s}),\,\delta\right)\cap\mathcal{E},$

where
$\mathcal{B}((t,\bm{s}),\,r)=\left\\{(t_{1},\bm{s}_{1})\in\mathbb{R}^{d+1}:\left|\left|(t_{1}-t,\bm{s}_{1}-\bm{s})\right|\right|_{2}\leq
r\right\\}$ with $\left|\left|\cdot\right|\right|_{2}$ being the standard
Euclidean norm.

  2. (b)

For any $(t,\bm{s})\in\partial\mathcal{E}$, the boundary of $\mathcal{E}$, it
satisfies that $\frac{\partial}{\partial
t}p(t,\bm{s})=\frac{\partial}{\partial s_{j}}p(t,\bm{s})=0$ and
$\frac{\partial^{2}}{\partial s_{j}^{2}}\mu(t,\bm{s})=0$ for all $j=1,...,d$.

  3. (c)

For any $\delta>0$, the Lebesgue measure of the set
$\partial\mathcal{E}\oplus\delta$ satisfies
$\left|\partial\mathcal{E}\oplus\delta\right|\leq A_{1}\cdot\delta$ for some
absolute constant $A_{1}>0$, where
$\partial\mathcal{E}\oplus\delta=\left\\{\bm{z}\in\mathbb{R}^{d+1}:\inf_{\bm{x}\in\partial\mathcal{E}}\left|\left|\bm{z}-\bm{x}\right|\right|_{2}\leq\delta\right\\}$.

Assumption A5(a) is adopted from the boundary condition (Assumption X) in Fan
and Guerre (2015), whose primitive and stronger version also appeared as
Assumption (A4) in Ruppert and Wand (1994). This is a relatively mild
condition that holds when the boundary $\partial\mathcal{E}$ is smooth or only
contains non-smooth vertices from some regular structures, such as $d+1$
dimensional cubes and convex cones. Indeed, Assumption A5(a) is valid as long
as any ball centered at a point in $\mathcal{E}$ near a vertex of
$\partial\mathcal{E}$ has its radius shrunk linearly when approaching the
vertex. More examples and discussions about Assumption A5(b) can be found in
Ruppert and Wand (1994); Fan and Guerre (2015). The main purpose of imposing
this support condition is to ensure that there are enough observations near
the boundary points so that the rate of convergence for our local polynomial
estimator $\widehat{\beta}_{2}(t,\bm{s})$ remains the same order at the
boundary points as the interior points of $\mathcal{E}$. Assumption A5(b)
regularizes the slope of $p(t,\bm{s})$ and the curvature of $\mu(t,\bm{s})$ at
the boundary point $(t,\bm{s})\in\partial\mathcal{E}$. Specifically, at the
boundary $\partial\mathcal{E}$, the joint density $p(t,\bm{s})$ needs to be
flat, while $\mu(t,\bm{s})$ should embrace zero curvatures. Similar to
Assumption A4, the partial derivatives $\frac{\partial}{\partial
t}p(t,\bm{s}),\frac{\partial}{\partial s_{j}}p(t,\bm{s})$ are defined by
computing them at a nearby interior point $(t^{\prime},\bm{s}^{\prime})$ and
taking the limit
$(t^{\prime},\bm{s}^{\prime})\to(t,\bm{s})\in\partial\mathcal{E}$. This
condition is another key requirement for the bias term from
$\widehat{\beta}_{2}(t,\bm{s})$ in the local polynomial regression (10) to
remain the same rate of convergence at the boundary points as the interior
points of $\mathcal{E}$. Finally, Assumption A5(c) also regularizes the
boundary $\partial\mathcal{E}$ so that it will not have any fractal or other
peculiar structures leading to an infinite perimeter of $\partial\mathcal{E}$.

To establish the (uniform) consistency of our local polynomial regression
estimator (10) and Nadaraya-Watson conditional CDF estimator (11), we rely on
the following regularity conditions on the kernel functions.

######  Assumption A6 (Regular kernel and VC-type conditions).

  1. (a)

The functions $K_{T}:\mathbb{R}\to[0,\infty)$ and
$K_{S}:\mathbb{R}^{d}\to[0,\infty)$ are compactly supported and Lispchitz
continuous kernels such that
$\int_{\mathbb{R}}K_{T}(t)\,dt=\int_{\mathbb{R}^{d}}K_{S}(\bm{s})\,d\bm{s}=1$,
$K_{T}(t)=K_{T}(-t)$, and $K_{S}$ is radially symmetric with $\int\bm{s}\cdot
K_{S}(\bm{s})d\bm{s}=\bm{0}$. In addition, for all $j=1,2,...,$ and
$\ell=1,...,d$, it holds that

$\displaystyle\kappa_{j}^{(T)}:=\int_{\mathbb{R}}u^{j}K_{T}(u)\,du<\infty,\quad\nu_{j}^{(T)}:=\int_{\mathbb{R}}u^{j}K_{T}^{2}(u)\,du<\infty,$
$\displaystyle\kappa_{j,\ell}^{(S)}:=\int_{\mathbb{R}^{d}}u_{\ell}^{j}K_{S}(\bm{u})\,d\bm{u}<\infty,\quad\text{
and
}\quad\nu_{j,k}^{(S)}:=\int_{\mathbb{R}^{d}}u_{\ell}^{j}K_{S}^{2}(\bm{u})\,d\bm{u}<\infty.$

Finally, both $K_{T}$ and $K_{S}$ are second-order kernels, _i.e._ ,
$\kappa_{2}^{(T)}>0$ and $\kappa_{2,\ell}^{(S)}>0$ for all $\ell=1,...,d$.

  2. (b)

Let
$\mathcal{K}_{q,d}=\Big{\\{}(y,\bm{z})\mapsto\left(\frac{y-t}{h}\right)^{\ell}\left(\frac{z_{i}-s_{i}}{b}\right)^{k_{1}}\left(\frac{z_{j}-s_{j}}{b}\right)^{k_{2}}K_{T}\left(\frac{y-t}{h}\right)K_{S}\left(\frac{\bm{z}-\bm{s}}{b}\right):(t,\bm{s})\in\mathcal{T}\times\mathcal{S};i,j=1,...,d;\ell=0,...,2q;k_{1},k_{2}=0,1;h,b>0\Big{\\}}$.
It holds that $\mathcal{K}_{q,d}$ is a bounded VC-type class of measurable
functions on $\mathbb{R}^{d+1}$.

  3. (c)

The function $\bar{K}_{T}:\mathbb{R}\to[0,\infty)$ is a second-order,
Lipschitz continuous, and symmetric kernel with a compact support, _i.e._ ,
$\int_{\mathbb{R}}\bar{K}_{T}(t)\,dt=1$, $\bar{K}_{T}(t)=\bar{K}_{T}(-t)$, and
$\int_{\mathbb{R}}t^{2}\bar{K}_{T}(t)\,dt\in(0,\infty)$.

  4. (d)

Let
$\bar{\mathcal{K}}=\left\\{y\mapsto\bar{K}_{T}\left(\frac{y-t}{\hslash}\right):t\in\mathcal{T},\hslash>0\right\\}$.
It holds that $\bar{\mathcal{K}}$ is a bounded VC-type class of measurable
functions on $\mathbb{R}$.

Assumption A6(a,c) are indeed not the regularity conditions but rather
properties of those commonly used kernel functions. Assumption A6(b,d) are
critical conditions for the uniform consistency of kernel-based function
estimators (Giné and Guillou, 2002; Einmahl and Mason, 2005) and can be
satisfied by a wide range of kernel functions, including Gaussian and
Epanechnikov kernels. For example, it is satisfied when the kernel is a
composite function between a polynomial in $d\geq 1$ variables and a real-
valued function of bounded variation; see Lemma 22 in Nolan and Pollard
(1987). Given that the kernel functions are bounded, we can always take
constant envelope functions for $\mathcal{K}_{q,d}$ and $\bar{\mathcal{K}}$ in
Assumption A6(b,d).

Finally, we introduce some notations that will be used in proofs of the
consistency results of $\widehat{\beta}_{2}(t,\bm{s})$ from the local
polynomial regression (10) (Lemmas 2 and 3) and the asymptotic linearity
results of our integral estimator and localized derivative estimator (Lemma
5). Under the notations in Assumption A6(a), we define a matrix

$\bm{M}_{q}=\begin{pmatrix}\left(\kappa_{i+j-2}^{(T)}\right)_{1\leq i,j\leq
q+1}&\bm{0}\\\
\bm{0}&\left(\kappa_{2,i-q-1}^{(S)}\mathbbm{1}_{\\{i=j\\}}\right)_{q+1<i,j\leq
q+1+d}\end{pmatrix}\in\mathbb{R}^{(q+1+d)\times(q+1+d)}.$ (14)

Notice that $\bm{M}_{q}$ only depends on the kernel functions $K_{T},K_{S}$ in
the local polynomial regression. For any
$(t,\bm{s})\in\mathcal{T}\times\mathcal{S}$, we also define the functions
$\bm{\Psi}_{t,\bm{s}},\bm{\psi}_{t,\bm{s}}:\mathbb{R}\times\mathbb{R}\times\mathbb{R}^{d}\to\mathbb{R}^{q+1+d}$
as:

$\displaystyle\begin{split}\bm{\Psi}_{t,\bm{s}}(y,z,\bm{v})&=\begin{bmatrix}\left(y\cdot\left(\frac{z-t}{h}\right)^{j-1}K_{T}\left(\frac{z-t}{h}\right)K_{S}\left(\frac{\bm{v}-\bm{s}}{b}\right)\right)_{1\leq
j\leq q+1}\\\
\left(y\cdot\left(\frac{v_{j-q-1}-s_{j-q-1}}{b}\right)K_{T}\left(\frac{z-t}{h}\right)K_{S}\left(\frac{\bm{v}-\bm{s}}{b}\right)\right)_{q+1<j\leq
q+1+d}\end{bmatrix}\\\
&=y\cdot\begin{bmatrix}\left(\left(\frac{z-t}{h}\right)^{j-1}K_{T}\left(\frac{z-t}{h}\right)K_{S}\left(\frac{\bm{v}-\bm{s}}{b}\right)\right)_{1\leq
j\leq q+1}\\\
\left(\left(\frac{v_{j-q-1}-s_{j-q-1}}{b}\right)K_{T}\left(\frac{z-t}{h}\right)K_{S}\left(\frac{\bm{v}-\bm{s}}{b}\right)\right)_{q+1<j\leq
q+1+d}\end{bmatrix}\equiv y\cdot\bm{\psi}_{t,\bm{s}}(z,\bm{v}).\end{split}$
(15)

###  4.2 Consistency of the Integral Estimator

Before establishing the consistency of our proposed integral estimator (8) and
localized derivative estimator (12), we first discuss the pointwise and
uniform consistency results of $\widehat{\beta}_{2}(t,\bm{s})$ from the local
polynomial regression (10) as building blocks.

######  Lemma 2 (Pointwise convergence of $\widehat{\beta}_{2}(t,\bm{s})$).

Suppose that Assumptions A3, A4, A5, and A6(a) hold under the confounding
model (2). Let $\widehat{\beta}_{2}(t,\bm{s})$ be the second element of
$\widehat{\bm{\beta}}(t,\bm{s})\in\mathbb{R}^{q+1}$ defined by (10) and
$\beta_{2}(t,\bm{s})=\frac{\partial}{\partial t}\mu(t,\bm{s})$. Then, for any
fixed $(t,\bm{s})\in\mathcal{E}$ and integer $q>0$, as
$h,b,\frac{\max\\{h,b\\}^{4}}{h}\to 0$, and $nh^{3}b^{d}\to\infty$, we have
that

$\displaystyle\widehat{\beta}_{2}(t,\bm{s})-\beta_{2}(t,\bm{s})$
$\displaystyle=\begin{cases}O\left(h^{q+1}+b^{2}+\frac{b^{4}}{h}\right)+O_{P}\left(\sqrt{\frac{1}{nh^{3}b^{d}}}\right)&\text{
if }q\text{ is odd and }(t,\bm{s})\in\mathcal{E}^{\circ},\\\
O\left(h^{q}+b^{2}+\frac{b^{4}}{h}\right)+O_{P}\left(\sqrt{\frac{1}{nh^{3}b^{d}}}\right)&\text{
if }q\text{ is even and }(t,\bm{s})\in\mathcal{E}^{\circ},\\\
O\left(h^{q}+\frac{\max\\{h,b\\}^{4}}{h}\right)+O_{P}\left(\sqrt{\frac{1}{nh^{3}b^{d}}}\right)&\text{
if }q\text{ is an integer and }(t,\bm{s})\in\partial\mathcal{E}.\end{cases}$

The complete statement of the asymptotic expressions for the conditional
variances and bias of $\widehat{\beta}_{2}(t,\bm{s})$ and the proof of Lemma 2
are given in Section B.2. Notice that our (conditional) rates of convergence
of $\widehat{\beta}_{2}(t,\bm{s})$ in Lemma 2 align with the standard results
in the literature of local polynomial regression (Ruppert and Wand, 1994; Fan
and Gijbels, 1996; Lu, 1996). The extra bias rate
$O\left(\frac{b^{4}}{h}\right)$ or $O\left(\frac{\max\\{h,b\\}^{4}}{h}\right)$
comes from a higher order term in the Tyler’s expansion and will be
asymptotically negligible when $h\asymp b$ as $n\to\infty$; see (26) and (28)
in Section B.2 for an example. More importantly, when we utilize the (partial)
local quadratic regression (_i.e._ , $q=2$) to obtain
$\widehat{\beta}_{2}(t,\bm{s})$, the order of $h\asymp b$ that optimally
trades off the bias and variance in Lemma 2 will be
$O\left(n^{-\frac{1}{d+7}}\right)$, which also leads to the optimal rate for
derivative estimation established in Stone (1980, 1982). Conventional
bandwidth selection methods generally will not yield $h,b$ with this optimal
order. However, the validity of bootstrap inference in Section 4.3 requires
bandwidths $h,b$ to be of a smaller order $O\left(n^{-\frac{1}{d+5}}\right)$
so that the asymptotic bias is negligible compared to the stochastic variation
(Bjerve et al., 1985; Hall, 1992); see Lemma 5 below. Since
$O\left(n^{-\frac{1}{d+5}}\right)$ is an optimal bandwidth order for
estimating the regression function $\mu(t,\bm{s})$, we can then apply any off-
the-shelf bandwidth selection method of multivariate local polynomial
regression (Wand and Jones, 1994; Yang and Tschernig, 1999; Li and Racine,
2004) to obtain our bandwidth parameters $h\asymp b$; see also Remark 4 below.

We also strengthen the pointwise rate of convergence to the uniform one as
follows.

######  Lemma 3 (Uniform convergence of $\widehat{\beta}_{2}(t,\bm{s})$).

Let $q>0$. Suppose that Assumptions A3, A4, A5, and A6(a,b) hold under the
confounding model (2). Let $\widehat{\beta}_{2}(t,\bm{s})$ be the second
element of $\widehat{\bm{\beta}}(t,\bm{s})\in\mathbb{R}^{q+1}$ and
$\beta_{2}(t,\bm{s})=\frac{\partial}{\partial t}\mu(t,\bm{s})$. Then, as
$h,b,\frac{\max\\{h,b\\}^{4}}{h}\to 0$ and
$\frac{nh^{3}b^{d}}{|\log(hb^{d})|},\frac{|\log(hb^{d})|}{\log\log
n}\to\infty$, we have that

$\sup_{(t,\bm{s})\in\mathcal{E}}\left|\widehat{\beta}_{2}(t,\bm{s})-\beta_{2}(t,\bm{s})\right|=O\left(h^{q}+b^{2}+\frac{\max\\{b,h\\}^{4}}{h}\right)+O_{P}\left(\sqrt{\frac{|\log(hb^{d})|}{nh^{3}b^{d}}}\right).$

The proof of Lemma 3 can be found in Section B.3. At a high level, the uniform
rate of convergence for the bias term
$\mathbb{E}\left[\widehat{\beta}_{2}(t,\bm{s})\right]-\beta_{2}(t,\bm{s})$
remains the same as in Lemma 2 under our regularity conditions. To handle the
stochastic variation term
$\widehat{\beta}_{2}(t,\bm{s})-\mathbb{E}\left[\widehat{\beta}_{2}(t,\bm{s})\right]$,
we approximate it (up to a scaled factor $\sqrt{nh^{3}b^{d}}$) by an empirical
process so that the upper bound for its uniform rate of convergence follows
from the results in Einmahl and Mason (2005). This uniform rate of convergence
in Lemma 3 is not only useful for deriving the asymptotic behaviors of our
integral estimator (8) and localized derivative estimator (12) in Theorem 4
below but can also facilitate our study of the bootstrap consistency in
Section 4.3.

With Lemma 3, we now present the uniform consistency results for our integral
estimator (8) and localized derivative estimator (12).

######  Theorem 4 (Convergence of $\widehat{\theta}_{C}(t)$ and
$\widehat{m}_{\theta}(t)$).

Let $q>0$ and $\mathcal{T}^{\prime}\subset\mathcal{T}$ be a compact set so
that $p_{T}(t)$ is uniformly bounded away from 0 within
$\mathcal{T}^{\prime}$. Suppose that Assumptions A1, A2, A3, A4, A5, and A6
hold. Then, as $h,b,\hslash,\frac{\max\\{h,b\\}^{4}}{h}\to 0$ and
$\frac{nh^{3}b^{d}}{|\log(hb^{d})|},\frac{|\log(hb^{d})|}{\log\log
n},\frac{n\hslash}{|\log\hslash|},\frac{|\log\hslash|}{\log\log n}\to\infty$,
we know that

$\displaystyle\sup_{t\in\mathcal{T}^{\prime}}\left|\widehat{\theta}_{C}(t)-\theta_{C}(t)\right|$
$\displaystyle=O\left(h^{q}+b^{2}+\frac{\max\\{b,h\\}^{4}}{h}\right)+O_{P}\left(\sqrt{\frac{|\log(hb^{d})|}{nh^{3}b^{d}}}+\hslash^{2}+\sqrt{\frac{|\log\hslash|}{n\hslash}}\right)$

and

$\displaystyle\sup_{t\in\mathcal{T}^{\prime}}\left|\widehat{m}_{\theta}(t)-m(t)\right|$
$\displaystyle=O_{P}\left(\frac{1}{\sqrt{n}}\right)+O\left(h^{q}+b^{2}+\frac{\max\\{b,h\\}^{4}}{h}\right)$
$\displaystyle\quad+O_{P}\left(\sqrt{\frac{|\log(hb^{d})|}{nh^{3}b^{d}}}+\hslash^{2}+\sqrt{\frac{|\log\hslash|}{n\hslash}}\right).$

The pointwise rates of convergence for $\widehat{\theta}_{C}(t)$ and
$\widehat{m}_{\theta}(t)$ will remain the same as stated in Theorem 4 due to
the randomness of observations near the boundary of density support
$\mathcal{E}$; see the proof of Theorem 4 in Section B.4 for detailed
arguments. We restrict the uniform consistency results in Theorem 4 to a
compact set $\mathcal{T}^{\prime}\subset\mathcal{T}$ in order to avoid the
density decay of $p_{T}(t)$ near the boundary of its support $\mathcal{T}$. If
$p_{T}(t)$ is uniformly bounded away from 0 in its support, we can take
$\mathcal{T}^{\prime}=\mathcal{T}$. The (uniform) rate of convergence of
$\widehat{m}_{\theta}(t)$ consists of two parts. The first part
$O_{P}\left(\frac{1}{\sqrt{n}}\right)$ comes from the simple sample average of
$Y_{i},i=1,...,n$ and is asymptotically negligible. The second dominant part
is due to the integral component

$\widehat{\Delta}_{h,b}(t)=\frac{1}{n}\sum_{i=1}^{n}\int_{\widetilde{t}=T_{i}}^{\widetilde{t}=t}\widehat{\theta}_{C}(\widetilde{t})d\widetilde{t},$
(16)

whose rate of convergence is determined by $\widehat{\theta}_{C}(t)$.
Furthermore, the (uniform) rate of convergence of $\widehat{\theta}_{C}(t)$
comprises the rate
$O\left(h^{q}+b^{2}+\frac{\max\\{b,h\\}^{4}}{h}\right)+O_{P}\left(\sqrt{\frac{|\log(hb^{d})|}{nh^{3}b^{d}}}\right)$
for estimating the partial derivative $\frac{\partial}{\partial
t}\mu(t,\bm{s})$ and the rate
$O(\hslash^{2})+O_{P}\left(\sqrt{\frac{|\log\hslash|}{n\hslash}}\right)$ for
estimating the conditional CDF $\mbox{$\mathrm{P}$}(\bm{s}|t)$.

###  4.3 Validity of the Bootstrap Inference

Before proving the consistency of our bootstrap procedure in Section 3.3, we
first derive the asymptotic linearity of our integral estimator (8) and
localized derivative estimator (12) as intermediate results.

######  Lemma 5 (Asymptotic linearity).

Let $q\geq 2$ in the local polynomial regression for estimating
$\frac{\partial}{\partial t}\mu(t,\bm{s})$ and
$\mathcal{T}^{\prime}\subset\mathcal{T}$ be a compact set so that $p_{T}(t)$
is uniformly bounded away from 0 within $\mathcal{T}^{\prime}$. Suppose that
Assumptions A1, A2, A3, A4, A5, and A6 hold. Then, if $h\asymp b\asymp
n^{-\frac{1}{\gamma}}$ and $\hslash\asymp n^{-\frac{1}{\varpi}}$ for some
$\gamma\geq\varpi>0$ such that $\frac{nh^{d+5}}{\log n}\to c_{1}$ and
$\frac{n\hslash^{5}}{\log n}\to c_{2}$ for some finite number $c_{1},c_{2}\geq
0$ and $\frac{\log n}{n\hslash^{2}},\frac{h^{d+3}\log
n}{\hslash},\frac{h^{d+3}}{\hslash^{2}}\to 0$ as $n\to\infty$, then for any
$t\in\mathcal{T}^{\prime}$, we have that

$\sqrt{nh^{3}b^{d}}\left[\widehat{\theta}_{C}(t)-\theta_{C}(t)\right]=\mathbb{G}_{n}\bar{\varphi}_{t}+o_{P}(1)\quad\text{
and
}\quad\sqrt{nh^{3}b^{d}}\left[\widehat{m}_{\theta}(t)-m(t)\right]=\mathbb{G}_{n}\varphi_{t}+o_{P}(1),$

where
$\bar{\varphi}_{t}(Y,T,\bm{S})=\mathbb{E}_{(T_{i_{3}},\bm{S}_{i_{3}})}\left[\frac{\bm{e}_{2}^{T}\bm{M}_{q}^{-1}\bm{\Psi}_{t,\bm{S}_{i_{3}}}\left(Y,T,\bm{S}\right)}{\sqrt{hb^{d}}\cdot
p(t,\bm{S}_{i_{3}})\cdot
p_{T}(t)}\cdot\frac{1}{\hslash}\bar{K}_{T}\left(\frac{t-T_{i_{3}}}{\hslash}\right)\right]$
and

$\displaystyle\begin{split}\varphi_{t}\left(Y,T,\bm{S}\right)&=\mathbb{E}_{T_{i_{2}}}\left[\int_{T_{i_{2}}}^{t}\bar{\varphi}_{\widetilde{t}}(Y,T,\bm{S})\,d\widetilde{t}\right]\\\
&=\mathbb{E}_{T_{i_{2}}}\left\\{\int_{T_{i_{2}}}^{t}\mathbb{E}_{(T_{i_{3}},\bm{S}_{i_{3}})}\left[\frac{\bm{e}_{2}^{T}\bm{M}_{q}^{-1}\bm{\Psi}_{\widetilde{t},\bm{S}_{i_{3}}}\left(Y,T,\bm{S}\right)}{\sqrt{hb^{d}}\cdot
p(\widetilde{t},\bm{S}_{i_{3}})\cdot
p_{T}(\widetilde{t})}\cdot\frac{1}{\hslash}\bar{K}_{T}\left(\frac{\widetilde{t}-T_{i_{3}}}{\hslash}\right)\right]d\widetilde{t}\right\\}.\end{split}$
(17)

Furthermore, we have the following uniform results as:

$\displaystyle\left|\sqrt{nh^{3}b^{d}}\sup_{t\in\mathcal{T}^{\prime}}\left|\widehat{\theta}_{C}(t)-\theta_{C}(t)\right|-\sup_{t\in\mathcal{T}^{\prime}}|\mathbb{G}_{n}\bar{\varphi}_{t}|\right|=O_{P}\left(\frac{\log
n}{\sqrt{n\hslash}}+\sqrt{nh^{d+7}}+\sqrt{\frac{h^{d+3}\log
n}{\hslash}}\right)$

and

$\displaystyle\left|\sqrt{nh^{3}b^{d}}\sup_{t\in\mathcal{T}^{\prime}}\left|\widehat{m}_{\theta}(t)-m(t)\right|-\sup_{t\in\mathcal{T}^{\prime}}|\mathbb{G}_{n}\varphi_{t}|\right|$
$\displaystyle=O_{P}\left(\sqrt{nh^{d+7}}+\sqrt{\frac{\log
n}{n\hslash^{2}}}+\sqrt{\frac{h^{d+3}\log
n}{\hslash}}+\sqrt{\frac{h^{d+3}}{\hslash^{2}}}\right).$

The proof of Lemma 5 is in Section B.5, in which our key argument is to write
$\widehat{m}_{\theta}(t)$ and $\widehat{\theta}_{C}(t)$ in the form of
V-statistics (Shieh, 2014). We make two remarks to this crucial lemma.

######  Remark 3 (Non-degeneracy and validity of pointwise confidence
intervals).

We prove in Lemma 13 of Section B.6 that the variances
$\text{Var}\left[\varphi_{t}(Y,T,\bm{S})\right]$ and
$\text{Var}\left[\bar{\varphi}_{t}(Y,T,\bm{S})\right]$ of the influence
functions in (17) are positive for each $t\in\mathcal{T}^{\prime}$ as long as
$\mathrm{Var}(\epsilon)=\sigma^{2}>0$ in model (2). Thus, the asymptotic
linearity results in Lemma 5 are non-degenerate. Similarly, one can show that
the bootstrap estimates $\widehat{m}_{\theta}^{*}(t)$ and
$\widehat{\theta}_{C}^{*}(t)$ are also asymptotically linear given the
observed data $\left\\{(Y_{i},T_{i},\bm{S}_{i})\right\\}_{i=1}^{n}$. It
indicates that the pointwise bootstrap confidence intervals for $m(t)$ and
$\theta(t)$ in Section 3.3 are asymptotically valid; see Lemma 23.3 in van der
Vaart (1998) and related results in Arcones and Gine (1992); Tang and Westling
(2024).

######  Remark 4 (Bandwidth selection).

In order for those remainder terms of the asymptotically linear forms in Lemma
5 to be of the order $o_{P}(1)$, we can choose the bandwidths $h\asymp b$ to
be of the order $O\left(n^{-\frac{1}{d+5}}\right)$ and $\hslash$ to be of the
order $O\left(n^{-\frac{1}{5}}\right)$, both of which match the outputs by the
usual bandwidth selection methods for nonparametric regression (Wasserman,
2006; Schindler, 2011). In other words, we can tune $h,b,\hslash$ via the
standard bandwidth selectors for estimating the regression function
$\mu(t,\bm{s})$ and conditional CDF $p(\bm{s}|t)$ without any explicit
undersmoothing.

Apart from asymptotic linearity, we establish couplings between
$\sqrt{nh^{3}b^{d}}\sup_{t\in\mathcal{T}^{\prime}}\left|\widehat{m}_{\theta}(t)-m(t)\right|$
and $\sup_{t\in\mathcal{T}^{\prime}}\left|\mathbb{G}_{n}(\varphi_{t})\right|$
(and similarly, for $\widehat{\theta}_{C}(t)$) in Lemma 5. These coupling
results serve as key ingredients for deriving the Gaussian approximations for
$\widehat{m}_{\theta}(t)$ and $\widehat{\theta}_{C}(t)$ in Theorem 6 below and
their bootstrap consistencies. Consider two function classes

$\mathcal{F}=\left\\{(v,x,\bm{z})\mapsto\varphi_{t}(v,x,\bm{z}):t\in\mathcal{T}^{\prime}\right\\}\;\text{
and
}\;\mathcal{F}_{\theta}=\left\\{(v,x,\bm{z})\mapsto\bar{\varphi}_{t}(v,x,\bm{z}):t\in\mathcal{T}^{\prime}\right\\}$
(18)

with $\varphi_{t},\bar{\varphi}_{t}$ defined in (17) respectively. Let
$\mathbb{B},\bar{\mathbb{B}}$ be two Gaussian processes indexed by
$\mathcal{F}$ and $\mathcal{F}_{\theta}$ respectively with zero means and
covariance functions

$\displaystyle\mathrm{Cov}(\mathbb{B}(f_{1}),\mathbb{B}(f_{2}))=\mathbb{E}\left[f_{1}(\bm{U})\cdot
f_{2}(\bm{U})\right]\;\text{ and
}\;\mathrm{Cov}(\bar{\mathbb{B}}(g_{1}),\bar{\mathbb{B}}(g_{2}))=\mathbb{E}\left[g_{1}(\bm{U})\cdot
g_{2}(\bm{U})\right]$

with $\bm{U}=(Y,T,\bm{S})$ for any $f_{1},f_{2}\in\mathcal{F}$ and
$g_{1},g_{2}\in\mathcal{F}_{\theta}$.

######  Theorem 6 (Gaussian approximation).

Let $q\geq 2$ in the local polynomial regression for estimating
$\frac{\partial}{\partial t}\mu(t,\bm{s})$ and
$\mathcal{T}^{\prime}\subset\mathcal{T}$ be a compact set so that $p_{T}(t)$
is uniformly bounded away from 0 within $\mathcal{T}^{\prime}$. Suppose that
Assumptions A1, A2, A3, A4, A5, and A6 hold. If $h\asymp b\asymp
n^{-\frac{1}{\gamma}}$ and $\hslash\asymp n^{-\frac{1}{\varpi}}$ for some
$\gamma\geq\varpi>0$ such that $\frac{nh^{d+5}}{\log n}\to c_{1}$ and
$\frac{n\hslash^{5}}{\log n}\to c_{2}$ for some finite number $c_{1},c_{2}\geq
0$ and $\frac{n\hslash^{2}}{\log n},\frac{\hslash}{h^{d+3}\log n},\hslash
n^{\frac{1}{4}},\frac{\hslash^{2}}{h^{d+3}}\to\infty$ as $n\to\infty$, then
there exist Gaussian processes $\mathbb{B},\bar{\mathbb{B}}$ such that

$\displaystyle\sup_{u\geq
0}\left|\mbox{$\mathrm{P}$}\left(\sqrt{nh^{3}b^{d}}\cdot\sup_{t\in\mathcal{T}^{\prime}}\left|\widehat{m}_{\theta}(t)-m(t)\right|\leq
u\right)-\mbox{$\mathrm{P}$}\left(\sup_{f\in\mathcal{F}}|\mathbb{B}(f)|\leq
u\right)\right|=O\left(\left(\frac{\log^{5}n}{nh^{d+3}}\right)^{\frac{1}{8}}\right),$
$\displaystyle\sup_{u\geq
0}\left|\mbox{$\mathrm{P}$}\left(\sqrt{nh^{3}b^{d}}\cdot\sup_{t\in\mathcal{T}^{\prime}}\left|\widehat{\theta}_{C}(t)-\theta_{C}(t)\right|\leq
u\right)-\mbox{$\mathrm{P}$}\left(\sup_{g\in\mathcal{F}_{\theta}}|\bar{\mathbb{B}}(g)|\leq
u\right)\right|=O\left(\left(\frac{\log^{5}n}{nh^{d+3}}\right)^{\frac{1}{8}}\right),$

where $\mathcal{F},\mathcal{F}_{\theta}$ are defined in (18).

The proof of Theorem 6 is in Section B.6. It demonstrates that the
distributions of
$\sqrt{nh^{3}b^{d}}\cdot\sup_{t\in\mathcal{T}^{\prime}}\left|\widehat{m}_{\theta}(t)-m(t)\right|$
and
$\sqrt{nh^{3}b^{d}}\cdot\sup_{t\in\mathcal{T}^{\prime}}\left|\widehat{\theta}_{C}(t)-\theta_{C}(t)\right|$
can be asymptotically approximated by the suprema of two separate Gaussian
processes respectively. To establish the consistency of our bootstrap
inference procedure in Section 3.3, it remains to show that the bootstrap
versions of the differences
$\sqrt{nh^{3}b^{d}}\cdot\sup_{t\in\mathcal{T}^{\prime}}\left|\widehat{m}_{\theta}^{*}(t)-\widehat{m}_{\theta}(t)\right|$
and
$\sqrt{nh^{3}b^{d}}\cdot\sup_{t\in\mathcal{T}^{\prime}}\left|\widehat{\theta}_{C}^{*}(t)-\widehat{\theta}_{C}(t)\right|$,
conditioning on the observed data
$\bm{U}_{n}=\left\\{(Y_{i},T_{i},\bm{S}_{i})\right\\}_{i=1}^{n}$, can be
asymptotically approximated by the suprema of the same Gaussian processes,
which are summarized in the following theorem.

######  Theorem 7 (Bootstrap consistency).

Let $\mathbb{U}_{n}=\left\\{(Y_{i},T_{i},\bm{S}_{i})\right\\}_{i=1}^{n}$ be
the observed data. Under the same setup of Theorem 6, we have that

$\displaystyle\sup_{u\geq
0}\left|\mbox{$\mathrm{P}$}\left(\sqrt{nh^{3}b^{d}}\cdot\sup_{t\in\mathcal{T}^{\prime}}\left|\widehat{m}_{\theta}^{*}(t)-\widehat{m}_{\theta}(t)\right|\leq
u\Big{|}\mathbb{U}_{n}\right)-\mbox{$\mathrm{P}$}\left(\sup_{f\in\mathcal{F}}|\mathbb{B}(f)|\leq
u\right)\right|=O_{P}\left(\left(\frac{\log^{5}n}{nh^{d+3}}\right)^{\frac{1}{8}}\right)$
$\displaystyle\sup_{u\geq
0}\left|\mbox{$\mathrm{P}$}\left(\sqrt{nh^{3}b^{d}}\cdot\sup_{t\in\mathcal{T}^{\prime}}\left|\widehat{\theta}_{C}^{*}(t)-\widehat{\theta}_{C}(t)\right|\leq
u\Big{|}\mathbb{U}_{n}\right)-\mbox{$\mathrm{P}$}\left(\sup_{g\in\mathcal{F}_{\theta}}|\bar{\mathbb{B}}(g)|\leq
u\right)\right|=O_{P}\left(\left(\frac{\log^{5}n}{nh^{d+3}}\right)^{\frac{1}{8}}\right),$

where $\widehat{m}_{\theta}^{*}(t)$ and $\widehat{\theta}_{C}^{*}(t)$ are the
integral estimator (8) and localized derivative estimator (12) based on a
bootstrap sample
$\mathbb{U}_{n}^{*}=\left\\{\left(Y_{i}^{*},T_{i}^{*},\bm{S}_{i}^{*}\right)\right\\}_{i=1}^{n}$.

The proof of Theorem 7 is in Section B.7. These results, together with Theorem
6, imply the asymptotic validity of bootstrap uniform confidence bands in
Section 3.3; see Corollary 8 below with its proof in Section B.8.

######  Corollary 8 (Uniform confidence band).

Under the setup of Theorem 7, we have that

$\displaystyle\mbox{$\mathrm{P}$}\left(\theta(t)\in\left[\widehat{\theta}_{C}(t)-\bar{\xi}_{1-\alpha}^{*},\widehat{\theta}_{C}(t)+\bar{\xi}_{1-\alpha}^{*}\right]\text{
for all
}t\in\mathcal{T}^{\prime}\right)=1-\alpha+O\left(\left(\frac{\log^{5}n}{nh^{d+3}}\right)^{\frac{1}{8}}\right),$
$\displaystyle\mbox{$\mathrm{P}$}\left(m(t)\in\left[\widehat{m}_{\theta}(t)-\xi_{1-\alpha}^{*},\widehat{m}_{\theta}(t)+\xi_{1-\alpha}^{*}\right]\text{
for all
}t\in\mathcal{T}^{\prime}\right)=1-\alpha+O\left(\left(\frac{\log^{5}n}{nh^{d+3}}\right)^{\frac{1}{8}}\right).$

##  5 Experiments

In this section, we evaluate the finite-sample performances of our proposed
integral estimator (8) and localized derivative estimator (12) through several
simulation studies. We also apply them to analyzing the causal effects of fine
particulate matter on the cardiovascular mortality rate as a case study.

###  5.1 Parameter Setup

Throughout the experiments, we use the Epanechnikov kernel for $K_{T}$ and
$K_{S}$ (with the product kernel technique) in the (partial) local quadratic
regression (10). To choose its default bandwidth parameters $h,b>0$, we
implement the rule-of-thumb method in Appendix A of Yang and Tschernig (1999)
as:

$h_{ROT}=C_{h}\left[\frac{\left(\nu_{0}^{(T)}\right)^{2}\widehat{R}\left(T_{(n)}-T_{(1)}\right)}{4\left(\kappa_{2}^{(T)}\right)^{2}\widehat{C}_{\mu}}\right]^{\frac{1}{5}}n^{-\frac{1}{d+5}},\;\bm{b}_{ROT}=C_{b}\left[\frac{d\left(\nu_{0,1}^{(S)}\right)^{2d}\widehat{R}\left(\bm{S}_{(n)}-\bm{S}_{(1)}\right)}{4n\left(\kappa_{2,1}^{(S)}\right)^{2}\widehat{C}_{\mu}}\right]^{-\frac{1}{d+5}},$

where $\nu_{0}^{(T)},\nu_{0,1}^{(S)},\kappa_{2}^{(T)},\kappa_{2,1}^{(S)}>0$
are defined in Assumption A6(a),
$\widehat{R}=\frac{1}{n-5}\sum_{i=1}^{n}\left[Y-\widehat{m}(T_{i},\bm{S}_{i})\right]^{2}$
with $\widehat{m}(t,\bm{s})=\widehat{\beta}_{1}(t,\bm{s})$ from the local
polynomial regression (10) with $q=4$ and $K_{T}=K_{S}\equiv 1$ (_i.e._ ,
global (partial) forth-order polynomial fitting), $\bm{S}_{(n)}-\bm{S}_{(1)}$
are the differences between coordinatewise maxima and minima that approximate
the integration ranges, and $\widehat{C}_{\mu,T},\widehat{C}_{\mu,\bm{S}}$ are
the estimated density-weighted curvatures (or second-order partial
derivatives) of $\mu$. We obtain
$\widehat{C}_{\mu,T},\widehat{C}_{\mu,\bm{S}}$ by fitting coordinatewise
global fourth-order polynomials and averaging the estimates over the observed
data $\left\\{(T_{i},\bm{S}_{i}),i=1,...,n\right\\}$; see also Section 4.3 in
García Portugués (2023). Here, we set $\bm{b}_{ROT}\in[0,\infty)^{d}$ as a
bandwidth vector instead of a scalar $b>0$ so that it can be more adaptive to
the scale of each covariate. In addition, unless stated otherwise, we set the
scaling constants to be $C_{h}=7,C_{b}=3$ for all the experiments. Finally, to
maintain fair comparisons between our proposed estimators and the usual
regression adjustment estimators $\widehat{m}_{RA}(t)$ in (6) and
$\widehat{\theta}_{RA}(t)$ in Remark 1, we obtain $\widehat{\mu}(t,\bm{s})$
and $\widehat{\beta}_{2}(t,\bm{s})$ using the (partial) local quadratic
regression (10) with the same choices of bandwidth parameters $h,b$.

As for the conditional CDF estimator (11), we leverage the Gaussian kernel for
$\bar{K}_{T}$ and set its default bandwidth parameter $\hslash$ to the normal
reference rule in Chacón et al. (2011); Chen et al. (2016) as
$\hslash_{NR}=\left(\frac{4}{3n}\right)^{\frac{1}{5}}\widehat{\sigma}_{T}$,
where $\widehat{\sigma}_{T}$ is the sample standard deviation of
$\\{T_{i},i=1,...,n\\}$. This rule is obtained by assuming a normal
distribution of $p_{T}(t)$ and optimizing the asymptotic bias-variance trade-
off, which could potentially oversmooth the conditional CDF estimator
(Sheather, 2004). However, since our proposed estimators of $m(t)$ and
$\theta(t)$ are not very sensitive to the choice of $\hslash$ and the order of
$\hslash_{NR}$ aligns with our theoretical requirement (see Remark 4), we
would stick to $\hslash_{NR}$ for $\hslash$ in the subsequent analysis.

As for the bootstrap inference, we set the resampling time $B=1000$ and the
significance level $\alpha=0.05$, _i.e._ , targeting at 95% confidence
intervals and uniform bands for inferring the true dose-response curve $m(t)$
and its derivative $\theta(t)$.

###  5.2 Simulation Studies

We consider three different model settings under the additive confounding
model (3) for our simulation studies. For each model setting, we generate
$\\{(Y_{i},T_{i},\bm{S}_{i})\\}_{i=1}^{n}$ with $n=2000$ i.i.d. observations.

$\bullet$ Single Confounder Model: The data-generating model is described in
(1).

$\bullet$ Linear Confounding Model: We consider the following linear effect
model with

$Y=T+6S_{1}+6S_{2}+\epsilon,\quad
T=2S_{1}+S_{2}+E,\quad\text{and}\quad\bm{S}=(S_{1},S_{2})\sim\text{Uniform}[-1,1]^{2}\subset\mathbb{R}^{2},$
(19)

where $E\sim\text{Uniform}[-0.5,0.5]$ and $\epsilon\sim\mathcal{N}(0,1)$. The
marginal supports of $T$ and $\bm{S}$ are $\mathcal{T}=[-3.5,3.5]$ and
$\mathcal{S}=[-1,1]^{2}$ respectively, and the support of conditional density
$p(t|\bm{s})$ for any $\bm{s}\in\mathcal{S}$ is much narrower than
$\mathcal{T}$.

Figure 2: Simulation results under the linear confounding model (19). Left:
The estimated dose-response curves by usual regression adjustment and our
proposed estimators overlain with the true one $m(t)$. Right: The estimated
derivatives of the dose-response curve by usual regression adjustment and our
proposed estimators overlain with the true one $\theta(t)$. The middle and
right panels also present the 95% confidence intervals and/or uniform
confidence bands from our proposed estimators as shaded regions and dashed
lines, respectively.

$\bullet$ Nonlinear Confounding Model: We consider the following nonlinear
effect model with

$\displaystyle\begin{split}&Y=T^{2}+T+10Z+\epsilon,\quad T=\cos\left(\pi
Z^{3}\right)+\frac{Z}{4}+E,\quad Z=4S_{1}+S_{2},\\\ &\quad\text{ and
}\quad\bm{S}=(S_{1},S_{2})\sim\text{Uniform}[-1,1]^{2}\subset\mathbb{R}^{2},\end{split}$
(20)

where $E\sim\text{Uniform}[-0.1,0.1]$ and $\epsilon\sim\mathcal{N}(0,1)$.
Again, due to the nonlinear confounding effect $\cos(\pi Z)$ and small
treatment variation $E$, the positivity condition (A0) fails on the generated
data from (20). We also note that the debiased estimator of $m(t)$ proposed by
Takatsu and Westling (2022) (and similarly, Kennedy et al. 2017) is incapable
of recovering $m(t)$ under this data model, because some of their pseudo-
outcomes have invalid values caused by nearly zero estimated conditional
densities of $p(t|\bm{s})$.

Figure 3: Simulation results under the nonlinear confounding model (20). See
Fig 2 for the description of each panel.

Because of the presence of confounding variables, simply regressing $Y_{i}$’s
to $T_{i}$’s in the generated data from each of the above models will give
rise to biased estimates of $m(t)$. We apply both the proposed integral
estimator $\widehat{m}_{\theta}(t)$ and localized derivative estimator
$\widehat{\theta}_{C}(t)$ with bootstrap inferences in Section 3.3 to the
generated data respectively. For comparisons, we also implement the usual
regression adjustment estimators using the (partial) local quadratic
regression (10) under the same choices of bandwidth parameters. The results
are shown in Fig 1, Fig 2, and Fig 3. Note that the regression adjustment
estimators are extremely unstable, especially around the boundary of the
marginal support $\mathcal{T}$. Due to the failure of (A0), fine-tuning the
bandwidth parameters $h,b$ will not ameliorate the stability and consistency
of regression adjustment estimators. On the contrary, our proposed estimators
can recover the true dose-response curves and its derivatives with relatively
narrow pointwise confidence intervals and uniform bands.

###  5.3 Case Study: Effect of PM2.5 on Cardiovascular Mortality Rate

Air pollution, especially fine particulate matter with diameters less than 2.5
$\unit{}$ (PM2.5), is known to contribute to an increase in cardiovascular
diseases (Brook et al., 2010; Krittanawong et al., 2023). Some recent studies
also identify a positive association between the PM2.5 level ($\unit{}/m^{3}$)
and the county-level cardiovascular mortality rate (CMR; deaths/100,000
person-years) in the United States after controlling for socioeconomic factors
(Wyatt et al., 2020a).

To showcase the applicability of our proposed integral and its derivative
estimators, we apply them to analyzing the PM2.5 and CMR data in Wyatt et al.
(2020b). The data contain average annual CMR for outcome $Y$ and PM2.5
concentration for treatment/exposure $T$ from 1990 to 2010 within $n=2132$
counties, which were obtained from the US National Center for Health
Statistics and Community Multiscale Air Quality modeling system, respectively.
Our covariate vector $\bm{S}\in\mathbb{R}^{10}$ comprises two parts of
confounding variables. The first part consists of two spatial confounding
variables, latitude and longitude of each county, that helps incorporate the
spatial dependence. The second part includes eight county-level socioeconomic
factors acquired from the US census: population in 2000, civilian unemployment
rate, median household income, percentage of female households with no spouse,
percentage of vacant housing units, percentage of owner occupied housing
units, percentage of high school graduates or above, and percentage of
households below poverty. For each county in a given year, we use the closest
records from a US census year (1990, 2000, or 2010) as the values of its
socioeconomic variables and then average the values of $Y,T,\bm{S}$ for each
county over these 21 years as the final data. Finally, similar to Takatsu and
Westling (2022), we focus on the values of PM2.5 between 2.5 $\unit{}/m^{3}$
and 11.5 $\unit{}/m^{3}$ to avoid boundary effects.

We apply our proposed integral estimator (8) and localized derivative
estimator (12) to the final data with the choices of bandwidth parameters in
Section 5.1. To take the magnitude differences of treatment $T$ and covariates
in $\bm{S}$ into consideration, we multiply the (coordinatewise) standard
deviations of $T$ and $\bm{S}$ on the data to $h_{ROT},\bm{b}_{ROT}$
respectively and set $C_{h}=7,C_{b}=20$ to smooth out the estimated
derivatives. To study the effects of confounding, we also consider regressing
$Y$ on $T$ only via local quadratic regression estimators as well as fitting
$Y$ on $T$ and spatial locations only via our proposed estimators. The results
in the left panel of Fig 4 demonstrate that the effects of confounding. Before
controlling any confounding variables, the fitted curve is not monotonic and
peaks at around 8 $\unit{}/m^{3}$. After controlling for spatial locations, it
becomes flatter but still decreases when PM2.5 level is above 8
$\unit{}/m^{3}$. Only when we incorporate both spatial and socioeconomic
covariates can the estimated relationship between PM2.5 and CMR becomes
monotonically increasing. The estimated changing rate of CMR with respect to
PM2.5 in the right panel of Fig 4 reveals the same conclusion from a different
angle. More interestingly, after we adjust for all available confounding
variables, the changing rate of CMR and its 95% confidence intervals flatten
out and are unanimously above 0 when the PM2.5 level is below 9
$\unit{}/m^{3}$. This indicates a strong signal of a positive association
between the PM2.5 level and CMR. However, when the PM2.5 level is above 10
$\unit{}/m^{3}$, it is unclear whether PM2.5 still positively contributes to
cardiovascular mortality, possibly due to other competing risks (Leiser et
al., 2019). Our findings generally align with the conclusions in Wyatt et al.
(2020a); Takatsu and Westling (2022), but the superiority of our proposed
estimators comes from two aspects. First, we estimate not only the
relationship between PM2.5 and CMR but also its derivative function, providing
new insights into their association. Second, our proposed estimators are
capable of resolving the potential violation of the positivity condition (A0)
and thus lead to convincing conclusions in this observational study.

Figure 4: Estimated relationships between the PM2.5 concentration and CMR or
its changing rate at the county level. Left: The estimated CMR with respect to
the PM2.5 concentration. Right: The estimated changing rates of CMR with
respect to the PM2.5 concentration. We also present the 95% confidence
intervals and uniform confidence bands as shaded regions and dashed lines
respectively for each regression scenario.

##  6 Discussion

In summary, this paper studies nonparametric inference on the dose-response
curve and its derivative function via innovative integral and localized
derivative estimators. The main advantage of our proposed estimators is that
they can consistently recover and infer the dose-response curve and its
derivative without assuming the unrealistically strong positivity condition
(A0) under the context of continuous treatments. We establish the consistency
and bootstrap validity of the proposed estimators when nuisance functions are
estimated by kernel smoothing methods but without requiring explicit
undersmoothing on bandwidth parameters. Simulation studies and real-world
applications demonstrate the effectiveness of our proposed estimators in
addressing the failure of (A0). There are several future directions that can
further advance the impacts of our proposed estimators.

1\. Estimation of nuisance parameters/functions: The proposed integral
estimator (8) and localized derivative estimator (12) require us to estimate
two nuisance functions, the partial derivative $\frac{\partial}{\partial
t}\mu(t,\bm{s})$ and the conditional CDF $P(\bm{s}|t)$. When we focus on the
kernel smoothing methods, the finite-sample performances of our integral
estimator (8) and localized derivative estimator (12) depend on the choices of
three bandwidth parameters $h,b,\hslash>0$. We only consider the rule-of-thumb
and normal reference bandwidth selection methods in this paper, but it would
be of research interest to study how the bandwidth choices can be improved
through the plug-in rule (Ruppert et al., 1995) or cross-validation (Li and
Racine, 2004). More broadly, other estimators, such as regression splines for
$\frac{\partial}{\partial t}\mu(t,\bm{s})$ (Friedman, 1991; Zhou and Wolfe,
2000) and nearest neighbor-type or local logistic approaches for $P(\bm{s}|t)$
(Stute, 1986; Hall et al., 1999), are also be applied to estimating the
nuisance functions and may lead to better performances of our proposed
estimators $\widehat{\theta}_{C}(t)$ and $\widehat{m}_{\theta}(t)$.

2\. IPTW and doubly robust estimators: Our proposed estimators (8) and (12)
are based on the idea of regression adjustment. Naturally, it would also be
interesting to investigate how to leverage our integral and localized
techniques to resolve the positivity requirement in the existing generalized
propensity score-based approaches (Hirano and Imbens, 2004; Imai and van Dyk,
2004) and doubly robust methods (Kennedy et al., 2017; Westling et al., 2020;
Colangelo and Lee, 2020; Semenova and Chernozhukov, 2021; Bonvini and Kennedy,
2022) for estimating the dose-response curve and its derivative function.

3\. Additive model diagnostics: Assumption A2 is vital to the consistency of
our proposed estimators, and it is unclear if this assumption holds beyond the
additive confounding model (3). However, it is possible to utilize our
estimators to test the correctness of the additive structure in (3). Under the
positivity condition (A0), one can test (3) by estimating the absolute
difference of $\left|\mathbb{E}\left[\frac{\partial}{\partial
t}\mu(t,\bm{S})\right]-\mathbb{E}\left[\frac{\partial}{\partial
t}\mu(t,\bm{S})\big{|}T=t\right]\right|$ from 0 via the regression adjustment
estimator $\widehat{\theta}_{RA}(t)$ in Remark 1 and our estimator
$\widehat{\theta}_{C}(t)$ in (12). Other testing procedures, such as the
marginal integration regression method (Linton and Nielsen, 1995), are also
applicable. When the positivity condition (A0) is violated, the partial
derivative $\beta_{2}(t,\bm{s})=\frac{\partial}{\partial t}\mu(t,\bm{s})$ is
still independent of $\bm{s}$ for any $t\in\mathcal{T}$ under (3), and our
estimator $\widehat{\beta}_{2}(t,\bm{s})$ in (10) should not depend on
$\bm{s}$ asymptotically. It suggests that statistical quantities, such as
$\Theta(t)=\sup_{\bm{s}\in\mathcal{S}(t)}\left|\widehat{\beta}_{2}(t,\bm{s})-\widehat{\theta}_{C}(t)\right|$
and $\Theta=\int_{\mathcal{T}}\Theta(t)\,dt$, should converge to 0 as
$n\to\infty$. Thus, one may derive the limiting distribution of $\Theta(t)$ or
$\Theta$ to carry out a procedure for statistically testing the additive
confounding model (3) without assuming (A0).

4\. Violation of ignorability (Assumption A1(b)): In observational studies,
there could be some unmeasured confounding variables that go beyond our
covariate vector $\bm{S}\in\mathcal{S}\subset\mathbb{R}^{d}$ and bias our
inference (VanderWeele, 2008). Hence, it is important to analyze the
sensitivities of the dose-response curve and its derivative with respect to
the violation of ignorability. In this scenario, the instrumental variable
model (Kilbertus et al., 2020) and the Riesz-Frechet representation technique
(Chernozhukov et al., 2022) might be useful for the analysis. As an
alternative, one can also incorporate an additional high-dimensional covariate
vector $\bm{Z}\in\mathbb{R}^{d^{\prime}}$ into our model (2) to avoid
unmeasured confounding variables. However, since nonparametric estimations on
$m(t)$ and $\theta(t)$ have slow rates of convergence when $d^{\prime}$ is
large, one may consider imposing a linear form on $\bm{Z}$ so that

$\mathbb{E}\left(Y|T=t,\bm{S}=\bm{s},\bm{Z}=\bm{z}\right)=m(t)+\eta(\bm{s})+\vartheta^{T}\bm{z}$

or utilize some orthogonalization techniques in Section 4.3 of Rothenhäusler
and Yu (2019) to conduct valid inference. We will leave them as future works.

## Acknowledgements

YZ is supported in part by YC’s NSF grant DMS-2141808 and AG’s NSF grant
DMS-2310578. YC is supported by NSF grants DMS-1952781, 2112907, 2141808, and
NIH U24-AG07212. AG is supported by NSF grant DMS-2310578.

## References

  * Arcones and Gine (1992) [author] Arcones, Miguel A.M. A. and Gine, EvaristE. (1992). On the Bootstrap of $U$ and $V$ Statistics. The Annals of Statistics 20 655 – 674.
  * Bahadur (1966) [author] Bahadur, R RajR. R. (1966). A note on quantiles in large samples. The Annals of Mathematical Statistics 37 577–580.
  * Bashtannyk and Hyndman (2001) [author] Bashtannyk, David MD. M. and Hyndman, Rob JR. J. (2001). Bandwidth selection for kernel conditional density estimation. Computational Statistics & Data Analysis 36 279–298.
  * Bjerve et al. (1985) [author] Bjerve, SteinarS., Doksum, Kjell AK. A. and Yandell, Brian SB. S. (1985). Uniform confidence bounds for regression based on a simple moving average. Scandinavian Journal of Statistics 159–169.
  * Bonvini and Kennedy (2022) [author] Bonvini, MatteoM. and Kennedy, Edward HE. H. (2022). Fast convergence rates for dose-response estimation. arXiv preprint arXiv:2207.11825.
  * Brook et al. (2010) [author] Brook, Robert D.R. D., Rajagopalan, SanjayS., Pope, C. ArdenC. A., Brook, Jeffrey R.J. R., Bhatnagar, AruniA., Diez-Roux, Ana V.A. V., Holguin, FernandoF., Hong, YulingY., Luepker, Russell V.R. V., Mittleman, Murray A.M. A., Peters, AnnetteA., Siscovick, DavidD., Smith, Sidney C.S. C., Whitsel, LaurieL. and Kaufman, Joel D.J. D. (2010). Particulate matter air pollution and cardiovascular disease: an update to the scientific statement from the American Heart Association. Circulation 121 2331–2378.
  * Buja et al. (1989) [author] Buja, AndreasA., Hastie, TrevorT. and Tibshirani, RobertR. (1989). Linear Smoothers and Additive Models. The Annals of Statistics 17 453 – 510.
  * Busso et al. (2014) [author] Busso, MatiasM., DiNardo, JohnJ. and McCrary, JustinJ. (2014). New evidence on the finite sample properties of propensity score reweighting and matching estimators. Review of Economics and Statistics 96 885–897.
  * Calonico et al. (2018) [author] Calonico, SebastianS., Cattaneo, Matias DM. D. and Farrell, Max HM. H. (2018). On the effect of bias estimation on coverage accuracy in nonparametric inference. Journal of the American Statistical Association 113 767–779.
  * Cattaneo et al. (2010) [author] Cattaneo, Matias DM. D., Crump, Richard KR. K. and Jansson, MichaelM. (2010). Robust data-driven inference for density-weighted average derivatives. Journal of the American Statistical Association 105 1070–1083.
  * Chacón et al. (2011) [author] Chacón, José EJ. E., Duong, TarnT. and Wand, MPM. (2011). Asymptotics for general multivariate kernel density derivative estimators. Statistica Sinica 807–840.
  * Chamberlain (1986) [author] Chamberlain, GaryG. (1986). Asymptotic efficiency in semi-parametric models with censoring. Journal of Econometrics 32 189–218.
  * Chen et al. (2015) [author] Chen, Yen-ChiY.-C., Genovese, Christopher R.C. R. and Wasserman, LarryL. (2015). Asymptotic theory for density ridges. The Annals of Statistics 43 1896–1928.
  * Chen et al. (2016) [author] Chen, Yen-ChiY.-C., Genovese, Christopher R.C. R. and Wasserman, LarryL. (2016). A comprehensive approach to mode clustering. Electronic Journal of Statistics 10 210 – 241.
  * Chen et al. (2017) [author] Chen, Yen-ChiY.-C., Genovese, Christopher RC. R. and Wasserman, LarryL. (2017). Density level sets: Asymptotics, inference, and visualization. Journal of the American Statistical Association 112 1684–1696.
  * Cheng and Chen (2019) [author] Cheng, GangG. and Chen, Yen-ChiY.-C. (2019). Nonparametric inference via bootstrapping the debiased estimator. Electronic Journal of Statistics 13 2194 – 2256.
  * Chernozhukov et al. (2014) [author] Chernozhukov, VictorV., Chetverikov, DenisD. and Kato, KengoK. (2014). Gaussian approximation of suprema of empirical processes. The Annals of Statistics 42 1564–1597.
  * Chernozhukov et al. (2022) [author] Chernozhukov, VictorV., Cinelli, CarlosC., Newey, WhitneyW., Sharma, AmitA. and Syrgkanis, VasilisV. (2022). Long story short: Omitted variable bias in causal machine learning Technical Report, National Bureau of Economic Research.
  * Colangelo and Lee (2020) [author] Colangelo, KyleK. and Lee, Ying-YingY.-Y. (2020). Double debiased machine learning nonparametric inference with continuous treatments. arXiv preprint arXiv:2004.03036.
  * Cole and Hernán (2008) [author] Cole, Stephen RS. R. and Hernán, Miguel AM. A. (2008). Constructing inverse probability weights for marginal structural models. American Journal of Epidemiology 168 656–664.
  * Cox (1958) [author] Cox, David RoxbeeD. R. (1958). Planning of Experiments. Wiley.
  * Crump et al. (2009) [author] Crump, Richard KR. K., Hotz, V JosephV. J., Imbens, Guido WG. W. and Mitnik, Oscar AO. A. (2009). Dealing with limited overlap in estimation of average treatment effects. Biometrika 96 187–199.
  * Dehejia and Wahba (1999) [author] Dehejia, Rajeev HR. H. and Wahba, SadekS. (1999). Causal effects in nonexperimental studies: Reevaluating the evaluation of training programs. Journal of the American statistical Association 94 1053–1062.
  * Díaz and van der Laan (2013) [author] Díaz, IvánI. and van der Laan, Mark JM. J. (2013). Targeted data adaptive estimation of the causal dose–response curve. Journal of Causal Inference 1 171–192.
  * Dvoretzky et al. (1956) [author] Dvoretzky, A.A., Kiefer, J.J. and Wolfowitz, J.J. (1956). Asymptotic Minimax Character of the Sample Distribution Function and of the Classical Multinomial Estimator. The Annals of Mathematical Statistics 27 642 – 669.
  * D’Amour et al. (2021) [author] D’Amour, AlexanderA., Ding, PengP., Feller, AviA., Lei, LihuaL. and Sekhon, JasjeetJ. (2021). Overlap in observational studies with high-dimensional covariates. Journal of Econometrics 221 644–654.
  * Efron (1979) [author] Efron, BradleyB. (1979). Bootstrap Methods: Another Look at the Jackknife. The Annals of Statistics 7 1 – 26.
  * Einmahl and Mason (2005) [author] Einmahl, UweU. and Mason, David M.D. M. (2005). Uniform in bandwidth consistency of kernel-type function estimators. The Annals of Statistics 33 1380 – 1403.
  * Fan and Gijbels (1996) [author] Fan, JianqingJ. and Gijbels, IreneI. (1996). Local polynomial modelling and its applications 66\. Chapman & Hall/CRC.
  * Fan and Guerre (2015) [author] Fan, YanqinY. and Guerre, EmmanuelE. (2015). Multivariate Local Polynomial Estimators: Boundary Properties and Uniform Asymptotic Linear Representation. Advances in Econometrics.
  * Fan et al. (1996) [author] Fan, JianqingJ., Gijbels, IrèneI., Hu, Tien-ChungT.-C. and Huang, Li-ShanL.-S. (1996). A study of variable bandwidth selection for local polynomial regression. Statistica Sinica 6 113–127.
  * Flores (2007) [author] Flores, CarlosC. (2007). Estimation of dose-response functions and optimal doses with a continuous treatment Technical Report, Department of Economics, University of Miami.
  * Freedman (1981) [author] Freedman, David AD. A. (1981). Bootstrapping regression models. The Annals of Statistics 9 1218–1228.
  * Friedman (1991) [author] Friedman, Jerome HJ. H. (1991). Multivariate adaptive regression splines. The Annals of Statistics 19 1–67.
  * García Portugués (2023) [author] García Portugués, EduardoE. (2023). Notes for nonparametric statistics. (Accessed on March 26, 2024).
  * Giessing (2023) [author] Giessing, AlexanderA. (2023). Gaussian and Bootstrap Approximations for Suprema of Empirical Processes. arXiv preprint arXiv:2309.01307.
  * Gilbert et al. (2023) [author] Gilbert, BrianB., Datta, AbhirupA., Casey, Joan A.J. A. and Ogburn, Elizabeth L.E. L. (2023). A causal inference framework for spatial confounding. arXiv preprint arXiv:2112.14946.
  * Gill and Robins (2001) [author] Gill, Richard DR. D. and Robins, James MJ. M. (2001). Causal inference for complex longitudinal data: the continuous case. Annals of Statistics 29 1785–1811.
  * Giné and Guillou (2002) [author] Giné, EvaristE. and Guillou, ArmelleA. (2002). Rates of strong uniform consistency for multivariate kernel density estimators. Annales de l’Institut Henri Poincare (B) Probability and Statistics 38 907–921.
  * Hall (1992) [author] Hall, PeterP. (1992). On bootstrap confidence intervals in nonparametric regression. The Annals of Statistics 20 695–711.
  * Hall et al. (1999) [author] Hall, PeterP., Wolff, Rodney CLR. C. and Yao, QiweiQ. (1999). Methods for estimating a conditional distribution function. Journal of the American Statistical Association 94 154–163.
  * Härdle and Stoker (1989) [author] Härdle, WolfgangW. and Stoker, Thomas MT. M. (1989). Investigating smooth multiple regression by the method of average derivatives. Journal of the American statistical Association 84 986–995.
  * Hastie and Loader (1993) [author] Hastie, TrevorT. and Loader, CliveC. (1993). Local regression: Automatic kernel carpentry. Statistical Science 120–129.
  * Hines et al. (2023) [author] Hines, OliverO., Diaz-Ordaz, KarlaK. and Vansteelandt, StijnS. (2023). Optimally weighted average derivative effects. arXiv preprint arXiv:2308.05456.
  * Hirano and Imbens (2004) [author] Hirano, KeisukeK. and Imbens, Guido WG. W. (2004). The propensity score with continuous treatments. Applied Bayesian modeling and causal inference from incomplete-data perspectives 226164 73–84.
  * Hirshberg and Wager (2020) [author] Hirshberg, David AD. A. and Wager, StefanS. (2020). Debiased inference of average partial effects in single-index models: Comment on wooldridge and zhu. Journal of Business & Economic Statistics 38 19–24.
  * Holmes et al. (2012) [author] Holmes, Michael PM. P., Gray, Alexander GA. G. and Isbell, Charles LeeC. L. (2012). Fast nonparametric conditional density estimation. arXiv preprint arXiv:1206.5278.
  * Huber et al. (2020) [author] Huber, MartinM., Hsu, Yu-ChinY.-C., Lee, Ying-YingY.-Y. and Lettry, LayalL. (2020). Direct and indirect effects of continuous treatments based on generalized propensity score weighting. Journal of Applied Econometrics 35 814–840.
  * Imai and van Dyk (2004) [author] Imai, KosukeK. and van Dyk, David AD. A. (2004). Causal inference with general treatment regimes: Generalizing the propensity score. Journal of the American Statistical Association 99 854–866.
  * Kammann and Wand (2003) [author] Kammann, EEE. and Wand, Matthew PM. P. (2003). Geoadditive models. Journal of the Royal Statistical Society Series C: Applied Statistics 52 1–18.
  * Keller and Szpiro (2020) [author] Keller, Joshua PJ. P. and Szpiro, Adam AA. A. (2020). Selecting a scale for spatial confounding adjustment. Journal of the Royal Statistical Society Series A: Statistics in Society 183 1121–1143.
  * Kennedy (2019) [author] Kennedy, Edward HE. H. (2019). Nonparametric causal effects based on incremental propensity score interventions. Journal of the American Statistical Association 114 645–656.
  * Kennedy et al. (2017) [author] Kennedy, Edward HE. H., Ma, ZongmingZ., McHugh, Matthew DM. D. and Small, Dylan SD. S. (2017). Nonparametric methods for doubly robust estimation of continuous treatment effects. Journal of the Royal Statistical Society Series B: Statistical Methodology 79 1229–1245.
  * Khan and Tamer (2010) [author] Khan, ShakeebS. and Tamer, ElieE. (2010). Irregular identification, support conditions, and inverse weight estimation. Econometrica 78 2021–2042.
  * Kilbertus et al. (2020) [author] Kilbertus, NikiN., Kusner, Matt JM. J. and Silva, RicardoR. (2020). A class of algorithms for general instrumental variable models. Advances in Neural Information Processing Systems 33 20108–20119.
  * Kong et al. (2010) [author] Kong, EfangE., Linton, OliverO. and Xia, YingcunY. (2010). Uniform Bahadur representation for local polynomial estimates of M-regression and its application to the additive model. Econometric Theory 26 1529–1564.
  * Krittanawong et al. (2023) [author] Krittanawong, ChayakritC., Qadeer, Yusuf KamranY. K., Hayes, Richard B.R. B., Wang, ZhenZ., Thurston, George D.G. D., Virani, SalimS. and Lavie, Carl J.C. J. (2023). PM2.5 and cardiovascular diseases: State-of-the-Art review. International Journal of Cardiology Cardiovascular Risk and Prevention 19 200217\.
  * Léger et al. (2022) [author] Léger, MaximeM., Chatton, ArthurA., Le Borgne, FlorentF., Pirracchio, RomainR., Lasocki, SigismondS. and Foucher, YohannY. (2022). Causal inference in case of near-violation of positivity: comparison of methods. Biometrical Journal 64 1389–1403.
  * Leiser et al. (2019) [author] Leiser, Claire LC. L., Smith, Ken RK. R., VanDerslice, James AJ. A., Glotzbach, Jason PJ. P., Farrell, Timothy WT. W. and Hanson, Heidi AH. A. (2019). Evaluation of the sex-and-age-specific effects of PM2. 5 on hospital readmission in the presence of the competing risk of mortality in the medicare population of Utah 1999–2009. Journal of Clinical Medicine 8 2114\.
  * Li and Racine (2004) [author] Li, QiQ. and Racine, JeffJ. (2004). Cross-validated local linear nonparametric regression. Statistica Sinica 485–512.
  * Linton and Nielsen (1995) [author] Linton, OliverO. and Nielsen, Jens PerchJ. P. (1995). A Kernel Method of Estimating Structured Nonparametric Regression Based on Marginal Integration. Biometrika 82 93–100.
  * Lu (1996) [author] Lu, Zhan-QianZ.-Q. (1996). Multivariate locally weighted polynomial fitting and partial derivative estimation. Journal of Multivariate Analysis 59 187–205.
  * Ma and Wang (2020) [author] Ma, XinweiX. and Wang, JingshenJ. (2020). Robust inference using inverse probability weighting. Journal of the American Statistical Association 115 1851–1860.
  * Manski (1990) [author] Manski, Charles FC. F. (1990). Nonparametric bounds on treatment effects. The American Economic Review 80 319–323.
  * Massart (1990) [author] Massart, PascalP. (1990). The tight constant in the Dvoretzky-Kiefer-Wolfowitz inequality. The Annals of Probability 18 1269–1283.
  * Neugebauer and van der Laan (2007) [author] Neugebauer, RomainR. and van der Laan, MarkM. (2007). Nonparametric causal effects based on marginal structural models. Journal of Statistical Planning and Inference 137 419–434.
  * Newey (1994) [author] Newey, Whitney KW. K. (1994). Kernel estimation of partial means and a general variance estimator. Econometric Theory 10 1–21.
  * Newey and Stoker (1993) [author] Newey, Whitney KW. K. and Stoker, Thomas MT. M. (1993). Efficiency of weighted average derivative estimators and index models. Econometrica 61 1199–1223.
  * Nolan and Pollard (1987) [author] Nolan, DeborahD. and Pollard, DavidD. (1987). $U$-Processes: Rates of Convergence. The Annals of Statistics 15 780 – 799.
  * Paciorek (2010) [author] Paciorek, Christopher JC. J. (2010). The importance of scale for spatial-confounding bias and precision of spatial regression estimators. Statistical Science 25 107-125.
  * Petersen et al. (2012) [author] Petersen, Maya LM. L., Porter, Kristin EK. E., Gruber, SusanS., Wang, YueY. and Van Der Laan, Mark JM. J. (2012). Diagnosing and responding to violations in the positivity assumption. Statistical methods in medical research 21 31–54.
  * Powell et al. (1989) [author] Powell, James LJ. L., Stock, James HJ. H. and Stoker, Thomas MT. M. (1989). Semiparametric estimation of index coefficients. Econometrica 1403–1430.
  * Robins (1986) [author] Robins, JamesJ. (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 1393–1512.
  * Robins et al. (2000) [author] Robins, James MJ. M., Hernan, Miguel AngelM. A. and Brumback, BabetteB. (2000). Marginal structural models and causal inference in epidemiology. Epidemiology 11 550–560.
  * Rothe (2017) [author] Rothe, ChristophC. (2017). Robust confidence intervals for average treatment effects under limited overlap. Econometrica 85 645–660.
  * Rothenhäusler and Yu (2019) [author] Rothenhäusler, DominikD. and Yu, BinB. (2019). Incremental causal effects. arXiv preprint arXiv:1907.13258.
  * Rubin (1974) [author] Rubin, Donald BD. B. (1974). Estimating causal effects of treatments in randomized and nonrandomized studies. Journal of Educational Psychology 66 688-701.
  * Rubin (1980) [author] Rubin, Donald BD. B. (1980). Randomization analysis of experimental data: The Fisher randomization test comment. Journal of the American Statistical Association 75 591–593.
  * Ruppert et al. (1995) [author] Ruppert, DavidD., Sheather, Simon JS. J. and Wand, Matthew PM. P. (1995). An effective bandwidth selector for local least squares regression. Journal of the American Statistical Association 90 1257–1270.
  * Ruppert and Wand (1994) [author] Ruppert, DavidD. and Wand, Matthew PM. P. (1994). Multivariate locally weighted least squares regression. The Annals of Statistics 22 1346–1370.
  * Schindler (2011) [author] Schindler, AnjaA. (2011). Bandwidth selection in nonparametric kernel estimation, PhD thesis, Göttingen, Georg-August Universität, Diss.
  * Schnell and Papadogeorgou (2020) [author] Schnell, PatrickP. and Papadogeorgou, GeorgiaG. (2020). Mitigating unobserved spatial confounding when estimating the effect of supermarket access on cardiovascular disease deaths. Annals of Applied Statistics 14 2069-2095.
  * Semenova and Chernozhukov (2021) [author] Semenova, ViraV. and Chernozhukov, VictorV. (2021). Debiased machine learning of conditional average treatment effects and other causal functions. The Econometrics Journal 24 264–289.
  * Shao (2003) [author] Shao, JunJ. (2003). Mathematical Statistics. Springer Science & Business Media.
  * Sheather (2004) [author] Sheather, Simon JS. J. (2004). Density estimation. Statistical Science 19 588–597.
  * Shieh (2014) [author] Shieh, Grace SG. S. (2014). U-and V-statistics. Wiley StatsRef: Statistics Reference Online.
  * Stone (1980) [author] Stone, Charles JC. J. (1980). Optimal rates of convergence for nonparametric estimators. The Annals of Statistics 8 1348–1360.
  * Stone (1982) [author] Stone, Charles JC. J. (1982). Optimal global rates of convergence for nonparametric regression. The Annals of Statistics 10 1040–1053.
  * Stute (1986) [author] Stute, WinfriedW. (1986). Conditional empirical processes. The Annals of Statistics 14 638–647.
  * Takatsu and Westling (2022) [author] Takatsu, KentaK. and Westling, TedT. (2022). Debiased inference for a covariate-adjusted regression function. arXiv preprint arXiv:2210.06448.
  * Tang and Westling (2024) [author] Tang, ZhouZ. and Westling, TedT. (2024). Consistency of the bootstrap for asymptotically linear estimators based on machine learning. arXiv preprint arXiv:2404.03064.
  * Thaden and Kneib (2018) [author] Thaden, HaukeH. and Kneib, ThomasT. (2018). Structural equation models for dealing with spatial confounding. The American Statistician 72 239–252.
  * van der Laan and Robins (2003) [author] van der Laan, Mark JM. J. and Robins, James MJ. M. (2003). Unified methods for censored longitudinal data and causality. Springer.
  * van der Vaart (1991) [author] van der Vaart, AadA. (1991). On differentiable functionals. The Annals of Statistics 19 178–204.
  * van der Vaart (1998) [author] van der Vaart, Aad W.A. W. (1998). Asymptotic Statistics. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press.
  * van der Vaart and Wellner (1996) [author] van der Vaart, A.A. and Wellner, J. A.J. A. (1996). Weak Convergence and Empirical Processes: With Applications to Statistics. Springer Series in Statistics. Springer.
  * VanderWeele (2008) [author] VanderWeele, Tyler JT. J. (2008). The sign of the bias of unmeasured confounding. Biometrics 64 702–706.
  * Wand and Jones (1994) [author] Wand, Matt PM. P. and Jones, M ChrisM. C. (1994). Kernel smoothing. CRC press.
  * Wasserman (2006) [author] Wasserman, LarryL. (2006). All of Nonparametric Statistics. Springer Science & Business Media.
  * Waud (1975) [author] Waud, D. R.D. R. (1975). Analysis of dose-response curves. In Smooth muscle (Edwin E.E. E. Daniel and David M.D. M. Paton, eds.) 471–506. Springer US, Boston, MA.
  * Westling et al. (2020) [author] Westling, TedT., Gilbert, PeterP. and Carone, MarcoM. (2020). Causal isotonic regression. Journal of the Royal Statistical Society Series B: Statistical Methodology 82 719–747.
  * Westreich and Cole (2010) [author] Westreich, DanielD. and Cole, Stephen RS. R. (2010). Invited commentary: positivity in practice. American Journal of Epidemiology 171 674–677.
  * Wiecha and Reich (2024) [author] Wiecha, NateN. and Reich, Brian JB. J. (2024). Two-stage Spatial Regression Models for Spatial Confounding. arXiv preprint arXiv:2404.09358.
  * Wu (1986) [author] Wu, Chien-Fu JeffC.-F. J. (1986). Jackknife, bootstrap and other resampling methods in regression analysis. The Annals of Statistics 14 1261–1295.
  * Wyatt et al. (2020a) [author] Wyatt, Lauren HL. H., Peterson, Geoffrey CLG. C., Wade, Timothy JT. J., Neas, Lucas ML. M. and Rappold, Ana GA. G. (2020a). The contribution of improved air quality to reduced cardiovascular mortality: Declines in socioeconomic differences over time. Environment international 136 105430\.
  * Wyatt et al. (2020b) [author] Wyatt, Lauren HL. H., Peterson, Geoffrey Colin LG. C. L., Wade, Tim JT. J., Neas, Lucas ML. M. and Rappold, Ana GA. 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\.
  * Yang and Ding (2018) [author] Yang, ShuS. and Ding, PengP. (2018). Asymptotic inference of causal effects with observational studies trimmed by the estimated propensity scores. Biometrika 105 487–493.
  * Yang and Tschernig (1999) [author] Yang, LijianL. and Tschernig, RolfR. (1999). Multivariate bandwidth selection for local linear regression. Journal of the Royal Statistical Society Series B: Statistical Methodology 61 793–815.
  * Zhou and Wolfe (2000) [author] Zhou, ShanggangS. and Wolfe, Douglas AD. A. (2000). On derivative estimation in spline regression. Statistica Sinica 10 93–108.

Contents

[sections] [sections]l0

The supplementary material contains the proofs of main results in the paper
and other auxiliary technical results.

##  Appendix A Nonparametric Bounds Under Zero Treatment Variations

In this section, we study nonparametric bounds on the dose-response curve
$m(t)$ and its derivative function $\theta(t)$ under the additive confounding
model (3) when the treatment variation $E$ has zero variance. We will leave
the derivations of nonparametric bounds under the general model (2) as a
future direction.

###  A.1 Nonparametric Bound on $m(t)$

Recall from Section 2.2 that when $\mathrm{Var}(E)=0$, the conditional mean
outcome function $\mu(t,\bm{s})$ can only be identified on the lower
dimensional surface
$\left\\{(t,\bm{s})\in\mathcal{T}\times\mathcal{S}:t=f(\bm{s})\right\\}$ so
that

$\mu\left(f(\bm{s}),\bm{s}\right)=\bar{m}\left(f(\bm{s})\right)+\eta(\bm{s})=m\left(f(\bm{s})\right)+\eta(\bm{s}).$
(21)

Given that $\mathrm{Var}(E)=0$, we can always fit the relation $T=f(\bm{S})$
from the observed data. Hence, we assume that the function $f$ is known for
any $\bm{s}\in\mathcal{S}$ in the sequel. To obtain a nonparametric bound for
all the possible values that $m(t)$ can take, we impose the following
assumption on the magnitude of the random effect function
$\eta:\mathcal{S}\to\mathbb{R}$ in (3); see Manski (1990) for related
discussions.

######  Assumption A7 (Bounded random effect).

Let $L_{f}(t)=\left\\{\bm{s}\in\mathcal{S}:f(\bm{s})=t\right\\}$ be a level
set of the function $f:\mathcal{S}\to\mathbb{R}$ at $t\in\mathcal{T}$. There
exists a constant $\rho_{1}>0$ such that

$\rho_{1}\geq\max\left\\{\sup_{t\in\mathcal{T}}\sup_{\bm{s}\in
L_{f}(t)}|\eta(\bm{s})|,\;\frac{\sup_{t\in\mathcal{T}}\sup_{\bm{s}\in
L_{f}(t)}\mu(f(\bm{s}),\bm{s})-\inf_{t\in\mathcal{T}}\inf_{\bm{s}\in
L_{f}(t)}\mu(f(\bm{s}),\bm{s})}{2}\right\\}.$

By (21) and the first lower on
$\rho_{1}\geq\sup_{t\in\mathcal{T}}\sup_{\bm{s}\in L_{f}(t)}|\eta(\bm{s})|$ in
Assumption A7, we know that

$\left|\mu(f(\bm{s}),\bm{s})-m(t)\right|=|\eta(\bm{s})|\leq\rho_{1}$

for any $\bm{s}\in L_{f}(t)$. It also implies that

$\displaystyle\begin{split}m(t)&\in\bigcap_{\bm{s}\in
L_{f}(t)}\left[\mu(f(\bm{s}),\bm{s})-\rho_{1},\,\mu(f(\bm{s}),\bm{s})+\rho_{1}\right]\\\
&=\left[\sup_{\bm{s}\in
L_{f}(t)}\mu(f(\bm{s}),\bm{s})-\rho_{1},\,\inf_{\bm{s}\in
L_{f}(t)}\mu(f(\bm{s}),\bm{s})+\rho_{1}\right],\end{split}$ (22)

which is the nonparametric bound on $m(t)$ that contains all the possible
values of $m(t)$ for any fixed $t\in\mathcal{T}$ when $\mathrm{Var}(E)=0$.
Notice that this bound (22) is well-defined and nonempty under the second
lower bound on $\rho_{1}$ in Assumption A7.

###  A.2 Nonparametric Bound on $\theta(t)$

We already know that $\mu(f(\bm{s}),\bm{s})$ and $f(\bm{s})$ can be identified
even when $\mathrm{Var}(E)=0$, so their derivatives with respect to $\bm{s}$
are also identifiable. Hence, we assume that the functions
$\nabla_{\bm{s}}\mu(f(\bm{s}),\bm{s})$ and $\nabla f(\bm{s})$ are known for
any $\bm{s}\in\mathcal{S}$ in the sequel. The nonparametric bound on the dose-
response curve $m(t)$ involves an upper bound on the random effect function
$\eta:\mathcal{S}\to\mathbb{R}$ (recall Assumption A7), while the
nonparametric bound on the derivative function $\theta(t)$ requires some
constraints on the gradient $\nabla\eta(\bm{s})$ of the random effect function
as follows.

######  Assumption A8 (Constraints on the random effect gradient).

Let $L_{f}(t)=\left\\{\bm{s}\in\mathcal{S}:f(\bm{s})=t\right\\}$ be a level
set of the function $f:\mathcal{S}\to\mathbb{R}$ at $t\in\mathcal{T}$. There
exist constants $\rho_{2},\rho_{3}>0$ such that

$\rho_{3}\leq\inf_{t\in\mathcal{T}}\inf_{\bm{s}\in
L_{f}(t)}\left|\left|\nabla\eta(\bm{s})\right|\right|_{\min}\leq\sup_{t\in\mathcal{T}}\sup_{\bm{s}\in
L_{f}(t)}\left|\left|\nabla\eta(\bm{s})\right|\right|_{\max}\leq\rho_{2}$

and

$\sup_{t\in\mathcal{T}}\sup_{\bm{s}\in
L_{f}(t)}\max_{j=1,...,d}\left[\frac{v_{j}(\bm{s})-\mathrm{sign}(g_{j}(\bm{s}))\cdot\rho_{2}}{g_{j}(\bm{s})}\right]\leq\inf_{t\in\mathcal{T}}\inf_{\bm{s}\in
L_{f}(t)}\min_{j=1,...,d}\left[\frac{v_{j}(\bm{s})+\mathrm{sign}(g_{j}(\bm{s}))\cdot\rho_{2}}{g_{j}(\bm{s})}\right],$

where
$\left|\left|\nabla\eta(\bm{s})\right|\right|_{\min}=\min_{j=1,...,d}\left|\frac{\partial}{\partial
s_{j}}\eta(\bm{s})\right|$,
$\left|\left|\nabla\eta(\bm{s})\right|\right|_{\max}=\max_{j=1,...,d}\left|\frac{\partial}{\partial
s_{j}}\eta(\bm{s})\right|$, $v_{j}(\bm{s})=\frac{\partial}{\partial
s_{j}}\mu(f(\bm{s}),\bm{s})$, and $g_{j}(\bm{s})=\frac{\partial}{\partial
s_{j}}f(\bm{s})$ for $j=1,...,d$.

The first inequality in Assumption A8 ensures that we can derive a finite and
meaningful nonparametric bound on $\theta(t)$ for any $t\in\mathcal{T}$. In
particular, a direct calculation shows that

$\displaystyle\nabla_{\bm{s}}\mu(f(\bm{s}),\bm{s})$
$\displaystyle=\nabla_{\bm{s}}m(f(\bm{s}))+\nabla\eta(\bm{s})$
$\displaystyle=m^{\prime}(f(\bm{s}))\cdot\nabla f(\bm{s})+\nabla\eta(\bm{s})$
$\displaystyle=\theta(t)\cdot\nabla f(\bm{s})+\nabla\eta(\bm{s})$

for all $\bm{s}\in L_{f}(t)$. This indicates that

$\displaystyle\left|\frac{\partial}{\partial
s_{j}}\mu(f(\bm{s}),\bm{s})-\theta(t)\cdot\frac{\partial}{\partial
s_{j}}f(\bm{s})\right|=\left|v_{j}(\bm{s})-\theta(t)\cdot
g_{j}(\bm{s})\right|=\left|\frac{\partial}{\partial
s_{j}}\eta(\bm{s})\right|\leq\rho_{2}$

for $j=1,...,d$ and therefore, a nonparametric bound on $\theta(t)$ is given
by

$\theta(t)\in\bigcap_{j=1}^{d}\left[\frac{v_{j}(\bm{s})-\mathrm{sign}(g_{j}(\bm{s}))\cdot\rho_{2}}{g_{j}(\bm{s})},\;\frac{v_{j}(\bm{s})+\mathrm{sign}(g_{j}(\bm{s}))\cdot\rho_{2}}{g_{j}(\bm{s})}\right].$
(23)

This bound is well-defined and nonempty under the second inequality in
Assumption A8. More importantly, the nonparametric bound (23) sheds light on
the critical role of the covariate effect function $f$ of $\bm{S}$ on
treatment $T$ in bounding the possible value of $\theta(t)$. Specifically,
when the variation of $f$ is higher (_i.e._ , $\nabla f(\bm{s})$ is larger in
its magnitude), the nonparametric bound (23) will be tighter.

##  Appendix B Proofs

This section consists of the proofs for theorems, lemmas, and propositions in
the main paper as well as other auxiliary results.

###  B.1 Proof of Proposition 1

######  Proposition 1 (Properties of the additive confounding model (3)).

Let $\theta_{M}(t)=\mathbb{E}\left[\frac{\partial}{\partial
t}\mu(t,\bm{S})\right]$ and
$\theta_{C}(t)=\mathbb{E}\left[\frac{\partial}{\partial
t}\mu(t,\bm{S})\big{|}T=t\right]$, where
$\mu(t,\bm{s})=\mathbb{E}\left[Y|T=t,\bm{S}=\bm{s}\right]$. Under the additive
confounding model (3) with $\mbox{$\mathbb{E}$}\left[\eta(\bm{S})\right]=0$,
the following results hold for all $t\in\mathcal{T}$:

  1. (a)

$m(t)=\bar{m}(t)$.

  2. (b)

$m(t)\neq\mbox{$\mathbb{E}$}(Y|T=t)=\bar{m}(t)+\mbox{$\mathbb{E}$}\left[\eta(\bm{S})|T=t\right]$
when $\mbox{$\mathbb{E}$}\left[\eta(\bm{S})|T=t\right]\neq 0$.

  3. (c)

$\theta(t)=\theta_{M}(t)=\theta_{C}(t)$.

  4. (d)

$\theta(t)\neq\frac{d}{dt}\mbox{$\mathbb{E}$}\left[\mu(t,\bm{S})|T=t\right]=\theta(t)+\frac{d}{dt}\mbox{$\mathbb{E}$}\left[\eta(\bm{S})|T=t\right]$
when $\frac{d}{dt}\mbox{$\mathbb{E}$}\left[\eta(\bm{S})|T=t\right]\neq 0$.

  5. (e)

$\mathbb{E}\left[\mu(T,\bm{S})\right]=\mathbb{E}\left[m(T)\right]$ even when
$\mbox{$\mathbb{E}$}\left[\eta(\bm{S})\right]\neq 0$.

The above results hold even if the treatment variation $E=0$ almost surely.

###### Proof of Proposition 1.

All the results follow from some simple calculations and the fact that

$\mu(t,\bm{s})=\mbox{$\mathbb{E}$}\left(Y|T=t,\bm{S}=\bm{s}\right)=\bar{m}(t)+\eta(\bm{s})$

under the additive confounding model (3).

(a) Given that $\mbox{$\mathbb{E}$}\left[\eta(\bm{S})\right]=0$, we note that

$\displaystyle m(t)$
$\displaystyle=\mbox{$\mathbb{E}$}\left[\mu(t,\bm{S})\right]=\mbox{$\mathbb{E}$}\left[\mbox{$\mathbb{E}$}\left(Y|T=t,\bm{S}\right)\right]=\mbox{$\mathbb{E}$}\left[\bar{m}(t)+\eta(\bm{S})\right]=\bar{m}(t).$

(b) By (a), we know that

$\displaystyle\mbox{$\mathbb{E}$}\left(Y|T=t\right)=\bar{m}(t)+\mbox{$\mathbb{E}$}\left[\eta(\bm{S})|T=t\right]\neq\bar{m}(t)=m(t)$

when $\mbox{$\mathbb{E}$}\left[\eta(\bm{S})|T=t\right]\neq 0$.

(c) First, we notice that

$\displaystyle\theta_{M}(t)\equiv\mbox{$\mathbb{E}$}\left[\frac{\partial}{\partial
t}\mu(t,\bm{S})\right]=\mbox{$\mathbb{E}$}\left[\bar{m}^{\prime}(t)\right]=\bar{m}^{\prime}(t)=m^{\prime}(t)=\theta(t).$

In addition, we also have that

$\displaystyle\theta_{C}(t)$
$\displaystyle\equiv\mbox{$\mathbb{E}$}\left[\frac{\partial}{\partial
t}\mu(t,\bm{S})\Big{|}T=t\right]=\mbox{$\mathbb{E}$}\left[\bar{m}^{\prime}(t)|T=t\right]=\bar{m}^{\prime}(t)=\theta(t).$

(d) From (c), we know that

$\frac{d}{dt}\mbox{$\mathbb{E}$}\left[\mu(t,\bm{S})|T=t\right]=\frac{d}{dt}\left[\bar{m}(t)+\mbox{$\mathbb{E}$}\left(\eta(\bm{S})|T=t\right)\right]=\bar{m}^{\prime}(t)+\frac{d}{dt}\mbox{$\mathbb{E}$}\left(\eta(\bm{S})|T=t\right)\neq\theta(t)$

when $\frac{d}{dt}\mbox{$\mathbb{E}$}\left[\eta(\bm{S})|T=t\right]\neq 0$.

(e) We know that the left hand side is equal to

$\mathbb{E}\left[\mu(T,\bm{S})\right]=\mathbb{E}\left[\bar{m}(T)+\eta(\bm{S})\right]=\mathbb{E}\left[\bar{m}(T)\right]+\mathbb{E}\left[\eta(\bm{S})\right],$

while the right hand side is given by

$\mathbb{E}\left[m(T)\right]=\mathbb{E}\left\\{\bar{m}(T)+\mathbb{E}\left[\eta(\bm{S})\right]\right\\}=\mathbb{E}\left[\bar{m}(T)\right]+\mathbb{E}\left[\eta(\bm{S})\right].$

This result thus holds even when $\mathbb{E}\left[\eta(\bm{S})\right]\neq 0$.

Finally, when deriving (a)-(e), we do not impose any condition on the
treatment variation $E$, the results hold no matter what distribution of $E$
is. ∎

###  B.2 Proof of Lemma 2

Before proving Lemma 2, we give a complete statement of the asymptotic
expressions for the conditional variance and bias of
$\widehat{\beta}_{2}(t,\bm{s})$ as follows.

######  Lemma 2 (Pointwise convergence of $\widehat{\beta}_{2}(t,\bm{s})$).

Suppose that Assumptions A3, A4, A5, and A6(a) hold under the confounding
model (2). Let $\widehat{\beta}_{2}(t,\bm{s})$ be the second element of
$\widehat{\bm{\beta}}(t,\bm{s})\in\mathbb{R}^{q+1}$,
$\beta_{2}(t,\bm{s})=\frac{\partial}{\partial t}\mu(t,\bm{s})$, and
$\mathbb{X}=\left\\{(T_{i},\bm{S}_{i})\right\\}_{i=1}^{n}$. Then, for any
$(t,\bm{s})\in\mathcal{E}^{\circ}$ and $q>0$, as $h,b,\frac{b^{4}}{h}\to 0$
and $nh^{3}b^{d}\to\infty$, we know that the asymptotic conditional covariance
of $\widehat{\beta}_{2}(t,\bm{s})\in\mathbb{R}$ is given by

$\displaystyle\mathrm{Var}\left[\widehat{\beta}_{2}(t,\bm{s})\big{|}\mathbb{X}\right]$
$\displaystyle=\frac{\sigma^{2}}{nh^{3}b^{d}\cdot
p(t,\bm{s})}\left[\bm{e}_{2}^{T}\bm{M}_{q}^{-1}\bm{M}_{q}^{*}\bm{M}_{q}^{-1}\bm{e}_{2}+O\left(\max\\{h,b\\}\right)+O_{P}\left(\sqrt{\frac{1}{nhb^{d}}}\right)\right],$

and the asymptotic conditional bias is given by

$\displaystyle\mathrm{Bias}\left[\widehat{\beta}_{2}(t,\bm{s})\big{|}\mathbb{X}\right]=\begin{cases}\frac{1}{p(t,\bm{s})}\left[h^{q+1}\tau_{q}^{\text{odd}}+b^{2}\tau_{q}^{*}+O\left(\frac{b^{4}}{h}\right)+o_{P}\left(\sqrt{\frac{1}{nhb^{d}}}\right)\right]&q\text{
is odd},\\\
\frac{1}{p(t,\bm{s})}\left[h^{q}\tau_{q}^{\text{even}}+b^{2}\tau_{q}^{*}+O\left(\frac{b^{4}}{h}\right)+o_{P}\left(\sqrt{\frac{1}{nhb^{d}}}\right)\right]&q\text{
is even}.\end{cases}$

Here, $\mathcal{E}^{\circ}$ is the interior of the support