# Nonparametric estimation of a covariate-adjusted counterfactual treatment
regimen response curve

Ashkan<EMAIL_ADDRESS>[ Luke
<EMAIL_ADDRESS>[
Brent A<EMAIL_ADDRESS>[
Mark J. van der<EMAIL_ADDRESS>[ Department of Biostatistics,
University of Rochesterpresep=, ]e1,e2,e3  Department of Biostatistics,
University of California at Berkeleypresep=, ]e4

###### Abstract

Flexible estimation of the mean outcome under a treatment regimen (i.e., value
function) is the key step toward personalized medicine. We define our target
parameter as a conditional value function given a set of baseline covariates
which we refer to as a stratum based value function. We focus on
semiparametric class of decision rules and propose a sieve based nonparametric
covariate adjusted regimen-response curve estimator within that class. Our
work contributes in several ways. First, we propose an inverse probability
weighted nonparametrically efficient estimator of the smoothed regimen-
response curve function. We show that asymptotic linearity is achieved when
the nuisance functions are undersmoothed sufficiently. Asymptotic and finite
sample criteria for undersmoothing are proposed. Second, using Gaussian
process theory, we propose simultaneous confidence intervals for the smoothed
regimen-response curve function. Third, we provide consistency and convergence
rate for the optimizer of the regimen-response curve estimator; this enables
us to estimate an optimal semiparametric rule. The latter is important as the
optimizer corresponds with the optimal dynamic treatment regimen. Some finite-
sample properties are explored with simulations.

,

,

,

cadlag functions ,

dynamic marginal structural models,

highly adaptive lasso,

kernel smoothing,

semiparametric decision rules,

undersmoothing,

###### keywords:

[class=MSC]

###### keywords:

appReferences

, and

##  1 Introduction

Marginal structural models have been widely used in causal inference as a tool
for estimating the mean outcome under a give decision rule (Robins, 1998) and
constructing optimal treatment regimens. This approach requires modeling the
potential outcome given a set of baseline covariates, including treatment and
possibly some subset of baseline information on a subject. Existing methods
typically impose a parametric causal model on the mean of the potential
outcome and draw inference under the assumption that the model is correctly
specified (Robins, Hernan and Brumback, 2000; Murphy et al., 2001; Joffe et
al., 2004; Orellana, Rotnitzky and Robins, 2010). The latter is particularly
restrictive when the conditioning set of treatment and baseline covariates
includes some continuous variables. Although model selection approaches have
been proposed (Van der Laan, Laan and Robins, 2003), it is likely that finite
sample size prevents recovering the true model and leads to poor out-of-sample
performance (Shinozaki and Suzuki, 2020). Neugebauer and van der Laan (2007)
proposed a nonparametric marginal structural model framework to improve the
causal interpretability of the estimated parameters of a potentially
misspecified parametric working model. However, when the working model is
misspecified, the causally interpretable parameters might not be the
parameters of interest and the resulting treatment regimen might be
suboptimal.

The parameters of the proposed marginal structural model can be estimated
using an inverse probability weighted loss functions where weights are the
conditional probability of following a particular decision rule. The
consistency and causal interpretability of the corresponding estimators rely
on the correctly specified weight function. To mitigate the chance of model
misspecification, one may decide to model the weight function
nonparametrically using data-adaptive techniques. However, this choice may
result in nonregular estimators of the primary parameters in the marginal
structural model. This lack of regularity is because the rate of convergence
for inverse probability weighted estimators rely on a rate of convergence of
the weight function which will be slower than the desired root-$n$ rate when
data-adaptive techniques are used. To overcome this concern, authors have used
so-called doubly robust estimators for the parameters in the working marginal
structural model; that is, estimators that are consistent as long as one of
either the treatment model or the outcome models are correctly specified
(Robins, 2004; Wahed and Tsiatis, 2004; Orellana, Rotnitzky and Robins, 2010;
Rosenblum, 2011; Petersen et al., 2014; Ertefaie et al., 2015). However,
doubly-robust nonparametric estimators may suffer from two important issues:
(1) the performance of the estimators depends on the modeling choice of
nuisance parameters and they can be irregular with large bias and a slow rate
of convergence (van der Laan and Luedtke, 2014; Benkeser et al., 2017); and
(2) the quality of the constructed decision rule still relies on correctly-
specified marginal structural model.

Alternative to the marginal structural models framework, several authors have
proposed methods to estimate nonparametric decision rules that circumvent the
need for specifying a parametric outcome model. Zhao et al. (2012) proposed an
inverse probability weighting approach to estimate a nonparametric decision
rule (Dudík, Langford and Li, 2011; Zhao et al., 2015). Doubly robust
augmented inverse probability weighting procedures have also been proposed
(Rubin and van der Laan, 2012; Zhang et al., 2012; Liu et al., 2018; Athey and
Wager, 2021). Due to the nonparametric nature of the estimated decision rule,
inference for the resulting value function (i.e., mean outcome under the
estimated rule) is challenging. Zhao et al. (2012) developed Fisher
consistency, the excess risk bound, and the convergence rate of the value
function to provide a theoretical guarantee for their proposed method. Athey
and Wager (2021) derived bounds for the regret function under the estimated
nonparametric treatment strategy and showed that the bounds decay as
$n^{-1/2}$ as long as the Vapnik–Chervonenkis (VC) dimension of the class of
decision rules does not grow too fast. Other recent germane contributions on
estimating optimal treatment regimes and conditional average treatment effects
include Swaminathan and Joachims (2015), Kallus (2018), Nie and Wager (2021),
Zhao et al. (2022), and references therein.

We define a stratum-based value function as the conditional expectation of the
potential outcome under a treatment regimen given a prespecified set of
baseline covariates $\bm{V}$. The value is a function of regimens and we are
interested in understanding the functional association between the mean
outcome and the regimens within a class which we refer to as the regimen-
response curve. We consider a current status data setting with two levels of
coarsening at random (i.e., treatment assignment and censoring). We overcome
the two important shortcomings of the existing methods by: (1) allowing the
regimen-response curve to be estimated nonparametrically; and (2) estimating
the weight functions of the inverse probability weighted loss function using a
data-adaptive technique. Recently, Ertefaie, Hejazi and van der Laan (2022)
proposed a nonparametric efficient inverse probability weighted estimator for
the average causal effect where the weights are estimated using undersmoothed
highly adaptive lasso. We will generalize below the same technique to estimate
conditional average causal effects via nonparametric marginal structural
models. An important byproduct of our nonparametric regimen-response curve
estimator is a method to construct semiparametric decision rules where the
parameters of the rule are allowed to be a nonparametric function of the set
baseline covariates (i.e., ${V}$-specific optimal decision rule). The
parameters are defined as the optimizer of the regimen-response curve which
allows us to provide a rate of convergence for the corresponding functional
parameter estimators using empirical process theory. Our theoretical
contributions can be summarized as follows:

  * •

Considering a kernel smoothed regimen-response curve function as the target
parameter, we show that the proposed inverse probability weighted estimator
(1) is nonparametrically efficient when the weight functions are properly
undersmoothed; (2) convergences to a Gaussian process at root-$n$ rate; and
(3) is uniformly consistent at a derived rate.

  * •

Considering the regimen-response curve function as the target parameter, we
derive the asymptotic normality of our estimator and show that the bias caused
by the kernel smoothing of the estimator vanishes under some mild assumptions
and properly undersmoothed kernel and weight functions.

  * •

Using a novel empirical process algorithm, we show that the minimizer of the
estimated regimen-response curve is a consistent estimator of the minimizer of
the true regimen-response curve and derive the convergence rate. This result
is of independent interest as it derives a rate of convergence for an
optimizer of a nonparametrically estimated function.

In contrast to the existing literature which targets the value function under
the estimated rule, our estimand is the stratum based value function. Deriving
the asymptotic properties of our estimand is more challenging because it will
be nonpathwise differentiable functional when $\bm{V}$ includes at least one
continuous variable (Bickel et al., 1998). Moreover, existing methods can
either estimate a parametric or nonparametric rules where in each case
flexibility or interpretability is sacrificed. Our proposed semiparametric
rule lands in between and provides a reasonable trade-off between flexibility
and interpretability. For example, clinicians might have strong opinion about
an important tailoring variable and want to ensure to capture the correct
functional form of that variable in the decision rule. In such settings, we
can include that variable in the set $\bm{V}$ and provide a stratum-based
decision rule where the decision rule will be a nonparametric function of
$\bm{V}$.

##  2 Preliminaries

###  2.1 Notation and Formulation

Let $T$ be a univariate failure time,
$\bm{W}\in\mathcal{W}\subset\mathbb{R}^{p}$ be a vector of all available
baseline covariates measured before treatment $A\in\\{0,1\\}$, and $C$ be a
censoring variable. Let $\bm{S}\in\mathcal{S}\subset\mathbb{R}^{q}$ be a
subvector of $\bm{W}$ and $\bm{V}=\bm{W}\setminus\bm{S}$ where
$\bm{V}\in\mathcal{V}\subset\mathbb{R}^{p-q}$. For simplicity of notation, we
denote $p-q$ as $r$. A decision rule $d^{\theta}\in\mathcal{D}$ is defined as
a function that maps values of $\bm{S}$ to a treatment option $\\{0,1\\}$. In
our formulation $\mathcal{D}$ is a class of deterministic rules indexed by a
vector of coefficients. Specifically, we consider
$d^{\theta}=I\\{\bm{S}^{\top}\theta>0\\}$ with
$\theta\in\Theta\subset\mathbb{R}^{q}$ where $I(.)$ is an indicator function.
Let $T^{\theta}$ be the potential outcome that would have been observed under
the decision rule $d^{\theta}$ for a given value of $\theta\in\Theta$. Suppose
that we have the full data
$X=\\{(T^{\theta}:d^{\theta}\in\mathcal{D}),\bm{W}\\}\sim
P_{X}\in\mathcal{M}^{F}$ where $\mathcal{M}^{F}$ is a nonparametric full data
model. Define
$E_{P_{X}}\\{I(T^{\theta}>t)\mid\bm{V}\\}=\Psi_{P_{X}}(\theta,\bm{V})\in\mathcal{F}_{\lambda}$
where $\mathcal{F}_{\lambda}$ is a class of cadlag functions with sectional
variation norm bounded by $\lambda$ and $t$ is a predetermined time point. The
full-data risk of $\Psi$ is defined as
$E_{P_{X}}\int_{\theta}\ell_{\theta}(\Psi,x)dF(\theta)$ where $\ell$ is the
log-likelihood loss function,
$\ell_{\theta}(\Psi)=I(T^{\theta}>t)\log\left(\frac{\Psi}{1-\Psi}\right)+\log(1-\Psi),$
and $dF(\theta)$ is a weight function/measure on the Euclidean set of all
$\theta$-values (e.g., multivariate normal with possibly large variance).
Define a V-specific optimal decision rule that maximizes the $t$-year survival
probability as
$\theta_{0}(\bm{V})=\operatorname*{argmax}_{\theta\in\Theta}\Psi_{P_{X}}(\theta,\bm{V})$.
Hence, we allow $\theta$ to be a functional parameter mapping $\mathcal{V}$
onto $\mathbb{R}^{q}$.

The observed data $\mathcal{O}=\\{Y=\min(T,C),\Delta=I(T<C),A,\bm{W}\\}$
follow some probability distribution $P_{0}$ that lies in some nonparametric
model $\mathcal{M}$. Suppose we observe $n$ independent, identically
distributed trajectories of $\mathcal{O}$. Let
$\Delta^{c}=\Delta+(1-\Delta)I\\{Y>t\\}$ and
$G:\mathcal{M}\rightarrow\mathcal{G}$ be a functional nuisance parameter where
$\mathcal{G}=\\{G(P):P\in\mathcal{M}\\}$. In our setting,
$G(P)\\{(A,\Delta^{c})\mid\bm{W}\\}$ denotes the joint treatment and censoring
mechanism under a distribution $P$. We refer to $G(P_{0})$ as $G_{0}$ and
$G(P)$ as $G$. Under coarsening at random assumption,
$G(P)\\{(A,\Delta^{c})\mid{X}\\}=G(P)\\{(A,\Delta^{c})\mid\bm{W}\\}$. Denote
$G^{a}\equiv G^{a}(P)(A\mid\bm{W})$ and $G^{c}\equiv
G^{c}(P)(\Delta^{c}\mid\bm{W},A)$. We denote the latter two functions under
$P_{0}$ as $G^{a}_{0}$ and $G^{c}_{0}$. Also define
$Q:\mathcal{M}\rightarrow\mathcal{Q}$ and
$\mathcal{Q}=\\{Q(P):P\in\mathcal{M}\\}$ where
$Q_{0}(\theta,\bm{W})=\mathbb{E}_{P_{0}}\\{I(Y^{\theta}>t)\mid\bm{W}\\}$. Let
$\prod$ be a projection operator in the Hilbert space $L_{0}^{2}(P)$ with
inner product $\langle h_{1},h_{2}\rangle=E_{P}(h_{1}h_{2})$.

The loss function $\int_{\theta}\ell_{\theta}(\Psi)dF(\theta)$ is based on the
full-data which is not observable, and thus, cannot be used to define our full
data functional parameter $\Psi_{P_{X}}$. However, under consistency
($T_{i}=T^{A_{i}}$) and no unmeasured confounder ($A\perp T^{a}\mid\bm{W}$ and
$\Delta^{c}\perp T\mid A,\bm{W}$) assumptions, it can be identified using the
observed data distribution $P_{0}$. Specifically, we define an inverse
probability weighting mapping of the full-data loss as

$\displaystyle
L_{G}(\Psi)=\int_{\theta}\left[\frac{\Delta^{c}I(A=d^{\theta})\xi(\theta,\bm{V})}{G\\{(A,\Delta^{c})\mid\bm{W}\\}}\\{I(Y>t)\log\left(\frac{\Psi}{1-\Psi}\right)+\log(1-\Psi)\\}\right]dF(\theta).$
(1)

Accordingly, define the regimen-response curve as
$\Psi_{P_{0}}=\operatorname*{argmin}_{\Psi\in\mathcal{F}_{\lambda}}E_{P_{0}}L_{G_{0}}(\Psi)$.
Following our convention, we denote $\Psi_{P_{0}}$ as $\Psi_{0}$. We define a
highly adaptive estimator of $\Psi_{P_{0}}$ as
$\Psi_{n}=\operatorname*{argmin}_{\Psi\in\mathcal{F}_{\lambda_{n}}}P_{n}L_{G_{n}}(\Psi)$
where $\lambda_{n}$ is a data adaptively selected $\lambda$ and $G_{n}$ is an
estimator of $G_{0}$. The choice of function $\xi$ is inconsequential to the
derivation of the efficient influence function of our target parameter but it
may impact the level of undersmoothing needed in finite sample. One may
consider $\xi(\theta,\bm{V})=1$ or
$\xi(\theta,\bm{V})=G\\{(A,\Delta^{c})\mid\bm{V}\\}$. The latter choice
corresponds to the stabilized weight.

###  2.2 Highly adaptive lasso

The highly adaptive lasso is a nonparametric regression approach that can be
used to estimate infinite dimensional functional parameters Benkeser and Van
Der Laan (2016); van der Laan (2017); van der Laan and Bibaut (2017). It forms
a linear combination of indicator basis functions to minimize the expected
value of a loss function under the constraint that the constraint that the
$L_{1}$-norm of the coefficient vector is bounded by a finite constant value.

Let $\mathcal{D}[0,\tau]$ be the Banach space of d-variate real valued cadlag
functions (right-continuous with left-hand limits) on a cube
$[0,\tau]\in\mathbb{R}^{d}$. For each function $f\in\mathcal{D}[0,\tau]$
define the supremum norm as $\|f\|_{\infty}=\sup_{w\in[0,\tau]}|f(w)|$. For
any subset $s$ of $\\{0,1,\ldots,d\\}$, we partition $[0,\tau]$ in
$\\{0\\}\\{\cup_{s}(0_{s},\tau_{s}]\\}$ and define the sectional variation
norm of such $f$ as

$\|f\|^{*}_{\zeta}=|f(0)|+\sum_{s\subset\\{1,\ldots,d\\}}\int_{0_{s}}^{\tau_{s}}|df_{s}(u_{s})|,$

where the sum is over all subsets of $\\{0,1,\ldots,d\\}$. For a given subset
$s\subset\\{0,1,\ldots,d\\}$, define $u_{s}=(u_{j}:j\in s)$ and $u_{-s}$ as
the complement of $u_{s}$. Then, $f_{s}:[0_{s},\tau_{s}]\rightarrow\mathbb{R}$
defined as $f_{s}(u_{s})=f(u_{s},0_{-s})$. Thus, $f_{s}(u_{s})$ is a section
of $f$ that sets the components in the complement of $s$ equal to zero and
varies only along the variables in $u_{r}$.

Let $G^{a}\equiv G(P)(A\mid W)$ and $G^{c}\equiv G(P)(\Delta^{c}\mid W,A)$
denote the treatment and censoring mechanism under an arbitrary distribution
$P\in\mathcal{M}$. Assuming that our nuisance functional parameter
$G^{a},G^{c}\in\mathcal{D}[0,\tau]$ has finite sectional variation norm, we
can represent $\text{logit }G^{a}$ as (Gill, van der Laan and Wellner, 1993)

$\displaystyle\text{logit }G^{a}(w):$ $\displaystyle=\text{logit
}G^{a}(0)+\sum_{s\subset\\{1,\ldots,d\\}}\int_{0_{s}}^{w_{s}}dG^{a}_{s}(u_{s})$
$\displaystyle=\text{logit
}G^{a}(0)+\sum_{s\subset\\{1,\ldots,d\\}}\int_{0_{s}}^{\tau_{s}}I(u_{s}\leq
w_{s})dG^{a}_{s}(u_{s}).$ (2)

The representation (2) can be approximated using a discrete measure that puts
mass on each observed $W_{s,i}$ denoted by $\alpha_{s,i}$. Let
$\phi_{s,i}(c_{s})=I(w_{i,s}\leq c_{s})$ where $w_{i,s}$ are support points of
$dG^{a}_{s}$. We then have

$\text{logit
}G^{a}_{\alpha}=\alpha_{0}+\sum_{s\subset\\{1,\ldots,d\\}}\sum_{i=1}^{n}\alpha_{s,i}\phi_{s,i},$

where
$|\alpha_{0}|+\sum_{s\subset\\{1,\ldots,d\\}}\sum_{i=1}^{n}|\alpha_{s,i}|$ is
an approximation of the sectional variation norm of $f$. The loss based highly
adaptive lasso estimator $\beta_{n}$ is defined as

$\alpha_{n}(\lambda)=\arg\min_{\alpha:|\alpha_{0}|+\sum_{s\subset\\{1,\ldots,d\\}}\sum_{i=1}^{n}|\alpha_{s,i}|<\lambda}P_{n}L(\text{logit
}G^{a}_{\alpha}).$

where $L(.)$ is a given loss function. Denote $G^{a}_{n,\lambda}\equiv
G^{a}_{\alpha_{n}(\lambda)}$ as the highly adaptive lasso estimate of $G^{a}$.
Similarly, we can define $G^{c}_{n,\lambda}$ as the highly adaptive lasso
estimate $G^{c}$. When the nuisance parameter correspond to a conditional
probability of a binary variable (e.g., propensity score with a binary
treatment or a binary censoring indicator), log-likelihood loss can be used.
In practice, $\lambda$ is unknown and must be obtained using data. We refer to
$\lambda_{n}$ as a value of $\lambda$ that is determined using data.

##  3 Dynamic marginal structural models with highly adaptive lasso

###  3.1 The challenge

The target functional parameter $\Psi_{0}$ is not a pathwise differentiable
function which makes the inference challenging. The existing literature
approximates this functional by pathwise differentiable function
$\Psi_{\beta}$, where $\beta$ is a vector of unknown parameters van der Laan
(2006); van der Laan and Petersen (2007); Murphy et al. (2001); Robins,
Orellana and Rotnitzky (2008); Orellana, Rotnitzky and Robins (2010);
Neugebauer and van der Laan (2007). The working model provides a user-
specified summary of the unknown true function relating the expectation of the
potential outcome under a strategy. Hence, the quality of the constructed
decision rule depends on how well the working model approximates the true
function. To mitigate this concern we propose a kernel smoothed highly
adaptive lasso estimator of $\Psi_{0}(\bm{v}_{0},\theta)$. Specifically, we
define an regimen-response curve estimator as

$\displaystyle\Psi_{nh}(\bm{v}_{0},\theta)=P_{n}U(P_{n},\bm{v}_{0},\Psi_{n})$
(3)

where
$U(P_{n},\bm{v}_{0},\Psi_{n})=(P_{n}K_{h,\bm{v}_{0}})^{-1}\\{K_{h,\bm{v}_{0}}\Psi_{n}(\theta)\\}$
and $\Psi_{n}$ is a highly adaptive lasso estimate of $\Psi_{0}$ obtained
based on the loss function (1). In this formulation,
$K_{h,\bm{v}_{0}}(v)=h^{-r}K\\{(v-\bm{v}_{0})h^{-1}\\}$ with
$K(v)=K(v_{1},\cdots,v_{r})$ where $K(v)$ is a kernel satisfying $\int
K(v)dv=1$. Moreover, we assume that $K(v)=\prod_{j=1}^{r}K_{u}(v_{j})$ is a
product kernel defined by a univariate kernel $K_{u}(.)$ and $K(v)$ is a
J-orthogonal kernel such that $\int K_{u}(v)dv=1$ and $\int K_{u}(v)v^{j}dv=0$
for $j=1,\cdots,J$.

We study the asymptotic behaviour of our estimator relative to
$\Psi_{0}(\bm{v_{0}},\theta)$ and a kernel smoothed parameter

$\Psi_{0h}(\bm{v}_{0},\theta)=P_{0}U(P_{0},\bm{v}_{0},\Psi_{0}),$

where $h$ is a fixed bandwidth. In general, when the nuisance parameter
$G_{0}$ is estimated using data adaptive regression techniques, the resulting
regimen-response curve estimator $\Psi_{nh}(\bm{v}_{0},\theta)$ won’t be
asymptotically linear which makes inference challenging. For example, consider
the target parameter $\Psi_{0h}(\bm{v}_{0},\theta)$, we have

$\displaystyle\Psi_{nh}(\bm{v}_{0},\theta)-\Psi_{0h}(\bm{v}_{0},\theta)=$
$\displaystyle(P_{n}-P_{0})U(P_{0},\bm{v}_{0},\Psi_{0})+$ $\displaystyle
P_{0}\\{U(P_{n},\bm{v}_{0},\Psi_{n})-U(P_{0},\bm{v}_{0},\Psi_{0})\\}+$
$\displaystyle(P_{n}-P_{0})\\{U(P_{n},\bm{v}_{0},\Psi_{n})-U(P_{0},\bm{v}_{0},\Psi_{0})\\}.$
(4)

The first term on the right-hand side of (4) is exactly linear and the third
term is negligible under Donsker condition of $\Psi_{n}$. The problem arises
because data adaptive regression techniques have a rate of convergence slower
than the desired root-$n$ rate, and thus, the bias of
$P_{0}\\{U(P_{n},\bm{v}_{0},\Psi_{n})-U(P_{0},\bm{v}_{0},\Psi_{0})\\}$
diverges to infinity.

###  3.2 The intuition

We show that when the nuisance parameters are properly undersmoothed, the
asymptotic linearity of our estimator $\Psi_{nh}(\bm{v}_{0},\theta)$ can be
retrieved. Specifically, in the proof of Theorem 4.2, we show that

$\displaystyle
P_{0}\\{U(P_{n},\bm{v}_{0},\Psi_{n})-U(P_{0},\bm{v}_{0},\Psi_{0})\\}=$
$\displaystyle(P_{n}-P_{0})\tilde{D}^{*}_{h,\bm{v_{0}}}(P_{0})+o_{p}(n^{-1/2})$
$\displaystyle\hskip
14.45377pt+P_{n}\frac{K_{h,\bm{v}_{0}}\Delta^{c}I(A=d^{\theta})}{G_{n}P_{n}K_{h,\bm{v}_{0}}}\left\\{I(Y\geq
t)-\Psi_{n}\right\\}$ (5) $\displaystyle\hskip
14.45377pt+P_{n}\frac{K_{h,\bm{v}_{0}}I(A=d^{\theta})}{G_{n}P_{n}K_{h,\bm{v}_{0}}}(Q_{n}-\Psi_{n})\left(\Delta^{c}-G_{n}^{c}\right)$
(6) $\displaystyle\hskip
14.45377pt+P_{n}\frac{K_{h,\bm{v}_{0}}(Q_{n}-\Psi_{n})}{G_{n}^{a}P_{n}K_{h,\bm{v}_{0}}}\left\\{I(A=d^{\theta})-G_{n}^{a}\right\\},$
(7)

where $\tilde{D}^{*}_{h,\bm{v_{0}}}(P_{0})$ includes certain components of the
canonical gradient of our statistical parameter $\Psi_{0h}(\bm{v_{0}},\theta)$
and thus, contributes to the efficient influence function of our estimator.
The last three terms do not converge to zero at the appropriate rate and thus
induce non-negligible bias. However, the terms (5), (6) and (7) resemble the
form of score equations for $I(Y\geq t)$, $\Delta^{c}$ and $I(A=d^{\theta})$.

We first discuss the score equations for $G^{.}$ where $.\in\\{a,c\\}$. We can
generate score $S_{g}(G_{n}^{.})$ for a nonparametric estimator $G_{n}^{.}$ by
paths $\\{G_{n,\epsilon}^{.g}(w)\\}$, such that

$\displaystyle\text{logit }G_{n,\epsilon}^{.g}(w)$ $\displaystyle=\text{logit
}G_{n}^{.}(0)\\{1+\epsilon
g(0)\\}+\sum_{s\subset\\{1,\ldots,d\\}}\int_{0_{s}}^{w_{s}}\phi_{s,u_{s}}(w)\\{1+\epsilon
g(s,u_{s})\\}d\text{logit }G_{n,s}^{.}(u),$

where $g$ is a uniformly bounded function on $[0,\tau]^{p}$. Accordingly when
a function is estimated using a highly adaptive lasso, the score functions can
be generated by a path $\\{1+\epsilon g(s,j)\\}\beta_{n,s,j}$ for a bounded
vector $g$ as

$S_{g}(G_{n,\lambda_{n}}^{.})=\frac{d}{d\text{logit
}G^{.}_{n,\lambda_{n}}}L(\text{logit
}G^{.}_{n,\lambda_{n}})\left\\{\sum_{(s,j)}g(s,j)\beta_{n,s,j}\phi_{s,j}\right\\},$
(8)

where $L(\cdot)$ is the log-likelihood loss function. For example, for the
treatment indicator,
$L(G^{a})=A\log\left(\frac{G^{a}}{1-G^{a}}\right)+\log(1-G^{a})$. The penalty
in the highly adaptive lasso imposes a restriction on the $g$ function.
Specifically, the path $\\{1+\epsilon g(s,j)\\}\beta_{n,s,j}$ can be generated
using those $g$ functions that do not change the $L_{1}$-norm of the
parameters. Let $\beta^{g}_{n,s,j}=\\{1+\epsilon g(s,j)\\}\beta_{n,s,j}$
denote the set of perturbed parameters. Then, for small enough $\epsilon$ such
that $\\{1+\epsilon g(s,j)\\}>0$,
$\|\beta^{g}_{n,s,j}\|_{1}=\|\beta^{g}_{n,s,j}\|_{1}\\{1+\epsilon
g(s,j)\\}=\|\beta^{g}_{n,s,j}\|_{1}+\epsilon g(s,j)\|\beta_{n,s,j}\|_{1}$.
Hence, the restriction is satisfied when the inner product of $g$ and the
vector $|\beta|$ is zero (i.e., $<g.|\beta|>=0$).

We now provide a thought experiment on how undersmoothing the nuisance
function estimates using a highly adaptive lasso can eliminate the terms (5),
(6) and (7). We first ignore the $L_{1}$-norm restriction on the choice of
function $g$ which is erroneous but contains the main idea of the approach.
Without the restriction, one can choose $g(s_{0},j_{0})=I(s=s_{0},j=j_{0})$.
The latter choice perturbs one coefficient at a time and corresponds to
maximum likelihood estimators. Then, for any $s_{0}$ and $j_{0}$, the score
function (19) becomes

$S_{I}(G_{n,\lambda_{n}}^{.})=\frac{d}{d\text{logit
}G^{.}_{n,\lambda_{n}}}L(\text{logit
}G^{.}_{n,\lambda_{n}})\left(\beta_{n,s_{0},j_{0}}\phi_{s_{0},j_{0}}\right).$

Because, $\beta_{n,s_{0},j_{0}}$ is a finite constant for all $s_{0}$ and
$j_{0}$, solving the above score equation is equivalent to solving

$\frac{d}{d\text{logit }G^{.}_{n,\lambda_{n}}}L(\text{logit
}G^{.}_{n,\lambda_{n}})\phi_{s_{0},j_{0}}.$

For example, for the treatment indicator, the score equation is given by
$(A-G_{n,\lambda_{n}}^{a})\phi_{s_{0},j_{0}}$. As we undersmooth the fit, we
solve more and more score equations (i.e., the number of score equations
increases with the number of features included in the model) and any linear
combination of those score equations will also be solved. Now let’s consider
the term (7) and let
$f=\frac{K_{h,\bm{v}_{0}}(Q_{n}-\Psi_{n})}{G_{n}^{a}P_{n}K_{h,\bm{v}_{0}}}$.
Assuming that $Q_{0}$, $\Psi_{0}$ and $G_{0}=G_{0}^{a}G_{0}^{c}$ are càdlàg
with finite sectional variation norm (Assumption 1 in Section 4), Gill, van
der Laan and Wellner (1993) showed that $f$ can be approximated as a linear
combination of indicator basis functions. Therefore, if we sufficiently
undersmooth $G_{n,\lambda}^{a}$ we will solve (7) up to a root-$n$ factor.
That is

$P_{n}\frac{K_{h,\bm{v}_{0}}(Q_{n}-\Psi_{n})}{G_{n}^{a}P_{n}K_{h,\bm{v}_{0}}}\left\\{I(A=d^{\theta})-G_{n}^{a}\right\\}=o_{p}(n^{-1/2}).$

The same argument can be made to show that the other terms (5) and (6) can be
made negligible when $\Psi_{n}$ and $G_{n}^{c}$ are properly undersmoothed.
Hence, the challenging term will be asymptotically linear with

$\displaystyle
P_{0}\\{U(P_{n},\bm{v}_{0},\Psi_{n})-U(P_{0},\bm{v}_{0},\Psi_{0})\\}=$
$\displaystyle(P_{n}-P_{0})\tilde{D}^{*}_{h,\bm{v_{0}}}(P_{0})+o_{p}(n^{-1/2}).$

Now, we study the actual case where the choice function $h$ is restricted to
those that do not change the $L_{1}-$norm of the coefficients. Lemma A.2 shows
than when (23), (24) and (25) are satisfied, our proposed estimator achieves
asymptotic linearity. As we undersmooth the fit, we start adding features with
small coefficients into the model. Conditions (23), (24) and (25) imply that
the $L_{1}-$norm must be increases until one of the corresponding score
equations is solved to a precision of $o_{p}(n^{-1/2})$.

Here we provide details for the score equations corresponding to the treatment
indicator where
$S_{g}(G_{n,\lambda_{n}}^{a})=(A-G_{n,\lambda_{n}}^{a})\left\\{\sum_{(s,j)}g(s,j)\beta_{n,s,j}\phi_{s,j}\right\\}$.
Let $r(g,G_{n,\lambda_{n}}^{a})=\sum_{(s,j)}g(s,j)\lvert\beta_{n,s,j}\rvert$.
For small enough $\epsilon$,

$\displaystyle\sum_{(s,j)}\lvert\\{1+\epsilon g(s,j)\\}\beta_{n,s,j}\rvert$
$\displaystyle=\sum_{(s,j)}\\{1+\epsilon g(s,j)\\}\lvert\beta_{n,s,j}\rvert$
$\displaystyle=\sum_{(s,j)}\lvert\beta_{n,s,j}\rvert+\epsilon
r(g,G_{n,\lambda_{n}}^{a}).$

Hence, for any $g$ satisfying $r(g,G_{n,\lambda_{n}}^{a})=0$ (i.e., $h$ does
not change the $L_{1}-$ norm of the coefficients), we have
$P_{n}S_{g}(G_{n,\lambda_{n}}^{a})=0$.

Let $D\\{f,G_{n,\lambda_{n}}^{a}\\}=f\cdot(A-G_{n,\lambda_{n}}^{a})$, where
$f$ is defined above, and let $\tilde{f}$ be an approximation of $f$ using the
basis functions that satisfy condition (24). Then,
$D\\{\tilde{f},G_{n,\lambda_{n}}^{a}\\}\in\\{S_{g}(G_{n,\lambda_{n}}^{a}):\lVert
g\rVert_{\infty}<\infty\\}$. Thus, there exists $g^{\star}$ such that
$D(\tilde{f},G_{n,\lambda_{n}}^{a})=S_{g^{\star}}(G_{n,\lambda_{n}}^{a})$;
however, for this particular choice of $g^{\star}$,
$r(g,G_{n,\lambda_{n}}^{a})$ may not be zero (i.e., the restriction on $g$
might be violated). Now, define $g$ such that $g(s,j)=g^{\star}(s,j)$ for
$(s,j)\neq(s^{\star},j^{\star})$; $\tilde{g}(s^{\star},j^{\star})$ is defined
such that

$\displaystyle\sum_{(s,j)\neq(s^{\star},j^{\star})}g^{\star}(s,j)\lvert\beta_{n,s,j}\rvert+g(s^{\star},j^{\star})\lvert\beta_{n,s^{\star},j^{\star}}\rvert=0.$
(9)

That is, $g$ matches $g^{\star}$ everywhere but for a single point
$(s^{\star},j^{\star})$, where it is forced to take a value such that
$r(g,G_{n,\lambda_{n}}^{a})=0$. As a result, for such a choice of $g$,
$P_{n}S_{h}(G_{n,\lambda_{n}}^{a})=0$ by definition. Below, we show that
$P_{n}S_{g}(G^{a}_{n,\lambda_{n}})-P_{n}D(\tilde{f},G_{n,\lambda_{n}}^{a})=o_{p}(n^{-1/2})$
which then implies that
$P_{n}D(\tilde{f},G_{n,\lambda_{n}}^{a})=o_{p}(n^{-1/2})$. We note that the
choice of $(s^{\star},j^{\star})$ is inconsequential.

$\displaystyle
P_{n}S_{g}(G_{n,\lambda_{n}}^{a})-P_{n}D\\{\tilde{f},G_{n,\lambda_{n}}^{a}\\}$
$\displaystyle=P_{n}S_{g}(G_{n,\lambda_{n}}^{a})-P_{n}S_{g^{*}}(G_{n,\lambda_{n}}^{a})$
$\displaystyle=P_{n}\left\\{\frac{d}{d\text{logit
}G^{a}_{n,\lambda_{n}}}L(\text{logit
}G^{a}_{n,\lambda_{n}})\left[\sum_{(s,j)}\left\\{g(s,j)-g^{\star}(s,j)\right\\}\beta_{n,s,j}\phi_{s,j}\right]\right\\}$
$\displaystyle=P_{n}\left[\frac{d}{d\text{logit
}G^{a}_{n,\lambda_{n}}}L(\text{logit
}G^{a}_{n,\lambda_{n}})\left\\{g(s^{\star},j^{\star})-g^{\star}(s^{\star},j^{\star})\right\\}\beta_{n,s^{\star},j^{\star}}\phi_{s^{\star},j^{\star}}\right]$
$\displaystyle=O_{p}\left(P_{n}\left[\frac{d}{d\text{logit
}G^{a}_{n,\lambda_{n}}}L(\text{logit
}G^{a}_{n,\lambda_{n}})\phi_{s^{\star},j^{\star}}\right]\right)$
$\displaystyle=o_{p}(n^{-1/2}).$

The details of the forth equality is given in the proof of Lemma A.2 and the
last equality follows from the assumption that
$\min_{(s,j)\in\mathcal{J}_{n}}\lVert P_{n}\frac{d}{d\text{logit
}G^{a}_{n,\lambda_{n}}}L(\text{logit
}G^{a}_{n,\lambda_{n}})(\phi_{s,j})\rVert=o_{p}(n^{-1/2})$ for $L(\cdot)$
being log-likelihood loss (i.e., condition (24)). As
$P_{n}S_{g}(G^{a}_{n,\lambda_{n}})=0$, it follows that
$P_{n}D(\tilde{f},G^{a}_{n,\lambda_{n}})=o_{p}(n^{-1/2})$. Using this result,
under mild assumptions, we showed that
$P_{n}D(f,G^{a}_{n,\lambda_{n}})=o_{p}(n^{-1/2})$ indicating that the term (7)
will be asymptotically negligible (see the proof of Lemma A.2 for details).

###  3.3 The estimator.

To improve the finite sample performance of our estimator we propose to
estimate the nuisance parameters using cross-fitting (Klaassen, 1987; Zheng
and van der Laan, 2011; Chernozhukov et al., 2017). We split the data at
random into $B$ mutually exclusive and exhaustive sets of size approximately
$nB^{-1}$. Let $P_{n,b}^{0}$ and $P_{n,b}^{1}$ denote the the empirical
distribution of a training and a validation sample, respectively. For a given
$\lambda$ and $h$, exclude a single (validation) fold of data and fit the
highly adaptive lasso estimator using data from the remaining $(B-1)$ folds;
use this model to estimate the nuisance parameters for samples in the holdout
(validation) fold. By repeating the process $B$ times, we will have estimates
of the nuisance parameters for all sample units. Accordingly, we define the
cross-fitted IPW estimator

$\displaystyle\Psi_{nh}^{\textsuperscript{CF}}(\bm{v}_{0},\theta)=B^{-1}\sum_{b=1}^{B}P_{n,b}^{1}U(P_{n,b}^{1},\bm{v}_{0},\Psi_{n,\lambda,b})$
(10)

where
$U(P_{n,b}^{1},\bm{v}_{0},\Psi_{n,\lambda,b})=(G_{n,\lambda,b}P_{n,b}^{1}K_{h,\bm{v}_{0}})^{-1}\\{K_{h,\bm{v}_{0}}\Delta^{c}I(A=d^{\theta})\Psi_{n,\lambda,b}\\}$,
$\Psi_{n,\lambda,b}$ and $G_{n,\lambda,b}$ are highly adaptive lasso estimates
of $\Psi_{0}$ and $G_{0}$ applied to the training sample for the bth sample
split for a given $\lambda$. The estimator uses a ($J-1$)-orthogonal kernel
with bandwidth $h$ centered at $\bm{v_{0}}$.

###  3.4 Data-adaptive bandwidth selector

The proposed estimator $\Psi_{nh}(\bm{v}_{0},\theta)$ approaches to
$\Psi_{0}(\bm{v}_{0},\theta)$ as $h_{n}\rightarrow 0$ at the cost of
increasing the variance of the estimator. Mean square error (MSE) provides an
ideal criterion for bias-variance trade-off. However, because the bias is
unknown the latter cannot be used directly. Here we propose an approach that
circumvent the need for knowing bias while targeting the optimal bias-variance
trade-off. Let $\sigma_{nh}$ be a consistent estimator of
$\sigma_{0}=[E\\{D^{*2}_{\bm{v}_{0}}(P_{0})\\}]^{1/2}$ where
$D^{*2}_{\bm{v}_{0}}(P_{0})$ is the efficient influence function for
$\Psi_{nh}(\bm{v}_{0},\theta)$ (see Theorem 4.2). The goal is to find an $h$
that minimizes the MSE or equivalently set the derivative of the MSE to zero.
That is

$\\{\Psi_{nh}(\bm{v}_{0},\theta)-\Psi_{0}(\bm{v}_{0},\theta)\\}\frac{\partial}{\partial
h}\Psi_{nh}(\bm{v}_{0},\theta)+\frac{\sigma_{nh}}{(nh^{r})^{1/2}}\frac{\partial}{\partial
h}\sigma_{nh}/(nh^{r})^{1/2}=0$

We know that the optimal (up to a constant) bias-variance trade-off implies
$\\{\Psi_{nh}(\bm{v}_{0},\theta)-\Psi_{0}(\bm{v}_{0},\theta)\\}\approx\sigma_{nh}/(nh^{r})^{1/2}$
(van der Laan, Bibaut and Luedtke, 2018). Hence, the optimal $h$ will also
solve one of the following equations

$\displaystyle\frac{\partial}{\partial
h}\Psi_{nh}(\bm{v}_{0},\theta)\pm\kappa\left(\frac{\partial}{\partial
h}\sigma_{nh}/(nh^{r})^{1/2}\right)=0,$ (11)

where $\kappa$ is a positive constant. Consider a scenario where
$\Psi_{nh}(\bm{v}_{0},\theta)$ shows an increasing trend as $h_{n}$ approaches
zero. This implies that
$\Psi_{0}(\bm{v}_{0},\theta)>\Psi_{nh}(\bm{v}_{0},\theta)-\kappa\sigma_{nh}/(nh^{r})^{1/2}$.
We then define our optimal finite sample bandwidth selector as

$\displaystyle
h_{n}=\operatorname*{argmax}_{h}\Psi_{nh}(\bm{v}_{0},\theta)-\kappa\sigma_{nh}/(nh^{r})^{1/2}.$

Similarly, when $\Psi_{nh}(\bm{v}_{0},\theta)$ shows decreasing trend as
$h_{n}$ approaches zero, we define
$h_{n}=\operatorname*{argmin}_{h}\Psi_{nh}(\bm{v}_{0},\theta)+\kappa\sigma_{nh}/(nh^{r})^{1/2}$.
A natural choice for the constant $\kappa$ is $(1-\alpha)$-quantile of the
standard normal distribution (i.e., $\zeta_{1-\alpha}$). Hence, $h_{n}$ is the
bandwidth value that minimizes the difference between the true value
$\Psi_{0}(\bm{v}_{0},\theta)$ and the corresponding lower (upper) confidence
bound.

###  3.5 Undersmoothing criteria

The asymptotic linearity results of Theorems 4.2 and 4.3 rely on estimating
the corresponding nuisance parameters using an undersmoothed highly adaptive
lasso. Specifically, Theorem 4.2 requires that

  * (a)

$\left|P_{n}D_{CAR}^{a}(P_{n},G_{n}^{a},Q_{n},\Psi_{n},\theta)\right|=o_{p}\left((nh^{r})^{-1/2}\right)$,

  * (b)

$\left|P_{n}D_{CAR}^{c}(P_{n},G_{n},Q_{n},\Psi_{n},\theta)\right|=o_{p}\left((nh^{r})^{-1/2}\right)$,

  * (c)

$\left|P_{n}D_{CAR}^{\Psi}(P_{n},G_{n},\Psi_{n},\theta)\right|=o_{p}\left((nh^{r})^{-1/2}\right)$,

where
$D_{CAR}^{a}(P_{n},G_{n}^{a},Q_{n},\Psi_{n},\theta)=\frac{K_{h,\bm{v}_{0}}}{P_{n}K_{h,\bm{v}_{0}}}(Q_{n}-\Psi_{n})\left(\frac{G_{n}^{a}-I(A=d^{\theta})}{G_{n}}\right)$,
$D_{CAR}^{c}(P_{n},G_{n},Q_{n},\Psi_{n},\theta)=\frac{K_{h,\bm{v}_{0}}I(A=d^{\theta})}{P_{n}K_{h,\bm{v}_{0}}}(Q_{n}-\Psi_{n})\left(\frac{G_{n}^{c}-\Delta^{c}}{G_{n}}\right)$
and
$D_{CAR}^{\Psi}(P_{n},G_{n},\Psi_{n},\theta)=\frac{K_{h,v}\Delta^{c}I(A=d^{\theta})}{G_{n}P_{n}K_{h,v}}\\{\Psi_{n}-I(Y>t)\\}$.
Motivated by our theoretical results, we propose the following practical
L1-norm bound selection criteria for $G_{n}^{a}$, $G_{n}^{c}$ and $\Psi_{n}$:

$\displaystyle\lambda_{n,G^{a}}$
$\displaystyle=\operatorname*{argmin}_{\lambda}\left|B^{-1}\sum_{b=1}^{B}P_{n,b}^{1}D_{CAR}^{a}(P_{n,b}^{1},G_{n,\lambda,b}^{a},G_{n,b}^{c},Q_{n,b},\Psi_{n,b},\theta)\right|,$
(12) $\displaystyle\lambda_{n,G^{c}}$
$\displaystyle=\operatorname*{argmin}_{\lambda}\left|B^{-1}\sum_{b=1}^{B}P_{n,b}^{1}D_{CAR}^{c}(P_{n,b}^{1},G_{n,b}^{a},G_{n,\lambda,b}^{c},Q_{n,b},\Psi_{n,b},\theta)\right|,$
(13) $\displaystyle\lambda_{n,\Psi}$
$\displaystyle=\min\left\\{\lambda:\left|B^{-1}\sum_{b=1}^{B}P_{n,b}^{1}D_{CAR}^{\Psi}(P_{n,b}^{1},G_{n,b}^{a},G_{n,b}^{c},\Psi_{n,\lambda,b},\theta)\right|\leq\frac{\sigma_{nh}}{(n\log
n)^{1/2}}\right\\},$ (14)

where $G_{n,b}^{a}$, $G_{n,b}^{c}$, $\Psi_{n,b}$ and $Q_{n,b}$ are cross-
validated highly adaptive lasso estimates of the corresponding nuisance
parameters with the L1-norm bound based on the global cross-validation
selector.

To achieve the Gaussian process convergence results in Theorem 4.3 (with a
fixed $h$), we would require stronger conditions. Specifically, the conditions
listed above must hold uniformly for any $\bm{v}\in\mathcal{V}$, that is,
$\sup_{v\in\mathcal{V}}\left|P_{n}D_{CAR}^{a}(P_{n},G_{n}^{a},Q_{n},\Psi_{n},\theta)\right|=o_{p}(n^{-1/2})$,
$\sup_{v\in\mathcal{V}}\left|P_{n}D_{CAR}^{c}(P_{n},G_{n},Q_{n},\Psi_{n},\theta)\right|=o_{p}(n^{-1/2})$
and
$\sup_{v\in\mathcal{V}}\left|P_{n}D_{CAR}^{\Psi}(P_{n},G_{n},\Psi_{n},\theta)\right|=o_{p}(n^{-1/2})$.
Accordingly, the practical criteria are given by

$\displaystyle\lambda_{n,G^{a}}^{u}$
$\displaystyle=\operatorname*{argmin}_{\lambda}\sup_{v\in\tilde{\mathcal{V}}}\left|B^{-1}\sum_{b=1}^{B}P_{n,b}^{1}D_{CAR}^{a}(P_{n,b}^{1},G_{n,\lambda,b}^{a},G_{n,b}^{c},Q_{n,b},\Psi_{n,b},\theta)\right|,$
(15) $\displaystyle\lambda_{n,G^{c}}^{u}$
$\displaystyle=\operatorname*{argmin}_{\lambda}\sup_{v\in\tilde{\mathcal{V}}}\left|B^{-1}\sum_{b=1}^{B}P_{n,b}^{1}D_{CAR}^{c}(P_{n,b}^{1},G_{n,b}^{a},G_{n,\lambda,b}^{c},Q_{n,b},\Psi_{n,b},\theta)\right|,$
(16) $\displaystyle\lambda_{n,\Psi}^{u}$
$\displaystyle=\min\left\\{\lambda:\sup_{v\in\tilde{\mathcal{V}}}\left|B^{-1}\sum_{b=1}^{B}P_{n,b}^{1}D_{CAR}^{\Psi}(P_{n,b}^{1},G_{n,b}^{a},G_{n,b}^{c},\Psi_{n,\lambda,b},\theta)\right|\leq\frac{\sigma_{nh}}{(n\log
n)^{1/2}}\right\\},$ (17)

where
$\tilde{\mathcal{V}}=(\tilde{\bm{v}}_{01},\cdots,\tilde{\bm{v}}_{0\alpha})$
includes $\alpha$ sample points of $\mathcal{V}$. The asymptotic linearity
result of Theorem 4.3 with a fixed bandwidth $h$ relies on similar
undersmoothing criteria as those listed in (a)–(c) but with $h$ being replaced
by one, and thus, the same practical undersmoothing criteria as(12)-(14) can
be considered.

In Theorems 4.2 and 4.3, we show that the proposed estimators are
asymptotically efficient when they solve the efficient influence function,
that is, $P_{n}D_{CAR}^{\diamond}=o_{p}\left((nh^{r})^{-1/2}\right)$ and
$P_{n}D_{CAR}^{\diamond}=o_{p}(n^{-1/2})$, respectively, for
$\diamond\in\\{a,c,\Psi\\}$. The $\operatorname*{argmin}$ criteria proposed in
this section correspond to the most efficient estimators for a given data.
Increasing $\lambda$ results in decreasing the empirical mean of pseudo score
functions
$\frac{K_{h,\bm{v}_{0}}}{P_{n}K_{h,\bm{v}_{0}}}(Q_{n}-\Psi_{n})\left(G_{n,\lambda}^{a}-I(A=d^{\theta})\right)$
and
$\frac{K_{h,\bm{v}_{0}}I(A=d^{\theta})}{P_{n}K_{h,\bm{v}_{0}}}(Q_{n}-\Psi_{n})\left(G_{n,\lambda}^{c}-\Delta^{c}\right)$
and increasing the variance of $(G_{n,\lambda}^{a})^{-1}$ and
$(G_{n,\lambda}^{c})^{-1}$. At a certain point in the grid of $\lambda$,
decreases in the empirical mean of the pseudo score functions are insufficient
for satisfying the $\operatorname*{argmin}$ conditions, which
$\left|P_{n}D_{CAR}^{\diamond}\right|$ starts increasing on account of
$(G_{n,\lambda}^{a})^{-1}$ and $(G_{n,\lambda}^{c})^{-1}$. Unlike
$\left|P_{n}D_{CAR}^{a}\right|$ and $\left|P_{n}D_{CAR}^{c}\right|$,
$\left|P_{n}D_{CAR}^{\Psi}\right|$ decreases as $\lambda$ increases which
motivates the undersmoothing selection criteria (14) and (17). Specifically,
the sectional variation norm $\lambda_{n,\Psi}$ and $\lambda_{n,\Psi}^{u}$ are
defined as the smallest value (larger than that chosen by the cross-validation
selector) for which the left-hand side of the condition is less than
$\sigma_{nh}/(n\log n)^{1/2}$ where $\sigma_{nh}$ be a consistent estimator of
$\sigma_{0}=[E\\{D^{*2}_{\bm{v}_{0}}(P_{0})\\}]^{1/2}$ with
$D^{*2}_{\bm{v}_{0}}(P_{0})$ being the efficient influence function for
$\Psi_{nh}(\bm{v}_{0},\theta)$. The cutpoint will lead to a bias that is of a
smaller magnitude of the standard error, and thus, won’t affect the coverage.

######  Remark 1.

The variance estimator $\sigma_{nh}^{2}$ is calculated using the efficient
influence function. One can, instead, use a conservative variance estimator
obtained by an influence function where the nuisance parameters are assumed to
be known. That is

$\displaystyle\sigma_{nh}^{2}$
$\displaystyle=P_{n}\left\\{\frac{K_{h,\bm{v}_{0}}}{P_{n}K_{h,\bm{v}_{0}}}\Psi_{n}(\theta)-\Psi_{0}(\bm{v}_{0},\theta)\right\\}^{2},$
(18)

where $\Psi_{0}(\bm{v}_{0},\theta)$ can be replaced by their corresponding
estimators. The variance estimator (18) can also be used in Section 3.4 where
we proposed an adaptive bandwidth selector.

##  4 Theoretical results

The asymptotic properties of our estimator relies on the following
assumptions.

######  Assumption 1 (Complexity of nuisance functions).

The functions $Q_{0}$, $\Psi_{0}$, $G_{0}^{a}$ and $G_{0}^{c}$ are càdlàg with
finite sectional variation norm.

Assumption 1 characterizes the global smoothness assumption of the true
functions which is much weaker than than local smoothness assumptions like
those characterized by an assumption utilizing Hölder balls (e.g., Robins et
al., 2008, 2017; Mukherjee, Newey and Robins, 2017). The assumption
facilitates the fast convergence rate of $n^{-1/3}\log(n)^{p/2}$ obtainable by
the highly adaptive lasso (regardless of dimensionality $p$) (van der Laan,
2017; van der Laan and Bibaut, 2017; Bibaut and van der Laan, 2019). The
latter rate is generally faster than the typical minimax rates that appear
under differing, but similarly positioned, assumptions regarding function
classes in nonparametric estimation. In fact, we expect that almost all true
conditional mean models that we have to deal with in real data will satisfy
Assumption 1. The sectional variation is defined in more detail in Appendix B.

######  Assumption 2.

(causal inference and censoring assumptions) For $i=1,\dotsc,n$,

  * (a)

_Consistency_ of the observed outcome: $T_{i}=T_{i}^{A_{i}}$;

  * (b)

_Unconfoundedness_ of the treatment and censoring mechanisms: $A_{i}\perp
T_{i}^{a}|\mathbf{W}_{i},~{}\forall a\in\\{0,1\\}$, and $\Delta_{i}^{c}\perp
T_{i}\mid A_{i},\bm{W}_{i}$;

  * (c)

_Positivity_ of the treatment and censoring mechanisms:
$\min_{a,\delta}G(A=a,\Delta^{c}=\delta\mid\mathbf{W})>\gamma$ for all
$\mathbf{W}\in\mathcal{W}$ where $\gamma$ is a small positive constant.

Assumption 2 lists standard causal identifiability assumptions. The
consistency assumption ensures that a subject’s potential outcome is not
affected by other subjects’ treatment levels and there are no different
versions of treatment. Assumption 2b is a version of the coarsening at random
assumption that imposes conditional independence between (1) treatment and
potential outcomes given the measured covariates; and (2) censoring indicator
and the observed outcome given the treatment and measured covariates.
Positivity states that every subject has some chance of receiving treatment
level $a$ with censoring status $\delta$ regardless of covariates. Using
Assumption 2 the full data risk of $\Psi$ can be identified using the observed
data as

$E_{P_{X}}\int_{\theta}\ell_{\theta}(\Psi,x)dF(\theta)=E_{P_{0}}L_{G_{0}}(\Psi),$

where $L_{G_{0}}(\Psi)$ and $\ell_{\theta}(\Psi,x)$ were defined in Section
2.1. Even if the Assumption 2a-b do not hold, it may still be of interest to
estimate $E_{P_{X}}\int_{\theta}\ell_{\theta}(\Psi,x)dF(\theta)$ using
$P_{n}L_{G_{0}}(\Psi)$ as an adjusted measure of association.

######  Remark 2.

The censoring at random assumption implies that $C\perp T\mid A,\bm{W}$ on
$C<T$ which is weaker than $C\perp T\mid A,\bm{W}$ (Van der Laan, Laan and
Robins, 2003, Chapter 1). However, because the additional restriction can not
be identified using the observed data, it will not impose any restriction on
the tangent space and thus, the statistical model remains nonparametric.

######  Assumption 3 (Undersmoothing).

The level of undersmoothing satisfies the following conditions:

  * (a)

Fixed $h$: Let ${f}_{\phi}^{\Psi}$, ${f}_{\phi}^{a}$ and ${f}_{\phi}^{c}$ be
the projections of $f^{\Psi}=\Delta^{c}I(A=d^{\theta})/G_{0}$,
$f^{a}=(Q_{0}-\Psi_{0})/G_{0}$ and
$f^{c}=I(A=d^{\theta})(Q_{0}-\Psi_{0})/G_{0}$ onto a linear span of basis
functions $\phi_{s,j}$ in $L^{2}(P)$, for $\phi_{s,j}$ satisfying conditions
(23), (24) and (25) of Lemma A.2. Then, $\lVert
f^{\Psi}-{f}_{\phi}^{\Psi}\rVert_{2,\mu}=O_{p}(n^{-1/4})$, $\lVert
f^{a}-{f}_{\phi}^{a}\rVert_{2,\mu}=O_{p}(n^{-1/4})$ and $\lVert
f^{c}-{f}_{\phi}^{c}\rVert_{2,\mu}=O_{p}(n^{-1/4})$.

  * (b)

Converging $h$: The functions $f_{h}^{\Psi}=K_{h,\bm{v}_{0}}f^{\Psi}$,
$f_{h}^{a}=K_{h,\bm{v}_{0}}f^{a}$ and $f_{h}^{c}=K_{h,\bm{v}_{0}}f^{c}$ are
such that $\tilde{f}^{\Psi}=h^{r}f_{h}^{\Psi}$, $\tilde{f}^{a}=h^{r}f_{h}^{a}$
and $\tilde{f}^{c}=h^{r}f_{h}^{c}$ are càdlàg with finite sectional variation
norm. Let $\tilde{f}_{\phi}^{\Psi}$, $\tilde{f}_{\phi}^{a}$ and
$\tilde{f}_{\phi}^{c}$ be projections of $\tilde{f}^{\Psi}$, $\tilde{f}^{a}$
and $\tilde{f}^{c}$ onto the linear span of the basis functions $\phi_{s,j}$
in $L^{2}(P)$, where $\phi_{s,j}$ satisfies condition (20), (21) and (20) of
Lemma A.1. Then,
$\lVert\tilde{f}^{\Psi}-\tilde{f}_{\phi}^{\Psi}\rVert_{2,\mu}=O_{p}(n^{-1/3})$,
$\lVert\tilde{f}^{a}-\tilde{f}_{\phi}^{a}\rVert_{2,\mu}=O_{p}(n^{-1/3})$ and
$\lVert\tilde{f}^{c}-\tilde{f}_{\phi}^{c}\rVert_{2,\mu}=O_{p}(n^{-1/3})$.

where $\mu$ is an appropriate measure (i.e., the Lebesgue or the counting
measure).

Assumption 3 is arguably the most important of the conditions to retain the
asymptotic linearity of $\Psi_{nh}(\bm{v}_{0},\theta)$. Assumption 3 (a)
states that when the bandwidth $h$ is fixed and the estimated coarsening
mechanisms (i.e., $G^{a}$ and $G^{c}$) and $\Psi_{n}$ are sufficiently
undersmoothed, the generated features can approximate functions $f^{a}$,
$f^{c}$ and $f^{\Psi}$ with $n^{-1/3}$ rate. Lemma A.2 in the supplementary
material provides theoretical undersmoothing conditions required. Assumption 3
(b) correspond to the case where the bandwidth is allows to converge to zero
and requires similar conditions as in part (a). Under Assumption 1,
${f}_{\phi}^{\Psi}$, ${f}_{\phi}^{a}$ and ${f}_{\phi}^{c}$ (and
$\tilde{f}^{\Psi}$, $\tilde{f}^{a}$ and $\tilde{f}^{c}$) will fall into the
class of càdlàg functions with finite sectional variation norm. Hence these
functions can be approximated using the highly adaptive lasso generated basis
functions at $n^{-1/3}$ rate (up to a $\log n$ factor). Importantly, the
required rate in Assumption 3 (a) is $O_{p}(n^{-1/4})$ which is slower than
the $O_{p}(n^{-1/3})$ rate obtained by the highly adaptive lasso estimator. In
Remark 4, we discuss that, under a certain smoothness assumption, the required
rate in Assumption 3 (b) can also be slowed down to $O_{p}(n^{-1/4})$.
Consequently, this key assumption is not particularly strong.

######  Assumption 4.

There exist constants $\kappa\geq 0$ and $C>0$ such that, for all $l>0$,
$pr(0<|\bm{S}^{\top}{{\theta}}_{0}|<l)\leq Cl^{\kappa}$.

Assumption 4 is the margin assumption of Audibert and Tsybakov (2007) which
can also be found in Tsybakov (2004). The assumption is needed to derive the
rate of convergence of
$\theta_{nh}(\bm{v})=\operatorname*{argmin}_{\theta\in\Theta}\Psi_{nh}(\bm{v},\theta)$
(Theorem 4.4). The extreme case of $\kappa=\infty$ corresponds to a case where
there is a margin around zero, that is,
$pr(0<|\bm{S}^{\top}{{\theta}}_{0}|<l)=0$. As shown in the proof of Theorem
4.4 the latter leads to the fastest rate of convergence for
$\theta_{nh}(\bm{v})$.

######  Theorem 4.1.

Let $\Psi$ be a functional parameter that is identified as
$\Psi_{0}=\operatorname*{argmin}_{\Psi\in\mathcal{F}_{\lambda}}P_{0}L_{G_{0}}(\Psi)$
where the loss function $L$ is defined in 1. Assume the functional parameter
space $\mathcal{F}_{\lambda}$ contained in all $p$-variate cadlag functions
with a bound on its variation norm $\lambda$. Let
$\Psi_{n}=\operatorname*{argmin}_{\Psi\in\mathcal{F}_{\lambda_{n}}}P_{n}L_{G_{n}}(\Psi)$
where $\lambda_{n}$ is the cross-validated selector of $\lambda$ and let
$d_{0}(\Psi_{n},\Psi_{0})=P_{0}L_{G_{0}}(\Psi_{n})-P_{0}L_{G_{0}}(\Psi_{0})$.
Also assume
$\tilde{\mathcal{F}}_{\tilde{\lambda}}=\\{L_{G_{0}}(\Psi)-L_{G_{0}}(\Psi_{0}):\Psi\in\mathcal{F}_{\lambda}\\}$
falls in to the class of cadlag functions with a bound on its variation norm
$\tilde{\lambda}$. Then, under Assumptions 1 and 2,
$d_{0}(\Psi_{n},\Psi_{0})=O_{p}(n^{-2/3}(\log n)^{4(p-1)/3})$.

Theorem 4.1 generalizes the results of van der Laan (2017) to settings where
the estimation of $\Psi_{0}$ involves estimation of nuisance parameters
$G_{0}$. The proof of the results requires techniques to handle the additional
complexities due to the presence of such nuisance parameters. The theorem
relies on Assumptions 1 without the need for undersmoothing the nuisance
function estimators.

######  Theorem 4.2.

Let $\Psi_{0}(\bm{v}_{0},\theta)=E\\{I(Y^{\theta}>t)\mid\bm{V}=\bm{v}_{0}\\}$
be $J$-times continuously differentiable at $\bm{v}_{0}$. Suppose the support
of $\mathbf{W}$ is uniformly bounded, i.e., $\mathcal{W}\subseteq[0,\tau]^{p}$
for some finite constant $\tau$. Let $\Psi_{n,\lambda_{n}^{\Psi}}$,
$G^{a}_{n,\lambda_{n}^{a}}$ and $G^{c}_{n,\lambda_{n}^{c}}$ be highly adaptive
lasso estimators of $\Psi_{0}$, $G_{0}^{a}$ and $G_{0}^{c}$ with $L_{1}$-norm
bounds equal to $\lambda_{n}^{\Psi}$, $\lambda_{n}^{a}$ and $\lambda_{n}^{c}$,
respectively. Assuming that

  * (i)

the bandwidth $h$ satisfies $h^{r}\rightarrow 0$ and
$nh^{r3/2}\rightarrow\infty$,

  * (ii)

Assumptions 1, 2 and and 3 hold,

we have

$\Psi_{nh}(\bm{v}_{0},\theta)-\Psi_{0}(\bm{v}_{0},\theta)=(P_{n}-P_{0})D_{\bm{v}_{0}}^{*}(P_{0})+h^{J}B_{0}(J,\bm{v}_{0})+o_{p}\left((nh^{r})^{-1/2}\right),$

where $\Psi_{nh}(\bm{v}_{0},\theta)$ is obtained using a ($J-1$)-orthogonal
kernel,
$G_{n,\lambda_{n}^{a},\lambda_{n}^{c}}=G^{a}_{n,\lambda_{n}^{a}}G^{c}_{n,\lambda_{n}^{c}}$
and $D_{\bm{v}_{0}}^{*}(P_{0})$ is the corresponding efficient influence
function defined in the proof of the theorem. Moreover, assuming that
$nh^{r+2J}\rightarrow
0$,$(nh^{r})^{1/2}\\{\Psi_{nh}(\bm{v}_{0},\theta)-\Psi_{0}(\bm{v}_{0},\theta)\\}$
converges to a mean zero normal random variable.

Theorem 4.2 gives conditions under which the proposed kernel estimator
$\Psi_{nh}(\bm{v}_{0},\theta)$ is consistent for
$\Psi_{0}(\bm{v}_{0},\theta)$, and it also gives the corresponding rate of
convergence. In general this result follows if the bandwidth decreases with
sample size sufficiently slowly, and if the nuisance functions $G$ is
estimated sufficiently well. The standard local linear kernel smoothing
requires that $nh^{3r}\rightarrow\infty$ to control the variance (Wasserman,
2006, Chapter 5). Using the entropy number of the class of càdlàg functions
with finite sectional variation norm, we showed that bandwidth can converge to
zero at much faster rate (i.e., $nh^{3r/2}\rightarrow\infty$) than the
standard results. Condition (ii) is the standard causal assumptions. Condition
(iii) indicates the level of undersmoothing of the nuisance parameter
estimators required to achieve asymptotic linearity of our estimator. Section
S6 in the supplementary material provides uniform convergence results with a
rate for our regimen-response curve estimator $\Psi_{nh}(\bm{v}_{0},\theta)$
for all $\bm{v}_{0}\in\mathcal{V}$.

######  Remark 3 (The optimal bandwidth rate).

Theorem 4.2 shows that in order to eliminate the bias term, we must
undersmooth the kernel such that $nh^{r+2J}\rightarrow 0$ (i.e.,
$h<n^{-\frac{1}{r+2J}}$). In combination with rate condition (i) in the
theorem, we have $n^{-\frac{2}{3r}}<h<n^{-\frac{1}{r+2J}}$. Note that for the
latter constraint to hold we must have $J>r/4$. Also, the theoretical
undersmoothing conditions in Lemma A.1 require $J>r$. Hence, the constraint
implies that the optimal bandwidth rate is $h_{n}=n^{-\frac{1}{r+2J}}$ with
$J>r$.

######  Remark 4 (Assumption 3 (b)).

We can allow slower rate of convergence of $O_{p}(n^{-1/4})$ in Assumption 3
(b) (i.e.,
$\lVert\tilde{f}^{\Psi}-\tilde{f}_{\phi}^{\Psi}\rVert_{2,\mu}=O_{p}(n^{-1/4})$,
$\lVert\tilde{f}^{a}-\tilde{f}_{\phi}^{a}\rVert_{2,\mu}=O_{p}(n^{-1/4})$ and
$\lVert\tilde{f}^{c}-\tilde{f}_{\phi}^{c}\rVert_{2,\mu}=O_{p}(n^{-1/4})$) at
the cost of requiring higher level of smoothness for
$\Psi_{0}(\bm{v}_{0},\theta)=E\\{I(Y^{\theta}>t)\mid\bm{V}=\bm{v}_{0}\\}$.
Specifically, the optimal bandwidth rate would be $h_{n}=n^{-\frac{1}{r+2J}}$
with $J>\frac{5}{2}r$.

######  Theorem 4.3.

Let
$\Psi_{0h}(\bm{v}_{0},\theta)=(P_{0}K_{h,\bm{v}_{0}})^{-1}\\{K_{h,\bm{v}_{0}}\Psi_{0}(\theta)\\}$.
Suppose the support of $\bm{W}$ is uniformly bounded, i.e.,
$\mathcal{W}\subseteq[0,\tau]^{p}$ for some finite constant $\tau$. Let
$\Psi_{n,\lambda_{n}^{\Psi}}$, $G^{a}_{n,\lambda_{n}^{a}}$ and
$G^{c}_{n,\lambda_{n}^{c}}$ be highly adaptive lasso estimators of $\Psi_{0}$,
$G_{0}^{a}$ and $G_{0}^{c}$ with $L_{1}$-norm bounds equal to
$\lambda_{n}^{\Psi}$, $\lambda_{n}^{a}$ and $\lambda_{n}^{c}$, respectively.
Under Assumptions 1, 2 and 3, the estimator $\Psi_{nh}(\bm{v}_{0},\theta)$
will be asymptotically linear. That is

$\Psi_{nh}(\bm{v}_{0},\theta)-\Psi_{0h}(\bm{v}_{0},\theta)=(P_{n}-P_{0})D_{h,\bm{v}_{0}}^{*}(P_{0})+o_{p}(n^{-1/2}),$

where
$G_{n,\lambda_{n}^{a},\lambda_{n}^{c}}=G^{a}_{n,\lambda_{n}^{a}}G^{c}_{n,\lambda_{n}^{c}}$
and $D_{h,\bm{v}_{0}}^{*}(P_{0})$ is the corresponding efficient influence
function defined in the proof of the theorem. Moreover,
$\sqrt{n}\\{\Psi_{nh}(\theta)-\Psi_{0h}(\theta)\\}$ converges weakly as a
random element of the cadlag function space endowed with the supremum norm to
a Gaussian process with covariance structure implied by the covariance
function
$\rho(\bm{v},\bm{v}^{\prime})=P_{0}D^{*}_{h,\bm{v}}(P_{0})D^{*}_{h,\bm{v}^{\prime}}(P_{0})$.

Theorem 4.3 strengthens the results of Theorem 4.2 when the bandwidth is
fixed. Specifically, it shows that for a fixed bandwidth $h>0$,
$\sqrt{n}\\{\Psi_{nh}(\theta)-\Psi_{0h}(\theta)\\}$ convergences to a Gaussian
process. This is an important results because in enables us to construct
simultaneous confidence intervals for the entire curve
$\Psi_{0h}(\bm{v},\theta),$ $\bm{v}\in\mathcal{V}$ and $\theta\in\Theta$.

######  Theorem 4.4.

Let
$\theta_{0}(\bm{v})=\operatorname*{argmin}_{\theta\in\Theta}\Psi_{0}(\bm{v},\theta)$
and
$\theta_{nh}(\bm{v})=\operatorname*{argmin}_{\theta\in\Theta}\Psi_{nh}(\bm{v},\theta)$
where both $\Psi_{0}$ and $\Psi_{nh}$ are defined in Theorem 4.2. Then, under
Assumptions 1-4,
$\|\theta_{nh}(\bm{v}_{0})-\theta_{0}(\bm{v}_{0})\|_{2}=O_{p}\left(n^{\frac{(r-4J)(2\kappa+4)}{(8J+4r)(3\kappa+8)}}\right)$.

Characterizing the minimizer (or in general optimizer) of the estimated
regimen-response curve is of interest as it forms an optimal individualized
decision rule. Under the margin assumption (i.e., Assumption 4) and using a
novel iterative empirical process theory, in Theorem 4.4, we derive the
convergence rate of a minimizer of a nonparametrically estimated function
$\Psi_{nh}(\bm{v},\theta)$.

##  5 Simulation Studies

We examine the performance of our estimator through two simulation studies.
Within these studies, we compare the results of our estimator against
theoretical values and random forest based estimation, demonstrating the
practical benefits of using undersmoothed highly adaptive lasso based
estimation when the data generating mechanism is unknown.

In both scenarios $W\sim\text{Unif}(-1,1)$ and $V\sim\text{Unif}(-1.3,1.3)$.
Additionally $\Delta|W,V\sim\text{Bernoulli}\\{\text{expit}(0.5W-0.5V)\\}$ and
$Y|A,W,V\sim\text{Bernoulli}\\{\text{expit}(0.5W+V^{2}+2VW-2AVW)\\}.$ For each
scenario we sample $n\in\\{240,480,960,2400\\}$ observations, applying our
estimator for each sample. The results displayed below give average effects
360 replicates of each scenario and sample size pair where we set
$\xi(\theta,\bm{V})=1$. The latter implies that $\Psi_{n}$ is obtained by
solving a set of simpler score equations (compared with
$\xi(\theta,\bm{V})=G\\{(A,\Delta^{c})\mid\bm{V}\\}$) which may call for more
undersmoothing in order to solve $D_{CAR}^{\Psi}(P_{n},G_{n},\Psi_{n})$.

The fundamental difference between the scenarios comes in the treatment
assignment mechanism. In Scenario 1, we have $A\sim\text{Bernoulli}(0.5)$,
giving us a randomized trial scenario, whereas in Scenario 2, we have
$A\sim\text{Bernoulli}\\{\text{expit}(0.7W+0.5V+0.5VW)\\}.$ The highly
adaptive lasso estimate of $G^{a}(A\mid W,V)$ requires solving fewer score
functions in Scenario 1 compared with Scenario 2. In fact the only score
equation that to be solved in the former is $P_{n}\\{A-G^{a}_{n}(A\mid
W,V)\\}=0$. Hence, higher level of undersmooting is required to solve
$D_{CAR}^{a}(P_{n},G_{n}^{a},Q_{n},\Psi_{n},\theta)$. In general, the simpler
the $G\\{(A,\Delta^{c})\mid W,V\\}$, the more undersmooting is needs to solve
$D_{CAR}^{a}(P_{n},G_{n}^{a},Q_{n},\Psi_{n},\theta)$ and
$D_{CAR}^{c}(P_{n},G_{n}^{a},Q_{n},\Psi_{n},\theta)$.

In all simulations we are targeting the estimands $\Psi_{0h}({v}_{0},\theta)$
and $\Psi_{0}({v}_{0},\theta)$ using the estimator $\Psi_{nh}({v}_{0},\theta)$
where ${v}_{0}=0.5$ and we consider the bandwidths $h\in
n^{-1/3}\times\\{0.5,1,1.5,2,2.5,3,4\\}$. Recall that the optimal bandwidth
rate is $n^{\frac{-1}{2J+r}}$ with $J>r$. In our simulation studies $r=1$ and
we are setting $J=r$ leading to a rate $n^{-1/3}$ which is faster than
$n^{-1/5}$ obtained by setting $J=2>r$. The choice of rate matters in cases
where the target parameter is $\Psi_{0}$ and the goal is to show that the
scaled biases remain constant and converges are nominal even under the faster
rate of $n^{-1/3}$ (Figures 2 and 3). Coverage results are also given for
$\Psi_{0}({v}_{0},\theta)$ where for each sample $h_{n}$ is chosen from the
possible bandwidths as described in Section 3.4. In all cases highly adaptive
lasso was implemented using the R package hal9001 considering quadratic basis
functions for up to 2-way interactions between all predictor variables.

Figure 1: Simulation studies Scenario 1: The target parameter of interest is
$\Psi_{0h}$. The plots show the scaled bias and coverage rate when the weight
functions are estimated using an undersmoothed highly adaptive lasso (HAL) and
a random forest (RF).

We estimate the weight functions (i.e., propensity score and censoring
indicator models) using an undersmoothed highly adaptive lasso (HAL) and a
random forest (RF). In both cases, $\Psi_{n}$ is a highly adaptive lasso
estimate of $\Psi_{0}$ obtained based on the loss function (1). Figure 1 shows
the scaled bias and the coverage rate of the estimators when the target
parameter is $\Psi_{0h}$. The results confirm our theoretical result presented
in Theorem 4.3. The scaled bias if the estimator with undersmoothed HAL is
considerably smaller than the one obtained using RF. Importantly, while the
coverage rate of undersmoothed HAL estimator remains close to the nominal rate
of 95%, the coverage rate of RF based estimator sharply declines as sample
size increases. To assess our theoretical result of Theorem 4.2, we plotted
the the scaled bias and the coverage rate of the estimators when the target
parameter is $\Psi_{0}$. Figure 2 shows that the undersmoothed HAL outperforms
RF based estimators. Figures S5-S6 in the Supplementary Material show that
similar results hold in Scenario 2. Finally, Figure 3 shows that our proposed
data-adaptive bandwidth selector performs very well when combined with
undersmoothed HAL estimators of weight functions. Specifically, even with
smaller sample size of $n=240$, our approach results in nearly 95% coverage
rate and as the sample size increases it quickly recovers the nominal rate.
Figure 4 displays $\Psi_{nh}(v_{0}=0.5,\theta)$ against $\theta$ values. As
the sample size increases the bias reduces across all the bandwidth values.
Figure S7 in the Supplementary Material shows the scaled bias of $\theta_{nh}$
when the weight functions are estimated using an undersmoothed highly adaptive
lasso. The plot shows that the scaled bias is relatively constant across
different sample sizes thereby confirming our result in Theorem 4.4 (by
setting $J=r=1$ and $\kappa\rightarrow\infty$).

Figure 2: Simulation studies Scenario 1: The target parameter of interest is
$\Psi_{0}$. The plots show the scaled bias and coverage rate when the weight
functions are estimated using an undersmoothed highly adaptive lasso (HAL) and
a random forest (RF). Figure 3: Simulation studies: The target parameter of
interest is $\Psi_{0}$. The dotted line shows the nominal rate. The optimal
bandwidth $h$ is selected using the proposed data adaptive approach. The plots
show the coverage rates in Scenarios 1 and 2 where the weight functions are
estimated using an undersmoothed highly adaptive lasso (HAL) and a random
forest (RF). Figure 4: Simulation studies Scenario 1: The target parameter of
interest is $\Psi_{0}$ (black solid line). The plots show $\Psi_{nh}$ for
different bandwidths when the weight functions are estimated using an
undersmoothed highly adaptive lasso and $v_{0}=0.5$.

##  6 Concluding Remarks

###  6.1 Multi-stage decision making

A natural and important extension of our work is to learn dynamic treatment
regimens in multi-stage problems. In settings with many decision points (e.g.,
mobile health), it would be beneficial to consider a stochastic class of
decision rules rather than a deterministic one (Kennedy, 2019; Luckett et al.,
2019).

###  6.2 Score based undersmoothing criteria

The practical undersmoothing criteria proposed in Section 3.5 requires
calculating the efficient influence function and deriving the relevant
$D_{CAR}$ terms. The latter can be challenging and tedious, for example, in
multi-stage settings (i.e., time-varying treatment) with many decision points.
This motivates the derivation of score based criteria which requires only the
score of the corresponding nuisance function. This will, in fact, resemble
closely the theoretical undersmoothing conditions listed in Lemmas A.1, A.1
and A.3 in Appendix and Supplementary materials. For example, for the
treatment assignment model, we can consider

$\lambda_{n,G^{a}}=\operatorname*{argmin}_{\lambda}B^{-1}\sum_{b=1}^{B}\left[\sum_{(s,j)\in\mathcal{J}_{n}}\frac{1}{\lVert\beta_{n,\lambda,b}\rVert_{L_{1}}}\bigg{\lvert}P_{n,b}^{1}\tilde{S}_{s,j}(\phi,G_{n,\lambda,b})\bigg{\rvert}\right],$
(19)

in which
$\lVert\beta_{n,\lambda}\rVert_{L_{1}}=\lvert\beta_{n,\lambda,0}\rvert+\sum_{s\subset\\{1,\ldots,d\\}}\sum_{j=1}^{n}\lvert\beta_{n,\lambda,s,j}\rvert$
is the $L_{1}$-norm of the coefficients $\beta_{n,\lambda,s,j}$ in the highly
adaptive lasso estimator $G_{n,\lambda}$ for a given $\lambda$, and
$\tilde{S}_{s,j}(\phi,G_{n,\lambda,b})=\phi_{s,j}(W)\\{A-G_{n,\lambda,b}(1\mid
W)\\}\\{G_{n,\lambda,b}(1\mid W)\\}^{-1}$. Note that
$S_{s,j}(\phi,G_{n,\lambda,b})$ involves the score equation for the propensity
score model (i.e., $\phi_{s,j}(W)\\{A-G_{n,\lambda,b}(1\mid W)\\}$) and the
inverse of the estimate propensity score. A more detailed discussion of the
score based criteria is provided in Section 4.2 of Ertefaie, Hejazi and van
der Laan (2022). However, the theoretical properties of these criteria and
their performance in regimen-curve estimation problems are unknown and merit
further research.

###  6.3 Near positively violation

In general, near positivity violation can negatively impact the performance of
causal estimators. In particular, it can lead to an increased bias and
variance in inverse probability weighted estimators and undersmoothing may
exacerbate the issue. One possible remedy is to consider truncation where
individuals whose estimated weight functions are smaller than a given cutpoint
are removed from the analyses (Crump et al., 2009). There are also other
methods to trim the population in an interpretable way (Traskin and Small,
2011; Fogarty et al., 2016). Another possible way to improve the performance
of our approach under near positivity violation is to fit a highly adaptive
lasso by enforcing the positivity $\min(G^{a}_{n},1-G^{a}_{n})>\gamma$ and
$\min(G^{c}_{n},1-G^{c}_{n})>\gamma$, for a positive constant $\gamma$. One
can even make $\gamma$ data dependent (i.e., consider $\gamma_{n}$). This is
an interesting direction from both methodological and practical point of view.

###  6.4 Variable selection

In our proposed method, we assume that the set of variables included in
regimen response-curve function and the decision rule (i.e., $\bm{V}$ and
$\bm{S}$, respectively) are a priori specified. Particularly, in health
research, variable selection in decision making is important as it may reduce
the treatment costs and burden on patients. Moreover, the quality of the
constructed rules can be severely hampered by inclusion of many extraneous
variables (Jones, Ertefaie and Strawderman, 2022). Therefore, it will be
interesting to systematically investigate how to perform variable selection in
regimen response-curve function and the decision rule and to provide
theoretical guarantees.

##  Appendix A Theoretical undersmoothing conditions

######  Lemma A.1 (Theoretical undersmoothing conditions when
$h_{n}\rightarrow 0$).

Let $\Psi\equiv\Psi(\theta,\bm{V})$, $G^{a}\equiv G(P)(A\mid\bm{W})$ and
$G^{c}\equiv G(P)(\Delta^{c}\mid\bm{W},A)$ denote $E(Y^{d^{\theta}}\mid V)$,
the treatment and censoring mechanism under an arbitrary distribution
$P\in\mathcal{M}$. Let $\Psi_{n,\lambda_{n}^{\Psi}}$,
$G^{a}_{n,\lambda_{n}^{a}}$ and $G^{c}_{n,\lambda_{n}^{c}}$ be the highly
adaptive lasso estimators of $\Psi$, $G^{a}$ and $G^{c}$ using $L_{1}$-norm
bound $\lambda_{n}^{\Psi}$, $\lambda_{n}^{a}$ and $\lambda_{n}^{c}$,
respectively. Choosing $\lambda_{n}^{\Psi}$, $\lambda_{n}^{a}$ and
$\lambda_{n}^{c}$ such that

$\displaystyle\min_{(s,j)\in\mathcal{J}_{n}^{\Psi}}{\bigg{\|}}P_{n}\frac{d}{d\Psi_{n,\lambda_{n}^{\Psi}}}L(\Psi_{n,\lambda_{n}^{\Psi}})(\phi_{s,j}){\bigg{\|}}$
$\displaystyle=o_{p}\left(n^{-1/2}h^{r/2}\right),$ (20)
$\displaystyle\min_{(s,j)\in\mathcal{J}_{n}^{a}}{\bigg{\|}}P_{n}\frac{d}{d\text{logit
}G_{n,\lambda_{n}^{a}}^{a}}L(\text{logit
}G_{n,\lambda_{n}^{a}}^{a})(\phi_{s,j}){\bigg{\|}}$
$\displaystyle=o_{p}\left(n^{-1/2}h^{r/2}\right),$ (21)
$\displaystyle\min_{(s,j)\in\mathcal{J}_{n}^{c}}{\bigg{\|}}P_{n}\frac{d}{d\text{logit
}G_{n,\lambda_{n}^{c}}^{c}}L(\text{logit
}G_{n,\lambda_{n}^{c}}^{c})(\phi_{s,j}){\bigg{\|}}$
$\displaystyle=o_{p}\left(n^{-1/2}h^{r/2}\right),$ (22)

where, in (20), $L(\cdot)$ is the loss function (1), and, in (21), $L(\cdot)$
is log-likelihood loss. Also, $\mathcal{J}_{n}^{.}$ is a set of indices for
the basis functions with nonzero coefficients in the corresponding model. Let
$D^{\Psi}(f_{h}^{\Psi},\Psi_{n})=f_{h}^{\Psi}\cdot\\{I(Y>t)-\Psi_{n}\\}$,
$D^{a}(f_{h}^{a},G_{n}^{a})=f_{h}^{a}\cdot(A-G_{n}^{a})$, and
$D^{c}(f^{c},G_{n}^{c})=f_{h}^{c}\cdot(\Delta^{c}-G_{n}^{c})$. The functions
$f_{h}^{\Psi}$, $f_{h}^{a}$ and $f_{h}^{c}$ are such that
$\tilde{f}^{\Psi}=h^{r}f_{h}^{\Psi}$, $\tilde{f}^{a}=h^{r}f_{h}^{a}$ and
$\tilde{f}^{c}=h^{r}f_{h}^{c}$ are càdlàg with finite sectional variation
norm. Let $\tilde{f}_{\phi}^{\Psi}$, $\tilde{f}_{\phi}^{a}$ and
$\tilde{f}_{\phi}^{c}$ be projections of $\tilde{f}^{\Psi}$, $\tilde{f}^{a}$
and $\tilde{f}^{c}$ onto the linear span of the basis functions $\phi_{s,j}$
in $L^{2}(P)$, where $\phi_{s,j}$ satisfies condition (20) and (21),
respectively. Then, under the optimal bandwidth rate $h=n^{\frac{-1}{r+2J}}$
and $J>r$, we have
$P_{n}D^{\Psi}({f}_{h}^{\Psi},\Psi_{n})=o_{p}\left((nh^{r})^{-1/2}\right)$,
$P_{n}D^{a}({f}_{h}^{a},G_{n}^{a})=o_{p}\left((nh^{r})^{-1/2}\right)$ and
$P_{n}D^{c}({f}_{h}^{c},G_{n}^{c})=o_{p}\left((nh^{r})^{-1/2}\right)$.

######  Lemma A.2 (Theoretical undersmoothing conditions for a fixed bandwidth
$h$).

Let $\Psi\equiv\Psi(\theta,\bm{V})$, $G^{a}\equiv G(P)(A\mid\bm{W})$ and
$G^{c}\equiv G(P)(\Delta^{c}\mid\bm{W},A)$ denote $E(Y^{d^{\theta}}\mid V)$,
the treatment and censoring mechanism under an arbitrary distribution
$P\in\mathcal{M}$. Let $\Psi_{n,\lambda_{n}^{\Psi}}$,
$G^{a}_{n,\lambda_{n}^{a}}$ and $G^{c}_{n,\lambda_{n}^{c}}$ be the highly
adaptive lasso estimators of $\Psi$, $G^{a}$ and $G^{c}$ using $L_{1}$-norm
bound $\lambda_{n}^{\Psi}$, $\lambda_{n}^{a}$ and $\lambda_{n}^{c}$,
respectively. Choosing $\lambda_{n}^{\Psi}$, $\lambda_{n}^{a}$ and
$\lambda_{n}^{c}$ such that

$\displaystyle\min_{(s,j)\in\mathcal{J}_{n}^{\Psi}}{\bigg{\|}}P_{n}\frac{d}{d\Psi_{n,\lambda_{n}^{\Psi}}}L(\Psi_{n,\lambda_{n}^{\Psi}})(\phi_{s,j}){\bigg{\|}}$
$\displaystyle=o_{p}\left(n^{-1/2}\right),$ (23)
$\displaystyle\min_{(s,j)\in\mathcal{J}_{n}^{a}}{\bigg{\|}}P_{n}\frac{d}{d\text{logit
}G_{n,\lambda_{n}^{a}}^{a}}L(\text{logit
}G_{n,\lambda_{n}^{a}}^{a})(\phi_{s,j}){\bigg{\|}}$
$\displaystyle=o_{p}\left(n^{-1/2}\right),$ (24)
$\displaystyle\min_{(s,j)\in\mathcal{J}_{n}^{c}}{\bigg{\|}}P_{n}\frac{d}{d\text{logit
}G_{n,\lambda_{n}^{c}}^{c}}L(\text{logit
}G_{n,\lambda_{n}^{c}}^{c})(\phi_{s,j}){\bigg{\|}}$
$\displaystyle=o_{p}\left(n^{-1/2}\right),$ (25)

where, in (23), $L(\cdot)$ is the loss function (1), and, in (24) and (25),
$L(\cdot)$ is log-likelihood loss. Also, $\mathcal{J}_{n}^{.}$ is a set of
indices for the basis functions with nonzero coefficients in the corresponding
model. Let $D^{\Psi}(f^{\Psi},\Psi_{n})=f^{\Psi}\cdot\\{I(Y>t)-\Psi_{n}\\}$,
$D^{a}(f^{a},G_{n}^{a})=f^{a}\cdot(A-G_{n}^{a})$, and
$D^{c}(f^{c},G_{n}^{c})=f^{c}\cdot(\Delta^{c}-G_{n}^{c})$. The functions
$f^{\Psi}$, $f^{a}$ and $f^{c}$ are càdlàg with finite sectional variation
norm. Let $\tilde{f}_{\phi}^{\Psi}$, $\tilde{f}_{\phi}^{a}$ and
$\tilde{f}_{\phi}^{c}$ be projections of $\tilde{f}^{\Psi}$, $\tilde{f}^{a}$
and $\tilde{f}^{c}$ onto the linear span of the basis functions $\phi_{s,j}$
in $L^{2}(P)$, where $\phi_{s,j}$ satisfies condition (23) and (24),
respectively. Then, we have
$P_{n}D^{\Psi}({f}^{\Psi},\Psi_{n})=o_{p}\left(n^{-1/2}\right)$,
$P_{n}D^{a}({f}^{a},G_{n}^{a})=o_{p}\left(n^{-1/2}\right)$ and
$P_{n}D^{c}({f}^{c},G_{n}^{c})=o_{p}\left(n^{-1/2}\right)$.

######  Lemma A.3 (Theoretical undersmoothing conditions for uniform
convergence with fixed bandwidth $h$).

Let $\Psi\equiv\Psi(\theta,\bm{V})$, $G^{a}\equiv G(P)(A\mid\bm{W})$ and
$G^{c}\equiv G(P)(\Delta^{c}\mid\bm{W},A)$ denote $E(Y^{d^{\theta}}\mid V)$,
the treatment and censoring mechanism under an arbitrary distribution
$P\in\mathcal{M}$. Let $\Psi_{n,\lambda_{n}^{\Psi}}$,
$G^{a}_{n,\lambda_{n}^{a}}$ and $G^{c}_{n,\lambda_{n}^{c}}$ be the highly
adaptive lasso estimators of $\Psi$, $G^{a}$ and $G^{c}$ using $L_{1}$-norm
bound $\lambda_{n}^{m}$, $\lambda_{n}^{a}$ and $\lambda_{n}^{c}$,
respectively. Choosing $\lambda_{n}^{\Psi}$, $\lambda_{n}^{a}$ and
$\lambda_{n}^{c}$ such that

$\displaystyle\min_{(s,j)\in\mathcal{J}_{n}^{\Psi}}\sup_{v\in{\mathcal{V}}}{\bigg{\|}}P_{n}\frac{d}{d\Psi_{n,\lambda_{n}^{\Psi}}}L(\Psi_{n,\lambda_{n}^{\Psi}})(\phi_{s,j}){\bigg{\|}}$
$\displaystyle=o_{p}(n^{-1/2}),$ (26)
$\displaystyle\min_{(s,j)\in\mathcal{J}_{n}^{a}}\sup_{v\in{\mathcal{V}}}{\bigg{\|}}P_{n}\frac{d}{d\text{logit
}G_{n,\lambda_{n}^{a}}^{a}}L(\text{logit
}G_{n,\lambda_{n}^{a}}^{a})(\phi_{s,j}){\bigg{\|}}$
$\displaystyle=o_{p}(n^{-1/2}),$ (27)
$\displaystyle\min_{(s,j)\in\mathcal{J}_{n}^{c}}\sup_{v\in{\mathcal{V}}}{\bigg{\|}}P_{n}\frac{d}{d\text{logit
}G_{n,\lambda_{n}^{c}}^{c}}L(\text{logit
}G_{n,\lambda_{n}^{c}}^{c})(\phi_{s,j}){\bigg{\|}}$
$\displaystyle=o_{p}(n^{-1/2}),$ (28)

where, in (26), $L(\cdot)$ is the loss function (1), and, in (27) and (28),
$L(\cdot)$ is log-likelihood loss. Also, $\mathcal{J}_{n}^{.}$ is a set of
indices for the basis functions with nonzero coefficients in the corresponding
model. Let $D^{\Psi}(f^{\Psi},\Psi_{n})=f^{\Psi}\cdot\\{I(Y>t)-\Psi_{n}\\}$,
$D^{a}(f^{a},G_{n}^{a})=f^{a}\cdot(A-G_{n}^{a})$, and
$D^{c}(f^{c},G_{n}^{c})=f^{c}\cdot(\Delta^{c}-G_{n}^{c})$. Here, $f^{\Psi}$,
$f^{a}$ and $f^{c}$ are càdlàg with finite sectional variation norm, and we
let $\tilde{f}^{\Psi}$, $\tilde{f}^{a}$ and $\tilde{f}^{c}$ be projections of
$f^{\Psi}$, $f^{a}$ and $f^{c}$ onto the linear span of the basis functions
$\phi_{s,j}$ in $L^{2}(P)$, where $\phi_{s,j}$ satisfies condition (26) and
(27), respectively. Assuming $\lVert
f^{\Psi}-\tilde{f}^{\Psi}\rVert_{2,P_{0}}=O_{p}(n^{-1/4})$, $\lVert
f^{a}-\tilde{f}^{a}\rVert_{2,P_{0}}=O_{p}(n^{-1/4})$ and $\lVert
f^{c}-\tilde{f}^{c}\rVert_{2,P_{0}}=O_{p}(n^{-1/4})$, it follows that
$P_{n}D^{\Psi}(\tilde{f}^{\Psi},\Psi_{n})=o_{p}(n^{-1/2})$,
$P_{n}D^{a}(\tilde{f}^{a},G_{n}^{a})=o_{p}(n^{-1/2})$,
$P_{n}D^{c}(\tilde{f}^{c},G_{n}^{c})=o_{p}(n^{-1/2})$
$P_{n}D^{\Psi}({f}^{\Psi},\Psi_{n})=o_{p}(n^{-1/2})$,
$P_{n}D^{a}({f}^{a},G_{n}^{a})=o_{p}(n^{-1/2})$ and
$P_{n}D^{c}({f}^{c},G_{n}^{c})=o_{p}(n^{-1/2})$.

##  Appendix B Sectional variation norm

Suppose that the function of interest falls into a càdlàg class of functions
on a p-dimensional cube $[0,\tau]\subset\mathcal{R}^{p}$ with finite sectional
variation norm where $\tau$ is the upper bound of all supports which is
assumed to be finite. The p-dimensional cube $[0,\tau]$ can be represented as
a union of lower dimensional cubes (i.e., $l$-dimensional with $l\leq p$) plus
the origin. That is $[0,\tau]=\\{\cup_{s}(0,\tau_{s}]\\}\cup\\{0\\}$ where
$\cup_{s}$ is over all subsets $s$ of $\\{1,2,\cdots,p\\}$.

Denote $\mathcal{D}[0,\tau]$ as the Banach space of $p$-variate real-valued
càdlàg functions on a cube $[0,\tau]\in\mathbb{R}^{p}$. For each function
$f\in\mathcal{D}[0,\tau]$, we define the $s-$th section $f$ as
$f_{s}(u)=f\\{u_{1}I(1\in s),\cdots,u_{d}\in I(p\in s)\\}$. By definition,
$f_{s}(u)$ varies along the variables in $u_{s}$ according to $f$ while
setting other variables to zero. We define the sectional variation norm of a
given $f$ as

$\lVert f\rVert^{*}_{\zeta}=\lvert
f(0)\rvert+\sum_{s\subset\\{1,\ldots,p\\}}\int_{0_{s}}^{\tau_{s}}\lvert
df_{s}(u)\rvert,$

where the sum is over all subsets of the coordinates $\\{0,1,\ldots,p\\}$. The
term $\int_{0_{s}}^{\tau_{s}}\lvert df_{s}(u)\rvert$ denotes the $s-$specific
variation norm. To clarify the definition, let’s consider a simple case of
$p=1$. In this case, the variation norm of function $f$ over $[0,\tau]$ is
defined as

$\lVert f\rVert^{*}_{\zeta}=\sup_{Q\in\mathcal{Q}}\sum_{i=1}^{n_{q}-1}\lvert
f(c_{i+1})-f(c_{i})\rvert$

where $c_{0}=0$, $Q$ is a partition of $[0,\tau]$ such as
$(0,c_{1}],(c_{1},c_{2}],\cdots,(c_{n-1},c_{n_{q}}=\tau]$ and $\mathcal{Q}$ is
the class of all possible partitions. For differentiable functions $f$, the
sectional variation norm is $\lVert f\rVert^{*}_{\zeta}=\int_{0}^{\tau}\lvert
f^{\prime}(u)\rvert du$.

For $p=2$, the variation norm of function $f$ over
$[0,\tau]\subset\mathcal{R}^{2}$ is defined as

$\lVert
f\rVert^{*}_{\zeta}=\sup_{Q\in\mathcal{Q}}\sum_{i=0}^{n_{q_{1}}-1}\sum_{j=0}^{n_{q_{2}}-1}\lvert
f(u_{i+1},v_{j+1})-f(u_{i},v_{j+1})-f(u_{i+1},v_{j})+f(u_{i},v_{j})\rvert$

where $u_{0}=v_{0}=0$, $Q$ is a partition of $[0,\tau]$ such as
$(0,u_{1}],(u_{1},u_{2}],\cdots,(u_{n_{q_{2}}-1},u_{n_{q_{1}}}=\tau]$,
$(0,v_{1}],(v_{1},v_{2}],\cdots,(v_{n_{q_{2}}-1},v_{n_{q_{2}}}=\tau]$ and
$\mathcal{Q}$ is the class of all possible partitions.

In practice, we choose the support points $u_{s,i}$ to be the observed values
of $\bm{W}$ denoted as $\tilde{w}_{s,1},\cdots,\tilde{w}_{s,n}$. Hence the
indicator basis functions will be
$\phi_{s,i}(w_{s})=I(w_{s}\geq\tilde{w}_{s,i})$. For simplicity lets consider
$p=1$, then, in a hypothetical example, the design matrix based on indicator
basis functions would be

$W$ | $I(W\geq 1.45)$ | $I(W\geq 0.84)$ | $I(W\geq 1.0)$ | $I(W\geq 0.16)$ | $\cdots$
---|---|---|---|---|---
1.45 | 1 | 1 | 0 | 1 | $\cdots$
0.84 | 0 | 1 | 0 | 1 | $\cdots$
1.0 | 1 | 1 | 1 | 1 | $\cdots$
0.16 | 0 | 0 | 0 | 1 | $\cdots$
$\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ |

The first author was supported by National Institute on Drug Abuse under award
number R01DA048764, and National Institute of Neurological Disorders and
Stroke under award numbers R33NS120240 and R61NS12024. The third author was
supported in part by NIH Grant R01AI074345.

Proofs of the theoretical results S1-S7 Supplements S1-S7 contain proofs of
our theoretical results. Additional figures Supplement S8 includes additional
figures for the simulation studies.

## References

  * Athey and Wager (2021) [author] Athey, SusanS. and Wager, StefanS. (2021). Policy learning with observational data. Econometrica 89 133–161.
  * Audibert and Tsybakov (2007) [author] Audibert, Jean-YvesJ.-Y. and Tsybakov, Alexandre BA. B. (2007). Fast learning rates for plug-in classifiers. The Annals of statistics 35 608–633.
  * Benkeser and Van Der Laan (2016) [author] Benkeser, DavidD. and Van Der Laan, MarkM. (2016). The highly adaptive lasso estimator. In Proceedings of the… International Conference on Data Science and Advanced Analytics. IEEE International Conference on Data Science and Advanced Analytics 2016 689\. NIH Public Access.
  * Benkeser et al. (2017) [author] Benkeser, DavidD., Carone, MarcoM., Laan, MJ Van DerM. V. D. and Gilbert, PBP. (2017). Doubly robust nonparametric inference on the average treatment effect. Biometrika 104 863–880.
  * Bibaut and van der Laan (2019) [author] Bibaut, Aurélien FA. F. and van der Laan, Mark JM. J. (2019). Fast rates for empirical risk minimization over c$\backslash$adl$\backslash$ag functions with bounded sectional variation norm. arXiv preprint arXiv:1907.09244.
  * Bickel et al. (1998) [author] Bickel, Peter JP. J., Klaassen, Chris AC. A., Bickel, Peter JP. J., Ritov, YY., Klaassen, JJ., Wellner, Jon AJ. A. and Ritov, YA’AcovY. (1998). Efficient and adaptive estimation for semiparametric models 2\. Springer New York.
  * Chernozhukov et al. (2017) [author] Chernozhukov, VictorV., Chetverikov, DenisD., Demirer, MertM., Duflo, EstherE., Hansen, ChristianC. and Newey, WhitneyW. (2017). Double/debiased/neyman machine learning of treatment effects. American Economic Review 107 261–65.
  * 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.
  * Dudík, Langford and Li (2011) [author] Dudík, MiroslavM., Langford, JohnJ. and Li, LihongL. (2011). Doubly robust policy evaluation and learning. arXiv preprint arXiv:1103.4601.
  * Ertefaie, Hejazi and van der Laan (2022) [author] Ertefaie, AshkanA., Hejazi, Nima SN. S. and van der Laan, Mark JM. J. (2022). Nonparametric inverse probability weighted estimators based on the highly adaptive lasso. Biometrics.
  * Ertefaie et al. (2015) [author] Ertefaie, AshkanA., Wu, TianshuangT., Lynch, Kevin GK. G. and Nahum-Shani, InbalI. (2015). Identifying a set that contains the best dynamic treatment regimes. Biostatistics 17 135–148.
  * Fogarty et al. (2016) [author] Fogarty, Colin BC. B., Mikkelsen, Mark EM. E., Gaieski, David FD. F. and Small, Dylan SD. S. (2016). Discrete optimization for interpretable study populations and randomization inference in an observational study of severe sepsis mortality. Journal of the American Statistical Association 111 447–458.
  * Gill, van der Laan and Wellner (1993) [author] Gill, Richard DR. D., van der Laan, Mark JM. J. and Wellner, Jon AJ. A. (1993). Inefficient estimators of the bivariate survival function for three models. Rijksuniversiteit Utrecht. Mathematisch Instituut.
  * Joffe et al. (2004) [author] Joffe, Marshall MM. M., Ten Have, Thomas RT. R., Feldman, Harold IH. I. and Kimmel, Stephen ES. E. (2004). Model selection, confounder control, and marginal structural models: review and new applications. The American Statistician 58 272–279.
  * Jones, Ertefaie and Strawderman (2022) [author] Jones, JeremiahJ., Ertefaie, AshkanA. and Strawderman, Robert LR. L. (2022). Valid post-selection inference in Robust Q-learning. arXiv preprint arXiv:2208.03233.
  * Kallus (2018) [author] Kallus, NathanN. (2018). Balanced policy evaluation and learning. Advances in neural information processing systems 31\.
  * 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.
  * Klaassen (1987) [author] Klaassen, Chris AJC. A. (1987). Consistent estimation of the influence function of locally asymptotically linear estimators. The Annals of Statistics 15 1548–1562.
  * Liu et al. (2018) [author] Liu, YingY., Wang, YuanjiaY., Kosorok, Michael RM. R., Zhao, YingqiY. and Zeng, DonglinD. (2018). Augmented outcome-weighted learning for estimating optimal dynamic treatment regimens. Statistics in medicine 37 3776–3788.
  * Luckett et al. (2019) [author] Luckett, Daniel JD. J., Laber, Eric BE. B., Kahkoska, Anna RA. R., Maahs, David MD. M., Mayer-Davis, ElizabethE. and Kosorok, Michael RM. R. (2019). Estimating dynamic treatment regimes in mobile health using v-learning. Journal of the American Statistical Association.
  * Mukherjee, Newey and Robins (2017) [author] Mukherjee, RajarshiR., Newey, Whitney KW. K. and Robins, James MJ. M. (2017). Semiparametric efficient empirical higher order influence function estimators. arXiv preprint arXiv:1705.07577.
  * Murphy et al. (2001) [author] Murphy, Susan AS. A., van der Laan, Mark JM. J., Robins, James MJ. M. and Group, Conduct Problems Prevention ResearchC. P. P. R. (2001). Marginal mean models for dynamic regimes. Journal of the American Statistical Association 96 1410–1423.
  * 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.
  * Nie and Wager (2021) [author] Nie, XinkunX. and Wager, StefanS. (2021). Quasi-oracle estimation of heterogeneous treatment effects. Biometrika 108 299–319.
  * Orellana, Rotnitzky and Robins (2010) [author] Orellana, LilianaL., Rotnitzky, AndreaA. and Robins, James MJ. M. (2010). Dynamic regime marginal structural mean models for estimation of optimal dynamic treatment regimes, part I: main content. The international journal of biostatistics 6\.
  * Petersen et al. (2014) [author] Petersen, MayaM., Schwab, JoshuaJ., Gruber, SusanS., Blaser, NelloN., Schomaker, MichaelM. and van der Laan, MarkM. (2014). Targeted maximum likelihood estimation for dynamic and static longitudinal marginal structural working models. Journal of causal inference 2 147–185.
  * Robins (1998) [author] Robins, JMJ. (1998). Marginal structural models. 1997 proceedings of the American Statistical Association, section on Bayesian statistical science (pp. 1–10). Retrieved from.
  * Robins (2004) [author] Robins, James MJ. M. (2004). Optimal structural nested models for optimal sequential decisions. In Proceedings of the second seattle Symposium in Biostatistics 189–326. Springer.
  * Robins, Hernan and Brumback (2000) [author] Robins, James MJ. M., Hernan, Miguel AngelM. A. and Brumback, BabetteB. (2000). Marginal structural models and causal inference in epidemiology.
  * Robins, Orellana and Rotnitzky (2008) [author] Robins, James M.J. M., Orellana, L.L. and Rotnitzky, A.A. (2008). Estimation and extrapolation of optimal treatment and testing strategies. Statistics in Medicine 27 4678–4721.
  * Robins et al. (2008) [author] Robins, JamesJ., Li, LinglingL., Tchetgen, EricE., van der Vaart, AadA. et al. (2008). Higher order influence functions and minimax estimation of nonlinear functionals. Probability and statistics: essays in honor of David A. Freedman 2 335–421.
  * Robins et al. (2017) [author] Robins, James MJ. M., Li, LinglingL., Mukherjee, RajarshiR., Tchetgen, Eric TchetgenE. T. and van der Vaart, AadA. (2017). Minimax estimation of a functional on a structured high-dimensional model. The Annals of Statistics 45 1951–1987.
  * Rosenblum (2011) [author] Rosenblum, MichaelM. (2011). Marginal structural models. In Targeted Learning 145–160. Springer.
  * Rubin and van der Laan (2012) [author] Rubin, Daniel BD. B. and van der Laan, Mark JM. J. (2012). Statistical issues and limitations in personalized medicine research with clinical trials. The international journal of biostatistics 8 18\.
  * Shinozaki and Suzuki (2020) [author] Shinozaki, TomohiroT. and Suzuki, EtsujiE. (2020). Understanding marginal structural models for time-varying exposures: pitfalls and tips. Journal of epidemiology JE20200226.
  * Swaminathan and Joachims (2015) [author] Swaminathan, AdithA. and Joachims, ThorstenT. (2015). Batch learning from logged bandit feedback through counterfactual risk minimization. The Journal of Machine Learning Research 16 1731–1755.
  * Traskin and Small (2011) [author] Traskin, MikhailM. and Small, Dylan SD. S. (2011). Defining the study population for an observational study to ensure sufficient overlap: a tree approach. Statistics in Biosciences 3 94–118.
  * Tsybakov (2004) [author] Tsybakov, Alexander BA. B. (2004). Optimal aggregation of classifiers in statistical learning. The Annals of Statistics 32 135–166.
  * Vaart and Wellner (1996) [author] Vaart, Aad WA. W. and Wellner, Jon AJ. A. (1996). Weak convergence and empirical processes: with applications to statistics. Springer.
  * van der Laan (2006) [author] van der Laan, Mark JM. J. (2006). Causal effect models for intention to treat and realistic individualized treatment rules.
  * van der Laan (2017) [author] van der Laan, MarkM. (2017). A generally efficient targeted minimum loss based estimator based on the highly adaptive lasso. The international journal of biostatistics 13\.
  * van der Laan and Bibaut (2017) [author] van der Laan, Mark JM. J. and Bibaut, Aurélien FA. F. (2017). Uniform consistency of the highly adaptive lasso estimator of infinite dimensional parameters. arXiv preprint arXiv:1709.06256.
  * van der Laan, Bibaut and Luedtke (2018) [author] van der Laan, Mark JM. J., Bibaut, AurélienA. and Luedtke, Alexander RA. R. (2018). CV-TMLE for nonpathwise differentiable target parameters. In Targeted Learning in Data Science 455–481. Springer.
  * Van der Laan, Laan and Robins (2003) [author] Van der Laan, Mark JM. J., Laan, MJM. and Robins, James MJ. M. (2003). Unified methods for censored longitudinal data and causality. Springer Science & Business Media.
  * van der Laan and Luedtke (2014) [author] van der Laan, Mark JM. J. and Luedtke, Alexander RA. R. (2014). Targeted learning of an optimal dynamic treatment, and statistical inference for its mean outcome. U.C. Berkeley Division of Biostatistics Working Paper Series.
  * van der Laan and Petersen (2007) [author] van der Laan, Mark JM. J. and Petersen, Maya LM. L. (2007). Statistical learning of origin-specific statically optimal individualized treatment rules. The international journal of biostatistics 3\.
  * van Der Laan and Rose (2018) [author] van Der Laan, Mark JM. J. and Rose, SherriS. (2018). Targeted Learning in Data Science: Causal Inference for Complex Longitudinal Studies. Springer.
  * Wahed and Tsiatis (2004) [author] Wahed, Abdus SA. S. and Tsiatis, Anastasios AA. A. (2004). Optimal estimator for the survival distribution and related quantities for treatment policies in two-stage randomization designs in clinical trials. Biometrics 60 124–133.
  * Wasserman (2006) [author] Wasserman, LarryL. (2006). All of nonparametric statistics. Springer Science & Business Media.
  * Zhang et al. (2012) [author] Zhang, BaqunB., Tsiatis, Anastasios AA. A., Laber, Eric BE. B. and Davidian, MarieM. (2012). A robust method for estimating optimal treatment regimes. Biometrics 68 1010–1018.
  * Zhao et al. (2012) [author] Zhao, YingqiY., Zeng, DonglinD., Rush, A JohnA. J. and Kosorok, Michael RM. R. (2012). Estimating individualized treatment rules using outcome weighted learning. Journal of the American Statistical Association 107 1106–1118.
  * Zhao et al. (2015) [author] Zhao, Ying-QiY.-Q., Zeng, DonglinD., Laber, Eric BE. B. and Kosorok, Michael RM. R. (2015). New statistical learning methods for estimating optimal dynamic treatment regimes. Journal of the American Statistical Association 110 583–598.
  * Zhao et al. (2022) [author] Zhao, QingyuanQ., Small, Dylan SD. S., Ertefaie, AshkanA. et al. (2022). Selective inference for effect modification via the lasso. Journal of the Royal Statistical Society Series B 84 382–413.
  * Zheng and van der Laan (2011) [author] Zheng, WenjingW. and van der Laan, Mark JM. J. (2011). Cross-validated targeted minimum-loss-based estimation. In Targeted Learning 459–474. Springer.

Supplement to “Nonparametric estimation of a covariate-adjusted counterfactual
treatment regimen response curve”

##  Appendix S1 Proof of Lemmas A.1, A.2 and A.3

###### Proof of Lemma A.1.

We first proof the result for the treatment indicator and $G^{a}$. Let
$\tilde{f}^{a}=h^{r}f_{h}^{a}$ and
$D^{a}\\{\tilde{f}^{a},G^{a}_{n,\lambda_{n}}\\}=\tilde{f}^{a}\cdot(A-G^{a}_{n,\lambda_{n}})$,
where $\tilde{f}^{a}$ is a càdlàg function with a finite sectional variation
norm, and let $\tilde{f}_{\phi}^{a}$ be an approximation of $\tilde{f}^{a}$
using the basis functions $\phi$ that satisfy condition (21).

$\displaystyle P_{n}D^{a}(f_{h}^{a},G^{a}_{n,\lambda_{n}})$
$\displaystyle=h^{-r}P_{n}D^{a}(\tilde{f}^{a},G^{a}_{n,\lambda_{n}})$
$\displaystyle=h^{-r}P_{n}D^{a}(\tilde{f}^{a}-\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})+h^{-r}P_{n}D^{a}(\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})$
$\displaystyle=h^{-r}(P_{n}-P_{0})D^{a}(\tilde{f}^{a}-\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})+h^{-r}P_{0}D^{a}(\tilde{f}^{a}-\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})+h^{-r}P_{n}D^{a}(\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})$
$\displaystyle=h^{-r}O_{p}(n^{-1/2}n^{-1/6})+h^{-r}O_{p}(n^{-1/3}n^{-1/3})+h^{-r}P_{n}D^{a}(\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}}).$

The last equality follows from
$\|\tilde{f}^{a}-\tilde{f}_{\phi}^{a}\|_{2}=O_{p}(n^{-1/3})$ and Lemma S1 as

$\displaystyle(P_{n}-P_{0})D^{a}(\tilde{f}^{a}-\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})$
$\displaystyle\leq\sup_{\|\tilde{f}^{a}-\tilde{f}_{\phi}^{a}\|_{2}\leq
n^{-1/3}}|(P_{n}-P_{0})D^{a}(\tilde{f}^{a}-\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})|$
$\displaystyle=O_{p}\\{n^{-1/2}\mathcal{E}(n^{-1/3},L^{2}(P))\\}\leq
O_{p}(n^{-1/2}n^{-1/6}).$

Moreover, by the convergence rate of the highly adaptive lasso estimate (i.e.,
$\lVert G^{a}_{0}-G^{a}_{n,\lambda_{n}}\rVert_{2,P_{0}}=O_{p}(n^{-1/3})$),

$P_{0}D^{a}(\tilde{f}^{a}-\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})=P_{0}(\tilde{f}^{a}-\tilde{f}_{\phi}^{a})(G_{0}^{a}-G^{a}_{n,\lambda_{n}})=O_{p}(n^{-1/3})O_{p}(n^{-1/3}).$

Therefore, the proof is complete if we show

$h^{-r}O_{p}(n^{-1/2}n^{-1/6})+h^{-r}O_{p}(n^{-1/3}n^{-1/3})+h^{-r}P_{n}D^{a}(\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})=o_{p}\left((nh^{r})^{-1/2}\right).$

Considering the optimal bandwidth rate of $h^{r}=n^{\frac{-r}{r+2J}}$, the
first term satisfies the desired rate when
$h^{-r}n^{-1/2}n^{-1/6}<n^{-1/2}h^{-r/2}$ which implies
$\frac{-1}{6}<\frac{-r}{2r+4J}$. The inequality is satisfied for $J>r$. We can
similarly show that $h^{-r}O_{p}(n^{-2/3})=o_{p}\left((nh^{r})^{-1/2}\right)$
for $J>r$. Hence, for $J>r$,

$P_{n}D^{a}(f_{h}^{a},G^{a}_{n,\lambda_{n}})=o_{p}\left((nh^{r})^{-1/2}\right)+h^{-r}P_{n}D^{a}(\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}}).$

In the following we show that
$h^{-r}P_{n}D^{a}(\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})$ is also of the
order $o_{p}\left((nh^{r})^{-1/2}\right)$.

For simplicity of notation, let $G^{a\dagger}=\text{logit }G^{a}$. Define a
set of score functions generated by a path $\\{1+\epsilon
g(s,j)\\}\beta_{n,s,j}$ for a uniformly bounded vector $g$ as

$S_{g}(G^{a}_{n,\lambda_{n}})=\frac{d}{dG^{a\dagger}_{n,\lambda_{n}}}L(G^{a\dagger}_{n,\lambda_{n}})\left\\{\sum_{(s,j)}g(s,j)\beta_{n,s,j}\phi_{s,j}\right\\}.$

When $L(\cdot)$ is the log-likelihood loss function,
$L(G)=A\log\left(\frac{G}{1-G}\right)+\log(1-G)$. Thus, $L(G^{a\dagger})=A\log
G^{a\dagger}+\log(G^{a\dagger}-1)$ and
$S_{g}(G^{a})=(A-G^{a}_{n,\lambda_{n}})\left\\{\sum_{(s,j)}g(s,j)\beta_{n,s,j}\phi_{s,j}\right\\}$.
Let $r(g,G^{a}_{n,\lambda_{n}})=\sum_{(s,j)}g(s,j)\lvert\beta_{n,s,j}\rvert$.
For small enough $\epsilon$,

$\displaystyle\sum_{(s,j)}\lvert\\{1+\epsilon g(s,j)\\}\beta_{n,s,j}\rvert$
$\displaystyle=\sum_{(s,j)}\\{1+\epsilon g(s,j)\\}\lvert\beta_{n,s,j}\rvert$
$\displaystyle=\sum_{(s,j)}\lvert\beta_{n,s,j}\rvert+\epsilon
r(g,G_{n,\lambda_{n}}).$

Hence, for any $g$ satisfying $r(g,G_{n,\lambda_{n}})=0$, we have
$P_{n}S_{g}(G^{a}_{n,\lambda_{n}})=0$.

Because
$D^{a}\\{\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}}\\}\in\\{S_{g}(G^{a}_{n,\lambda_{n}}):\lVert
g\rVert_{\infty}<\infty\\}$, there exists $g^{\star}$ such that
$D^{a}(\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})=S_{g^{\star}}(G^{a}_{n,\lambda_{n}})$.
However, for this particular choice of $g^{\star}$,
$r(g,G^{a}_{n,\lambda_{n}})$ may not be zero. Now, define $g$ such that
$g(s,j)=g^{\star}(s,j)$ for $(s,j)\neq(s^{\star},j^{\star})$;
$\tilde{g}(s^{\star},j^{\star})$ is defined such that

$\displaystyle\sum_{(s,j)\neq(s^{\star},j^{\star})}g^{\star}(s,j)\lvert\beta_{n,s,j}\rvert+g(s^{\star},j^{\star})\lvert\beta_{n,s^{\star},j^{\star}}\rvert=0.$
(S29)

That is, $g$ matches $g^{\star}$ everywhere but for a single point
$(s^{\star},j^{\star})$, where it is forced to take a value such that
$r(g,G_{n,\lambda_{n}})=0$. As a result, for such a choice of $g$,
$P_{n}S_{g}(G^{a}_{n,\lambda_{n}})=0$ by definition. Below, we show that
$P_{n}S_{g}(G^{a}_{n,\lambda_{n}})-P_{n}D^{a}\\{\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}}\\}=o_{p}\left(n^{-1/2}h^{r/2}\right)$
which then implies that
$h^{-r}P_{n}D^{a}\\{\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}}\\}=o_{p}\left((nh^{r})^{-1/2}\right)$.
We note that the choice of $(s^{\star},j^{\star})$ is inconsequential.

$\displaystyle
P_{n}S_{g}(G^{a}_{n,\lambda_{n}})-P_{n}D^{a}\\{\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}}\\}$
$\displaystyle=P_{n}S_{g}(G^{a}_{n,\lambda_{n}})-P_{n}S_{g^{*}}(G^{a}_{n,\lambda_{n}})$
$\displaystyle=P_{n}\left\\{\frac{d}{dG^{a\dagger}_{n,\lambda_{n}}}L(G^{a\dagger}_{n,\lambda_{n}})\left[\sum_{(s,j)}\left\\{g(s,j)-g^{\star}(s,j)\right\\}\beta_{n,s,j}\phi_{s,j}\right]\right\\}$
$\displaystyle=P_{n}\left[\frac{d}{dG^{a\dagger}_{n,\lambda_{n}}}L(G^{a\dagger}_{n,\lambda_{n}})\left\\{g(s^{\star},j^{\star})-g^{\star}(s^{\star},j^{\star})\right\\}\beta_{n,s^{\star},j^{\star}}\phi_{s^{\star},j^{\star}}\right]$
$\displaystyle=P_{n}\left[\frac{d}{dG^{a\dagger}_{n,\lambda_{n}}}L(G^{a\dagger}_{n,\lambda_{n}})\kappa(s^{\star},j^{\star})\phi_{s^{\star},j^{\star}}\right]$
$\displaystyle=\kappa(s^{\star},j^{\star})P_{n}\left[\frac{d}{dG^{a\dagger}_{n,\lambda_{n}}}L(G^{a\dagger}_{n,\lambda_{n}})\phi_{s^{\star},j^{\star}}\right],$

where the third equality follows from equation (S30) above with

$\kappa(s^{\star},j^{\star})=-\frac{\sum_{(s,j)\neq(s^{\star},j^{\star})}g^{\star}(s,j)\lvert\beta_{n,s,j}\rvert}{\lvert\beta_{n,s^{\star},j^{\star}}\rvert}\beta_{n,s^{\star},j^{\star}}-g^{\star}(s^{\star},j^{\star})\beta_{n,s^{\star},j^{\star}}.$

Moreover,

$\left\lvert\kappa(s^{\star},j^{\star})\right\rvert\leq\sum_{(s,j)}\lvert
g^{\star}(s,j)\beta_{n,s,j}\rvert.$

Assuming $\tilde{f}_{\phi}^{a}$ has finite sectional variation norm, the
$L_{1}$ norm of the coefficients approximating $\tilde{f}_{\phi}^{a}$ will be
finite which implies that $\sum_{(s,j)}\lvert
g^{\star}(s,j)\beta_{n,s,j}\rvert$ is finite, and thus
$\lvert\kappa(s^{\star},j^{\star})\rvert=O_{p}(1)$. Then,

$\displaystyle
P_{n}S_{g}(G^{a}_{n,\lambda_{n}})-P_{n}D^{a}(\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})$
$\displaystyle=O_{p}\left(P_{n}\left[\frac{d}{dG^{a\dagger}_{n,\lambda_{n}}}L(G^{a\dagger}_{n,\lambda_{n}})\phi_{s^{\star},j^{\star}}\right]\right)$
$\displaystyle=o_{p}\left(n^{-1/2}h^{r/2}\right),$

where the last equality follows from the assumption that
$\min_{(s,j)\in\mathcal{J}_{n}}\lVert
P_{n}\frac{d}{dG^{a\dagger}_{n,\lambda_{n}}}L(G^{a\dagger}_{n,\lambda_{n}})(\phi_{s,j})\rVert=o_{p}\left(n^{-1/2}h^{r/2}\right)$
for $L(\cdot)$ being log-likelihood loss. As
$P_{n}S_{{g}}(G^{a}_{n,\lambda_{n}})=0$, it follows that
$P_{n}D^{a}(\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})=o_{p}\left(n^{-1/2}h^{r/2}\right)$
and thus
$h^{-r}P_{n}D^{a}(\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})=o_{p}\left((nh^{r})^{-1/2}\right)$.
We just showed that
$P_{n}D^{a}(f_{h}^{a},G^{a}_{n,\lambda_{n}})=o_{p}\left((nh^{r})^{-1/2}\right)$.
Similarly we can show that
$P_{n}D^{c}(f_{h}^{c},G^{c}_{n,\lambda_{n}})=o_{p}\left((nh^{r})^{-1/2}\right)$
and $P_{n}D^{\Psi}(f_{h}^{\Psi},\Psi_{n})=o_{p}\left((nh^{r})^{-1/2}\right)$
which completes the proof. ∎

###### Proof of Lemma A.2.

The prove follows from the proof of Lemma A.1 with the bandwidth $h$ replaced
by a constant value (e.g., one). ∎

###### Proof of Lemma A.3.

We first proof the result for the treatment indicator and $G^{a}$.

$\displaystyle\sup_{v\in\mathcal{V}}\left|P_{n}D^{a}(f^{a},G^{a}_{n,\lambda_{n}})\right|$
$\displaystyle\leq\sup_{v\in\mathcal{V}}\left|P_{n}D^{a}(f^{a}-\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})\right|+\sup_{v\in\mathcal{V}}\left|P_{n}D^{a}(\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})\right|$
$\displaystyle\leq\sup_{v\in\mathcal{V}}\left|(P_{n}-P_{0})D^{a}(f^{a}-\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})\right|+\sup_{v\in\mathcal{V}}\left|P_{0}D^{a}(f^{a}-\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})\right|$
$\displaystyle+\sup_{v\in\mathcal{V}}\left|P_{n}D^{a}(\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})\right|$
$\displaystyle=O_{p}(n^{-1/2}n^{-1/6})+O_{p}(n^{-1/3}n^{-1/3})+\sup_{v\in\mathcal{V}}\left|P_{n}D^{a}(\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})\right|.$

The last equality follows from
$\sup_{v\in\mathcal{V}}\|\tilde{f}^{a}-\tilde{f}_{\phi}^{a}\|_{2}=O_{p}(n^{-1/3})$
and Lemma S1 as

$\displaystyle\sup_{v\in\mathcal{V}}(P_{n}-P_{0})D^{a}(\tilde{f}^{a}-\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})$
$\displaystyle\leq\sup_{\begin{subarray}{c}v\in\mathcal{V}\\\
\|\tilde{f}^{a}-\tilde{f}_{\phi}^{a}\|_{2}\leq
n^{-1/3}\end{subarray}}|(P_{n}-P_{0})D^{a}(\tilde{f}^{a}-\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})|$
$\displaystyle=O_{p}\\{n^{-1/2}\mathcal{E}(n^{-1/3},L^{2}(P))\\}\leq
O_{p}(n^{-1/2}n^{-1/6}).$

Moreover, by the convergence rate of the highly adaptive lasso estimate (i.e.,
$\sup_{v\in\mathcal{V}}\lVert
G^{a}_{0}-G^{a}_{n,\lambda_{n}}\rVert_{2,P_{0}}=O_{p}(n^{-1/3})$),

$\sup_{v\in\mathcal{V}}\left|P_{0}D^{a}(\tilde{f}^{a}-\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})\right|=\sup_{v\in\mathcal{V}}\left|P_{0}(\tilde{f}^{a}-\tilde{f}_{\phi}^{a})(G_{0}^{a}-G^{a}_{n,\lambda_{n}})\right|=O_{p}(n^{-1/3})O_{p}(n^{-1/3}).$

Therefore, the proof is complete if we show

$O_{p}(n^{-1/2}n^{-1/6})+O_{p}(n^{-1/3}n^{-1/3})+\sup_{v\in\mathcal{V}}\left|P_{n}D^{a}(\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})\right|=o_{p}(n^{-1/2}).$

$\sup_{v\in\mathcal{V}}\left|P_{n}D^{a}(f^{a},G^{a}_{n,\lambda_{n}})\right|=o_{p}(n^{-1/2})+\sup_{v\in\mathcal{V}}\left|P_{n}D^{a}(\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})\right|.$

In the following we show that
$\sup_{v\in\mathcal{V}}\left|P_{n}D^{a}(\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})\right|$
is also of the order $o_{p}(n^{-1/2})$.

For simplicity of notation, let $G^{a\dagger}=\text{logit }G^{a}$. Define a
set of score functions generated by a path $\\{1+\epsilon
g(s,j)\\}\beta_{n,s,j}$ for a uniformly bounded vector $g$ as

$S_{g}(G^{a}_{n,\lambda_{n}})=\frac{d}{dG^{a\dagger}_{n,\lambda_{n}}}L(G^{a\dagger}_{n,\lambda_{n}})\left\\{\sum_{(s,j)}g(s,j)\beta_{n,s,j}\phi_{s,j}\right\\}.$

When $L(\cdot)$ is the log-likelihood loss function,
$L(G)=A\log\left(\frac{G}{1-G}\right)+\log(1-G)$. Thus, $L(G^{a\dagger})=A\log
G^{a\dagger}+\log(G^{a\dagger}-1)$ and
$S_{g}(G^{a})=(A-G^{a}_{n,\lambda_{n}})\left\\{\sum_{(s,j)}g(s,j)\beta_{n,s,j}\phi_{s,j}\right\\}$.
Let $r(g,G^{a}_{n,\lambda_{n}})=\sum_{(s,j)}g(s,j)\lvert\beta_{n,s,j}\rvert$.
For small enough $\epsilon$,

$\displaystyle\sum_{(s,j)}\lvert\\{1+\epsilon g(s,j)\\}\beta_{n,s,j}\rvert$
$\displaystyle=\sum_{(s,j)}\\{1+\epsilon g(s,j)\\}\lvert\beta_{n,s,j}\rvert$
$\displaystyle=\sum_{(s,j)}\lvert\beta_{n,s,j}\rvert+\epsilon
r(g,G_{n,\lambda_{n}}).$

Hence, for any $g$ satisfying $r(g,G_{n,\lambda_{n}})=0$, we have
$P_{n}S_{g}(G^{a}_{n,\lambda_{n}})=0$.

Because
$D^{a}\\{\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}}\\}\in\\{S_{g}(G^{a}_{n,\lambda_{n}}):\lVert
g\rVert_{\infty}<\infty\\}$, there exists $g^{\star}$ such that
$D^{a}(\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})=S_{g^{\star}}(G^{a}_{n,\lambda_{n}})$.
However, for this particular choice of $g^{\star}$,
$r(g,G^{a}_{n,\lambda_{n}})$ may not be zero. Now, define $g$ such that
$g(s,j)=g^{\star}(s,j)$ for $(s,j)\neq(s^{\star},j^{\star})$;
$\tilde{g}(s^{\star},j^{\star})$ is defined such that

$\displaystyle\sum_{(s,j)\neq(s^{\star},j^{\star})}g^{\star}(s,j)\lvert\beta_{n,s,j}\rvert+g(s^{\star},j^{\star})\lvert\beta_{n,s^{\star},j^{\star}}\rvert=0.$
(S30)

That is, $g$ matches $g^{\star}$ everywhere but for a single point
$(s^{\star},j^{\star})$, where it is forced to take a value such that
$r(g,G_{n,\lambda_{n}})=0$. As a result, for such a choice of $g$,
$P_{n}S_{g}(G^{a}_{n,\lambda_{n}})=0$ by definition. Below, we show that
$P_{n}S_{g}(G^{a}_{n,\lambda_{n}})-P_{n}D^{a}\\{\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}}\\}=o_{p}(n^{-1/2})$
which then implies that
$P_{n}D^{a}\\{\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}}\\}=o_{p}(n^{-1/2})$.
We note that the choice of $(s^{\star},j^{\star})$ is inconsequential.

$\displaystyle
P_{n}S_{g}(G^{a}_{n,\lambda_{n}})-P_{n}D^{a}\\{\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}}\\}$
$\displaystyle=P_{n}S_{g}(G^{a}_{n,\lambda_{n}})-P_{n}S_{g^{*}}(G^{a}_{n,\lambda_{n}})$
$\displaystyle=P_{n}\left\\{\frac{d}{dG^{a\dagger}_{n,\lambda_{n}}}L(G^{a\dagger}_{n,\lambda_{n}})\left[\sum_{(s,j)}\left\\{g(s,j)-g^{\star}(s,j)\right\\}\beta_{n,s,j}\phi_{s,j}\right]\right\\}$
$\displaystyle=P_{n}\left[\frac{d}{dG^{a\dagger}_{n,\lambda_{n}}}L(G^{a\dagger}_{n,\lambda_{n}})\left\\{g(s^{\star},j^{\star})-g^{\star}(s^{\star},j^{\star})\right\\}\beta_{n,s^{\star},j^{\star}}\phi_{s^{\star},j^{\star}}\right]$
$\displaystyle=P_{n}\left[\frac{d}{dG^{a\dagger}_{n,\lambda_{n}}}L(G^{a\dagger}_{n,\lambda_{n}})\kappa(s^{\star},j^{\star})\phi_{s^{\star},j^{\star}}\right]$
$\displaystyle=\kappa(s^{\star},j^{\star})P_{n}\left[\frac{d}{dG^{a\dagger}_{n,\lambda_{n}}}L(G^{a\dagger}_{n,\lambda_{n}})\phi_{s^{\star},j^{\star}}\right],$

where the third equality follows from equation (S30) above with

$\kappa(s^{\star},j^{\star})=-\frac{\sum_{(s,j)\neq(s^{\star},j^{\star})}g^{\star}(s,j)\lvert\beta_{n,s,j}\rvert}{\lvert\beta_{n,s^{\star},j^{\star}}\rvert}\beta_{n,s^{\star},j^{\star}}-g^{\star}(s^{\star},j^{\star})\beta_{n,s^{\star},j^{\star}}.$

Moreover,

$\left\lvert\kappa(s^{\star},j^{\star})\right\rvert\leq\sum_{(s,j)}\lvert
g^{\star}(s,j)\beta_{n,s,j}\rvert.$

Assuming $\tilde{f}_{\phi}^{a}$ has finite sectional variation norm, the
$L_{1}$ norm of the coefficients approximating $\tilde{f}_{\phi}^{a}$ will be
finite which implies that $\sum_{(s,j)}\lvert
g^{\star}(s,j)\beta_{n,s,j}\rvert$ is finite, and thus
$\lvert\kappa(s^{\star},j^{\star})\rvert=O_{p}(1)$. Hence,

$\displaystyle\sup_{v\in\mathcal{V}}\left|P_{n}S_{g}(G^{a}_{n,\lambda_{n}})-P_{n}D^{a}(\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})\right|$
$\displaystyle=O_{p}\left(\sup_{v\in\mathcal{V}}\left|P_{n}\left[\frac{d}{dG^{a\dagger}_{n,\lambda_{n}}}L(G^{a\dagger}_{n,\lambda_{n}})\phi_{s^{\star},j^{\star}}\right]\right|\right)$
$\displaystyle=o_{p}(n^{-1/2}),$

where the last equality follows from the assumption that
$\min_{(s,j)\in\mathcal{J}_{n}}\sup_{v\in\mathcal{V}}\lVert
P_{n}\frac{d}{dG^{a\dagger}_{n,\lambda_{n}}}L(G^{a\dagger}_{n,\lambda_{n}})(\phi_{s,j})\rVert=o_{p}(n^{-1/2})$
for $L(\cdot)$ being log-likelihood loss. As
$P_{n}S_{{g}}(G^{a}_{n,\lambda_{n}})=0$, it follows that
$P_{n}D^{a}(\tilde{f}_{\phi}^{a},G^{a}_{n,\lambda_{n}})=o_{p}(n^{-1/2})$.
We just showed that $P_{n}D^{a}(f^{a},G^{a}_{n,\lambda_{n}})=o_{p}(n^{-1/2})$.
Similarly we can show that
$P_{n}D^{c}(f^{c},G^{c}_{n,\lambda_{n}})=o_{p}(n^{-1/2})$ and
$P_{n}D^{\Psi}(f^{\Psi},\Psi_{n})=o_{p}(n^{-1/2})$ which completes the proof.
∎

##  Appendix S2 Other Lemmas

######  Lemma S1 (Lemma 3.4.2 in Vaart and Wellner (1996)).

For a class of functions $\mathcal{F}$ with envelop
$\textbf{F}=\sup_{f\in\mathcal{F}}|f(x)|$ bounded from above by $M<\infty$, we
have

$E_{P}\sup_{f\in\mathcal{F},\|f\|_{P}\leq\delta}n^{1/2}(P_{n}-P_{0})(f)\leq\mathcal{E}(\delta,\mathcal{F},L^{2}(P))\left\\{1+\frac{\mathcal{E}(\delta,\mathcal{F},L^{2}(P))}{\delta^{2}n^{1/2}}M\right\\}.$

This implies that
$\sup_{f\in\mathcal{F},\|f\|_{P}\leq\delta}n^{1/2}(P_{n}-P_{0})(f)$ is bounded
in probability by the right-hand side.

######  Lemma S2 (Proposition 2 in Bibaut and van der Laan (2019)).

Let $p\geq 2$ and $M>0$. Denote $\mathcal{F}_{p,M}$ the class of cadlag
functions on $[0,1]^{p}$ with sectional variation norm smaller than $M$.
Suppose that $P_{0}$ is such that, for all $1\leq r\leq\infty$, for all real-
valued function $f$ on $[0,1]^{p}$, $\|f\|_{P_{0},b}=c(b)\|f\|_{\mu,b}$, for
some $c(b)>0$, and where $\mu$ is the Lebesgue measure. Then, for any
$0<\delta<1$, the bracketing entropy integral of $\mathcal{F}_{p,M}$ with
respect to $\|.\|_{P_{0},b}$ norm is bounded by

$\mathcal{E}_{[]}(\delta,\mathcal{F}_{p,M},\|.\|_{P_{0},b})\leq\\{C(b,p)M\delta\\}^{1/2}\mid\log(\delta/M)\mid^{p-1},$

where $C(b,p)$ is a constant that depends only on $b$ and $p$.

######  Lemma S3.

Consider a sequence of processes $\\{n^{1/2}(P_{n}-P_{0})H_{n,h}(f)\\}$ where

$H_{n,h}(f)=h^{r/2}K_{h,\bm{v}_{0}}(f_{n}-f_{0}).$

Assuming the following:

  * (i)

The function $f_{0}$ belongs to a cadlag class of functions with finite
sectional variation norm.

  * (ii)

The function $f_{n}$ is the highly adaptive lasso estimate of $f_{0}$.

  * (iii)

The bandwidth $h_{n}$ satisfies $h^{r}\rightarrow 0$ and
$nh^{r3/2}\rightarrow\infty$.

Then $n^{1/2}(P_{n}-P_{0})H_{n,h}(f_{n})=o_{p}(1)$.

###### Proof.

Let $\tilde{K}_{h,\bm{v}_{0}}(v)=h^{r}K_{h,\bm{v}_{0}}(v)$. Define
$h^{-r/2}\tilde{H}_{n,h}(f)=H_{n,h}(f)$. The function $\tilde{H}_{n}(f)$
belongs to a class of function $\mathcal{F}_{n}=\\{\tilde{H}_{n,h}(f):h\\}$.
We have

$\displaystyle\|\tilde{H}_{n}(f_{n})\|_{P_{0}}$
$\displaystyle\leq\|\tilde{K}_{h,\bm{v}_{0}}(v)\|_{P_{0}}\|(f_{n}-f_{0})\|_{P_{0}}$
$\displaystyle=O_{p}(h^{r/2})O_{p}(n^{-1/3})$

where the last equality follows from the kernel and the highly adaptive lasso
rate of convergence. Using the result of Lemma S2,
$n^{1/2}(P_{n}-P_{0})\tilde{H}_{n}(f_{n})=O_{p}(h^{r/4}n^{-1/6})$. Therefore,

$\displaystyle n^{1/2}(P_{n}-P_{0})H_{n}(f_{n})=O_{p}(h^{-r/4}n^{-1/6}).$

Thus, to obtain the desired rate we must have $h^{-r/4}n^{-1/6}=o(1)$ which
implies that $h^{r}$ has to converge to zero slower than $n^{-2/3}$ (i.e.,
$h^{r}>n^{-2/3}$). Note that under continuity and no smoothness assumptions,
the rate for the bandwidth is $h^{r}=n^{-1/3}$ which is much slower than
$n^{-2/3}$. Thus, under continuity, for every level of smoothness the choice
of optimal $h$ will satisfy the rate $h^{r}=n^{-2/3}$ (i.e., no condition is
imposed).

∎

######  Lemma S4.

Let the function $f_{n}$ is the highly adaptive lasso estimate of $f_{0}$
where $f_{0}$ belongs to a cadlag class of functions with finite sectional
variation norm. Then

$P_{n}\frac{K_{h,\bm{v}_{0}}f_{0}}{P_{0}K_{h,\bm{v}_{0}}}\left(f_{0}-f_{n}\right)=P_{n}\frac{K_{h,\bm{v}_{0}}f_{n}}{P_{n}K_{h,\bm{v}_{0}}}\left(f_{0}-f_{n}\right)+o_{p}\left((nh^{r})^{-1/2}\right).$

###### Proof.

We have

$\displaystyle
P_{n}\frac{K_{h,\bm{v}_{0}}f_{0}}{P_{0}K_{h,\bm{v}_{0}}}\left(f_{0}-f_{n}\right)$
$\displaystyle=P_{n}\frac{K_{h,\bm{v}_{0}}f_{n}}{P_{n}K_{h,\bm{v}_{0}}}\left(f_{0}-f_{n}\right)+P_{n}K_{h,\bm{v}_{0}}(f_{0}-f_{n})\left(\frac{f_{0}}{P_{0}K_{h,\bm{v}_{0}}}-\frac{f_{n}}{P_{n}K_{h,\bm{v}_{0}}}\right)$
$\displaystyle=P_{n}\frac{K_{h,\bm{v}_{0}}f_{n}}{P_{n}K_{h,\bm{v}_{0}}}\left(f_{0}-f_{n}\right)+(P_{n}-P_{0})K_{h,\bm{v}_{0}}(f_{0}-f_{n})\left(\frac{f_{0}}{P_{0}K_{h,\bm{v}_{0}}}-\frac{f_{n}}{P_{n}K_{h,\bm{v}_{0}}}\right)$
$\displaystyle\hskip
144.54pt+P_{0}K_{h,\bm{v}_{0}}(f_{0}-f_{n})\left(\frac{f_{0}}{P_{0}K_{h,\bm{v}_{0}}}-\frac{f_{n}}{P_{n}K_{h,\bm{v}_{0}}}\right)$
$\displaystyle=P_{n}\frac{K_{h,\bm{v}_{0}}f_{n}}{P_{n}K_{h,\bm{v}_{0}}}\left(f_{0}-f_{n}\right)+o_{p}\left((nh^{r})^{-1/2}\right).$

The last equality follows because
$P_{0}K_{h,\bm{v}_{0}}(f_{0}-f_{n})\left(\frac{f_{0}}{P_{0}K_{h,\bm{v}_{0}}}-\frac{f_{n}}{P_{n}K_{h,\bm{v}_{0}}}\right)=o_{p}\left((nh^{r})^{-1/2}\right)$
by the highly adaptive lasso rate of convergence. ∎

######  Lemma S5.

Let $P_{0}L_{G_{0}}(\Psi_{n},\Psi_{0})\equiv
P_{0}L_{G_{0}}(\Psi_{n})-P_{0}L_{G_{0}}(\Psi_{0})$ where $L_{G}$ is the
weighted squared-loss function defined in (1) and
$\mathcal{M}_{\lambda}=\\{m\in\mathcal{M}:\|m\|_{\nu}<\lambda\\}$ be the set
of cadlag functions with variation norm smaller than $\lambda$. Then,
uniformly over $\lambda$, we have

$\left|P_{0}\\{L_{G_{0}}(\Psi_{n},\Psi_{0})-L_{G_{n}}(\Psi_{n},\Psi_{0})\\}\right|\leq
C\left\\{d_{01}(G_{n},G_{0})\right\\}^{1/2}\left\\{P_{0}L_{G_{0}}(\Psi_{n},\Psi_{0})\right\\}^{1/2},$

where $C$ is a positive finite constant.

######  Lemma S6.

Suppose $\|\theta-\theta_{0}\|_{2}=O_{p}(\delta_{n})$ and
$f(\theta)=\frac{h^{r}K_{h,\bm{v}_{0}}I(A=d^{\theta})}{G_{0}(\theta)}Q_{0}(\theta)$.
Under Assumption 4,
$\|f(\theta)-f(\theta_{0})\|_{2}=O_{p}(h^{r/2}\delta_{n})$.

###### Proof.

We assume that $\bm{S}$ and $Q_{0}(\theta)$ are bounded uniformly by some
constants $C_{S}$ and $C_{Q}$. Moreover, by strong positivity assumption
$\min\\{G_{0}(\theta,\bm{W}),1-G_{0}(\theta,\bm{W})\\}>\gamma$ for all
$\bm{W}\in\mathcal{W}$ and $\theta\in\Theta$. By definition the norm can be
written as

$\displaystyle\|f(\theta)-f(\theta_{0})\|_{2}$
$\displaystyle=\left\|h^{r}K_{h,\bm{v}_{0}}\left\\{\frac{I(A=d^{\theta})}{G_{0}(\theta)}Q_{0}(\theta)-\frac{I(A=d^{\theta_{0}})}{G_{0}(\theta_{0})}Q_{0}(\theta_{0})\right\\}\right\|_{2}$
$\displaystyle\leq\|h^{r}K_{h,\bm{v}_{0}}\|_{2}\left\|(\gamma)^{-1}\max\\{|Q_{0}(\theta_{0})|,|Q_{0}(\theta)|\\}I(d^{\theta}\neq
d^{\theta_{0}})\right\|_{2}$
$\displaystyle=\|h^{r}K_{h,\bm{v}_{0}}\|_{2}\left\|(\gamma)^{-1}\max\\{|Q_{0}(\theta_{0})|,|Q_{0}(\theta)|\\}I(d^{\theta}\neq
d^{\theta_{0}})\\{I(|\bm{S}^{\top}\theta_{0}|>t)+I(|\bm{S}^{\top}\theta_{0}|<t)\\}\right\|_{2}$
$\displaystyle\leq\gamma^{-1}C_{Q}\|h^{r}K_{h,\bm{v}_{0}}\|_{2}\left\\{Ct^{\kappa/2}+\left\|I(|\bm{S}^{\top}\theta_{0}-\bm{S}^{\top}\theta|>t)\right\|_{2}\right\\}$
$\displaystyle\leq\gamma^{-1}C_{Q}\|h^{r}K_{h,\bm{v}_{0}}\|_{2}\left\\{Ct^{\kappa/2}+\frac{\|\bm{S}^{\top}\theta_{0}-\bm{S}^{\top}\theta\|_{2}}{t}\right\\}.$

The forth inequality follows from Assumption 4, the Cauchy-Schwarz and
$|\bm{S}^{\top}\theta_{0}-\bm{S}^{\top}\theta|>|\bm{S}^{\top}\theta_{0}|$. The
latter holds because because $I(d^{\theta}\neq d^{\theta_{0}})=1$ when (a)
$\bm{S}^{\top}{{\theta}}>0$ and $\bm{S}^{\top}\theta_{0}<0$; or (b)
$\bm{S}^{\top}\theta<0$ and $\bm{S}^{\top}\theta_{0}>0$. The former and latter
imply that
$\bm{S}^{\top}{\theta}_{0}-\bm{S}^{\top}\theta<\bm{S}^{\top}\theta_{0}<0$ and
$0<\bm{S}\theta_{0}<\bm{S}^{\top}\theta_{0}-\bm{S}^{\top}\theta$,
respectively, and thus,
$|\bm{S}^{\top}\theta_{0}-\bm{S}^{\top}\theta|>|\bm{S}^{\top}\theta_{0}|$. The
fifth inequality holds by Markov inequalities. The upper bound is minimized
when
$t=C_{t}\|\bm{S}^{\top}\theta_{0}-\bm{S}^{\top}\theta\|_{2}^{\frac{1}{\kappa/2+1}}$
for a constant $C_{t}$ which depends on $\kappa$ and $C$. Hence,

$\displaystyle\|f(\theta)-f(\theta_{0})\|_{2}$ $\displaystyle\leq
C^{\dagger}\|h^{r}K_{h,\bm{v}_{0}}\|_{2}\left\\{C\|\bm{S}^{\top}\theta_{0}-\bm{S}^{\top}\theta\|_{2}^{\frac{\kappa}{\kappa+2}}+\|\bm{S}^{\top}\theta_{0}-\bm{S}^{\top}\theta\|_{2}^{\frac{\kappa}{\kappa+2}}\right\\}$
$\displaystyle=O_{p}(h^{r/2}\delta^{\frac{\kappa}{\kappa+2}}).$ (S31)

Hence where there is a margin around zero, that is,
$pr(0<|\bm{S}^{\top}{{\theta}}_{0}|<l)=0$, the right hand side of (S2) reduces
to $O_{p}(h^{r/2}\delta)$. ∎

##  Appendix S3 Proof of Theorem 4.1: Rate of convergence of HAL with unknown
nuisance parameters

Define $P_{0}L_{G_{0}}(\Psi_{n},\Psi_{0})\equiv
P_{0}L_{G_{0}}(\Psi_{n})-P_{0}L_{G_{0}}(\Psi_{0})$. Then using the definition
of the loss-based dissimilarity, we have

$\displaystyle 0\leq d_{0}(\Psi_{n},\Psi_{0})$
$\displaystyle=P_{0}L_{G_{0}}(\Psi_{n},\Psi_{0})$
$\displaystyle=(P_{0}-P_{n})L_{G_{0}}(\Psi_{n},\Psi_{0})+P_{n}L_{G_{0}}(\Psi_{n},\Psi_{0})$
$\displaystyle=O_{p}(n^{-2/3}(\log
n)^{4(p-1)/3})+P_{n}L_{G_{n}}(\Psi_{n},\Psi_{0})+P_{n}\\{L_{G_{0}}(\Psi_{n},\Psi_{0})-L_{G_{n}}(\Psi_{n},\Psi_{0})\\}$
$\displaystyle\leq O_{p}(n^{-2/3}(\log
n)^{4(p-1)/3})+(P_{n}-P_{0})\\{L_{G_{0}}(\Psi_{n},\Psi_{0})-L_{G_{n}}(\Psi_{n},\Psi_{0})\\}$
$\displaystyle\hskip
216.81pt+P_{0}\\{L_{G_{0}}(\Psi_{n},\Psi_{0})-L_{G_{n}}(\Psi_{n},\Psi_{0})\\}$
$\displaystyle=P_{0}\\{L_{G_{0}}(\Psi_{n},\Psi_{0})-L_{G_{n}}(\Psi_{n},\Psi_{0})\\}+O_{p}(n^{-2/3}(\log
n)^{4(p-1)/3})$

Moreover, using Lemma S5 we have

$\displaystyle d_{0}(\Psi_{n},\Psi_{0})$ $\displaystyle\leq
C\left\\{d_{01}(G_{n},G_{0})\right\\}^{1/2}\left\\{d_{0}(\Psi_{n},\Psi_{0})\right\\}^{1/2}+O_{p}(n^{-2/3}(\log
n)^{4(p-1)/3})$ $\displaystyle=O_{p}(n^{-2/3}(\log
n)^{4(p-1)/3})+O_{p}(n^{-1/3}(\log
n)^{2(p-1)/3})\left\\{d_{0}(\Psi_{n},\Psi_{0})\right\\}^{1/2}$

where $C$ is a positive finite constant. The last equality follows from the
convergence rate of HAL minimum loss based estimator of $g_{0}$. Let
$r_{1}(n)=O_{p}(n^{-2/3}(\log n)^{4(p-1)/3})$ and
$r_{2}(n)=O_{p}(n^{-1/3}(\log n)^{2(p-1)/3})$, the above equality is
equivalent to

$\left\\{d_{0}(\Psi_{n},\Psi_{0})\right\\}^{1/2}\leq\frac{\\{4r_{1}(n)+r_{2}^{2}(n)\\}^{1/2}+r_{2}(n)}{2}$

which implies that $d_{0}(\Psi_{n},\Psi_{0})=O_{p}(n^{-2/3}(\log
n)^{4(p-1)/3})$.

##  Appendix S4 Proof of Theorem 4.2: Asymptotic linearity of the regimen-
response curve estimator

We represent the difference between $\Psi_{nh}(\bm{v}_{0},\theta)$ and
$\Psi_{0}(\bm{v}_{0},\theta)$ as

$\Psi_{nh}(\bm{v}_{0},\theta)-\Psi_{0}(\bm{v}_{0},\theta)=\\{\Psi_{0h}(\bm{v}_{0},\theta)-\Psi_{0}(\bm{v}_{0},\theta)\\}+\\{\Psi_{nh}(\bm{v}_{0},\theta)-\Psi_{0h}(\bm{v}_{0},\theta)\\}.$

Moreover, we can write

$\displaystyle\Psi_{nh}(\bm{v}_{0},\theta)-\Psi_{0h}(\bm{v}_{0},\theta)=$
$\displaystyle
P_{n}K_{h,\bm{v}_{0}}\Psi_{n}\left(\frac{1}{P_{n}K_{h,\bm{v}_{0}}}-\frac{1}{P_{0}K_{h,\bm{v}_{0}}}\right)$
(S32)
$\displaystyle+(P_{n}-P_{0})\frac{K_{h,\bm{v}_{0}}}{P_{0}K_{h,\bm{v}_{0}}}(\Psi_{n}-\Psi_{0})+(P_{n}-P_{0})\frac{K_{h,\bm{v}_{0}}}{P_{0}K_{h,\bm{v}_{0}}}(\Psi_{0})$
$\displaystyle+P_{0}\frac{K_{h,\bm{v}_{0}}}{P_{0}K_{h,\bm{v}_{0}}}(\Psi_{n}-\Psi_{0}).$

###  S4.1 Convergence rate of $\Psi_{0h}(v_{0},\theta)-\Psi_{0}(v_{0},\theta)$

We consider a ($J-1$)-orthogonal kernel with bandwidth $h$ centered at
$\bm{v}_{0}$. Then following Lemma 25.1 in van Der Laan and Rose (2018), as
$h\rightarrow 0$,

$\Psi_{0h}(\bm{v}_{0},\theta)-\Psi_{0}(\bm{v}_{0},\theta)=h^{J}B_{0}(J,\bm{v}_{0}),$

where

$B_{0}(J,\bm{v}_{0})=\sum_{\\{\eta\in\\{0,\cdots,J\\}^{d}:\sum_{l}\eta_{l}=J\\}}\Psi_{0}(\bm{v}_{0})^{\eta}\int_{s}k(s)\frac{\prod_{l}s_{l}^{\eta_{l}}}{\prod_{l}\eta_{l}!}ds.$

###  S4.2 Asymptotic behaviour of the first term on the RHS of (S32)

Because, $(P_{n}-P_{0})K_{h,\bm{v}_{0}}=O_{p}(n^{-1/2}h^{-r/2})$,

$\displaystyle
P_{n}K_{h,\bm{v}_{0}}\Psi_{n}\left(\frac{1}{P_{n}K_{h,\bm{v}_{0}}}-\frac{1}{P_{0}K_{h,\bm{v}_{0}}}\right)=$
$\displaystyle(P_{n}-P_{0})K_{h,\bm{v}_{0}}\Psi_{n}\left(\frac{1}{P_{n}K_{h,\bm{v}_{0}}}-\frac{1}{P_{0}K_{h,\bm{v}_{0}}}\right)$
$\displaystyle+P_{0}K_{h,\bm{v}_{0}}\Psi_{n}\left(\frac{1}{P_{n}K_{h,\bm{v}_{0}}}-\frac{1}{P_{0}K_{h,\bm{v}_{0}}}\right)$
$\displaystyle=$ $\displaystyle
P_{0}K_{h,\bm{v}_{0}}\Psi_{n}\left(\frac{-(P_{n}-P_{0})K_{h,\bm{v}_{0}}}{P_{n}K_{h,\bm{v}_{0}}P_{0}K_{h,\bm{v}_{0}}}\right)+o_{p}(n^{-1/2}h^{-r/2})$
$\displaystyle=$ $\displaystyle
P_{0}K_{h,\bm{v}_{0}}\Psi_{n}\left(\frac{-(P_{n}-P_{0})K_{h,\bm{v}_{0}}}{P_{0}^{2}K_{h,\bm{v}_{0}}}\right)$
$\displaystyle+P_{0}K_{h,\bm{v}_{0}}\Psi_{n}\left(\frac{-\\{(P_{n}-P_{0})K_{h,\bm{v}_{0}}\\}^{2}}{P_{n}K_{h,\bm{v}_{0}}P_{0}K_{h,\bm{v}_{0}}}\right)+o_{p}(n^{-1/2}h^{-r/2})$
$\displaystyle=\frac{-(P_{n}-P_{0})K_{h,\bm{v}_{0}}}{P_{0}K_{h,\bm{v}_{0}}}P_{0}\frac{K_{h,\bm{v}_{0}}}{P_{0}K_{h,\bm{v}_{0}}}(\Psi_{n}-\Psi_{0})+\frac{-(P_{n}-P_{0})K_{h,\bm{v}_{0}}}{P_{0}K_{h,\bm{v}_{0}}}\Psi_{0h}$
$\displaystyle+o_{p}(n^{-1/2}h^{-r/2})$
$\displaystyle=\frac{-(P_{n}-P_{0})K_{h,\bm{v}_{0}}}{P_{0}K_{h,\bm{v}_{0}}}\Psi_{0h}+o_{p}(n^{-1/2}h^{-r/2})$

The second and third equalities follows because
$\\{(P_{n}-P_{0})K_{h,\bm{v}_{0}}\\}^{2}=O_{p}(n^{-1}h^{-r})$, and thus, those
are of order $o_{p}(n^{-1/2}h^{-r/2})$. In the last equality we have also used
the consistency of $\Psi_{n}$.

###  S4.3 Asymptotic negligibility of the second term on the RHS of (S32)

We will show that the second order term
$(P_{n}-P_{0})K_{h,\bm{v}_{0}}(\Psi_{n}-\Psi_{0})=o_{p}\left((nh^{r})^{-1/2}\right)$.
Let $\tilde{K}_{h,\bm{v}_{0}}(v)=h^{r}K_{h,\bm{v}_{0}}(v)$. Define
$H_{n,h}(m)=K_{h,\bm{v}_{0}}(m-\Psi_{0})$ and
$h^{-r}\tilde{H}_{n,h}(m)=H_{n,h}(m)$. The function $\tilde{H}_{n}(m)$ belongs
to a class of function $\mathcal{F}_{n}=\\{\tilde{H}_{n,h}(m):h\\}$. We have

$\displaystyle\|\tilde{H}_{n}(\Psi_{n})\|_{P_{0}}$
$\displaystyle\leq\|\tilde{K}_{h,\bm{v}_{0}}(v)\|_{P_{0}}\|(\Psi_{n}-\Psi_{0})\|_{P_{0}}$
$\displaystyle=O_{p}(h^{r/2})O_{p}(n^{-1/3}).$

The rest follows from Lemma S3. Therefore,
$(P_{n}-P_{0})\frac{K_{h,\bm{v}_{0}}}{P_{0}K_{h,\bm{v}_{0}}}(\Psi_{n}-\Psi_{0})=o_{p}\left((nh^{r})^{-1/2}\right)$
when the bandwidth $h_{n}$ satisfies $h^{r}\rightarrow 0$ and
$nh^{r3/2}\rightarrow\infty$.

###  S4.4 Asymptotic behaviour of the forth term on the RHS of (S32)

The forth term can be written as

$\displaystyle
P_{0}\frac{K_{h,\bm{v}_{0}}}{P_{0}K_{h,\bm{v}_{0}}}(\Psi_{n}-\Psi_{0})=$
$\displaystyle
P_{0}\frac{K_{h,\bm{v}_{0}}\Delta^{c}I(A=d^{\theta})}{G_{0}P_{0}K_{h,\bm{v}_{0}}}(\Psi_{n}-\Psi_{0})$
$\displaystyle=$ $\displaystyle
P_{0}\frac{K_{h,\bm{v}_{0}}\Delta^{c}I(A=d^{\theta})}{G_{0}P_{0}K_{h,\bm{v}_{0}}}\\{\Psi_{n}-I(Y>t)\\}$
$\displaystyle=$
$\displaystyle(P_{n}-P_{0})\frac{K_{h,\bm{v}_{0}}\Delta^{c}I(A=d^{\theta})}{G_{0}P_{0}K_{h,\bm{v}_{0}}}\\{I(Y>t)-\Psi_{0}\\}$
$\displaystyle-
P_{n}\frac{K_{h,\bm{v}_{0}}\Delta^{c}I(A=d^{\theta})}{G_{0}P_{0}K_{h,\bm{v}_{0}}}\\{I(Y\geq
t)-\Psi_{n}\\}+o_{p}\left((nh^{r})^{-1/2}\right).$

The third equality follows from Section S4.3.

The last term on the RHS can be represented as

$\displaystyle
P_{n}\frac{K_{h,\bm{v}_{0}}\Delta^{c}I(A=d^{\theta})}{G_{0}P_{0}K_{h,\bm{v}_{0}}}\\{I(Y\geq
t)-\Psi_{n}\\}$
$\displaystyle=P_{n}\frac{K_{h,\bm{v}_{0}}\Delta^{c}I(A=d^{\theta})}{G_{n}P_{0}K_{h,\bm{v}_{0}}}\\{I(Y\geq
t)-\Psi_{n}\\}$
$\displaystyle+P_{n}\frac{K_{h,\bm{v}_{0}}\Delta^{c}I(A=d^{\theta})}{P_{0}K_{h,\bm{v}_{0}}}\\{I(Y\geq
t)-\Psi_{n}\\}\left(\frac{G_{n}-G_{0}}{G_{0}G_{n}}\right)$
$\displaystyle=P_{n}\frac{K_{h,\bm{v}_{0}}\Delta^{c}I(A=d^{\theta})}{G_{n}P_{0}K_{h,\bm{v}_{0}}}\\{I(Y\geq
t)-\Psi_{n}\\}+$
$\displaystyle(P_{n}-P_{0})\frac{K_{h,\bm{v}_{0}}\Delta^{c}I(A=d^{\theta})}{P_{0}K_{h,\bm{v}_{0}}}\\{I(Y\geq
t)-\Psi_{n}\\}\left(\frac{G_{n}-G_{0}}{G_{0}G_{n}}\right)$
$\displaystyle+P_{0}\frac{K_{h,\bm{v}_{0}}\Delta^{c}I(A=d^{\theta})}{P_{0}K_{h,\bm{v}_{0}}}\\{I(Y\geq
t)-\Psi_{n}\\}\left(\frac{G_{n}-G_{0}}{G_{0}G_{n}}\right).$

Using Lemma S3 and similar techniques used in Section S4.3, we have
$(P_{n}-P_{0})\frac{K_{h,\bm{v}_{0}}\Delta^{c}I(A=d^{\theta})}{P_{0}K_{h,\bm{v}_{0}}}\\{I(Y\geq
t)-\Psi_{n}\\}\left(\frac{G_{n}-G_{0}}{G_{0}G_{n}}\right)=o_{p}\left((nh^{r})^{-1/2}\right)$.

Let $Q_{0}=E\\{I(Y\geq t)|A=d^{\theta},\bm{W}\\}$. Then

$\displaystyle
P_{0}\frac{K_{h,\bm{v}_{0}}\Delta^{c}I(A=d^{\theta})}{P_{0}K_{h,\bm{v}_{0}}}\\{I(Y\geq
t)-\Psi_{n}\\}\left(\frac{G_{n}-G_{0}}{G_{0}G_{n}}\right)=$ $\displaystyle
P_{0}\frac{K_{h,\bm{v}_{0}}\Delta^{c}I(A=d^{\theta})}{P_{0}K_{h,\bm{v}_{0}}}(Q_{0}-\Psi_{n})\left(\frac{G_{n}-G_{0}}{G_{0}G_{n}}\right)$
$\displaystyle=$ $\displaystyle-
P_{0}\frac{K_{h,\bm{v}_{0}}G_{0}}{P_{0}K_{h,\bm{v}_{0}}}(Q_{0}-\Psi_{n})\left(\frac{G_{0}-G_{n}}{G_{0}^{2}}\right)$
$\displaystyle+P_{0}\frac{K_{h,\bm{v}_{0}}G_{0}}{P_{0}K_{h,\bm{v}_{0}}}(Q_{0}-\Psi_{n})\left(\frac{(G_{0}-G_{n})^{2}}{G_{n}G_{0}^{2}}\right).$

Using the rate of convergence of the highly adaptive lasso estimator
$P_{0}\frac{K_{h,\bm{v}_{0}}G_{0}}{P_{0}K_{h,\bm{v}_{0}}}(Q_{0}-\Psi_{n})\left(\frac{(G_{0}-G_{n})^{2}}{G_{n}G_{0}^{2}}\right)=o_{p}\left((nh^{r})^{-1/2}\right)$.
Moreover,

$\displaystyle
P_{0}\frac{K_{h,\bm{v}_{0}}(Q_{0}-\Psi_{n})}{P_{0}K_{h,\bm{v}_{0}}}\left(\frac{G_{0}-G_{n}}{G_{0}}\right)$
$\displaystyle=P_{0}\frac{K_{h,\bm{v}_{0}}(Q_{0}-\Psi_{n})}{P_{0}K_{h,\bm{v}_{0}}}\left(\frac{G_{0}^{c}G_{0}^{a}-G_{n}^{c}G_{n}^{a}}{G_{0}^{c}G_{0}^{a}}\right)$
$\displaystyle=P_{0}\frac{K_{h,\bm{v}_{0}}(Q_{0}-\Psi_{n})}{P_{0}K_{h,\bm{v}_{0}}}G_{0}^{a}\left(\frac{G_{0}^{c}-G_{n}^{c}}{G_{0}^{c}G_{0}^{a}}\right)+P_{0}\frac{K_{h,\bm{v}_{0}}(Q_{0}-\Psi_{n})}{P_{0}K_{h,\bm{v}_{0}}}G_{n}^{c}\left(\frac{G_{0}^{a}-G_{n}^{a}}{G_{0}^{c}G_{0}^{a}}\right)$
$\displaystyle=-(P_{n}-P_{0})\frac{K_{h,\bm{v}_{0}}I(A=d^{\theta})}{G_{0}^{a}P_{0}K_{h,\bm{v}_{0}}}(Q_{0}-\Psi_{0})\left(\frac{\Delta^{c}-G_{0}^{c}}{G_{0}^{c}}\right)$
$\displaystyle\hskip
14.45377pt-(P_{n}-P_{0})\frac{K_{h,\bm{v}_{0}}(Q_{0}-\Psi_{0})}{P_{0}K_{h,\bm{v}_{0}}}\left(\frac{I(A=d^{\theta})-G_{0}^{a}}{G_{0}^{a}}\right)$
$\displaystyle\hskip
14.45377pt+P_{n}\frac{K_{h,\bm{v}_{0}}I(A=d^{\theta})}{G_{0}^{a}P_{0}K_{h,\bm{v}_{0}}}(Q_{0}-\Psi_{n})\left(\frac{\Delta^{c}-G_{n}^{c}}{G_{0}^{c}}\right)$
$\displaystyle\hskip
14.45377pt+P_{n}\frac{K_{h,\bm{v}_{0}}(Q_{0}-\Psi_{n})}{P_{0}K_{h,\bm{v}_{0}}}G_{n}^{c}\left(\frac{I(A=d^{\theta})-G_{n}^{a}}{G_{0}^{c}G_{0}^{a}}\right)+o_{p}\left((nh^{r})^{-1/2}\right).$

The first equality follows because

$\displaystyle(P_{n}-P_{0})\frac{K_{h,\bm{v}_{0}}I(A=d^{\theta})}{G_{0}^{a}P_{0}K_{h,\bm{v}_{0}}}(Q_{0}-\Psi_{n})\left(\frac{\Delta^{c}-G_{n}^{c}}{G_{0}^{c}}\right)$
$\displaystyle=(P_{n}-P_{0})\frac{K_{h,\bm{v}_{0}}I(A=d^{\theta})}{G_{0}^{a}P_{0}K_{h,\bm{v}_{0}}}(Q_{0}-\Psi_{0})\left(\frac{\Delta^{c}-G_{0}^{c}}{G_{0}^{c}}\right)$
$\displaystyle\hskip 180.67499pt+o_{p}\left((nh^{r})^{-1/2}\right),$
$\displaystyle(P_{n}-P_{0})\frac{K_{h,\bm{v}_{0}}(Q_{0}-\Psi_{n})}{P_{0}K_{h,\bm{v}_{0}}}G_{n}^{c}\left(\frac{I(A=d^{\theta})-G_{n}^{a}}{G_{0}^{a}G_{0}^{c}}\right)$
$\displaystyle=(P_{n}-P_{0})\frac{K_{h,\bm{v}_{0}}(Q_{0}-\Psi_{0})}{P_{0}K_{h,\bm{v}_{0}}}\left(\frac{I(A=d^{\theta})-G_{0}^{a}}{G_{0}^{a}}\right)$
$\displaystyle\hskip 180.67499pt+o_{p}\left((nh^{r})^{-1/2}\right).$

Moreover, by Lemma S4, we have

$\displaystyle
P_{n}\frac{K_{h,\bm{v}_{0}}I(A=d^{\theta})}{G_{0}^{a}P_{0}K_{h,\bm{v}_{0}}}(Q_{0}-\Psi_{n})\left(\frac{\Delta^{c}-G_{n}^{c}}{G_{0}^{c}}\right)$
$\displaystyle=P_{n}\frac{K_{h,\bm{v}_{0}}I(A=d^{\theta})}{G_{n}^{a}P_{0}K_{h,\bm{v}_{0}}}(Q_{n}-\Psi_{n})\left(\frac{\Delta^{c}-G_{n}^{c}}{G_{n}^{c}}\right)+o_{p}\left((nh^{r})^{-1/2}\right)$
$\displaystyle
P_{n}\frac{K_{h,\bm{v}_{0}}(Q_{0}-\Psi_{n})}{P_{0}K_{h,\bm{v}_{0}}}G_{n}^{c}\left(\frac{I(A=d^{\theta})-G_{n}^{a}}{G_{0}^{c}G_{0}^{a}}\right)$
$\displaystyle=P_{n}\frac{K_{h,\bm{v}_{0}}(Q_{n}-\Psi_{n})}{P_{n}K_{h,\bm{v}_{0}}}\left(\frac{I(A=d^{\theta})-G_{n}^{a}}{G_{n}^{a}}\right)+o_{p}\left((nh^{r})^{-1/2}\right)$

Hence, if we undersmooth $G_{n}^{a}$ and $G_{n}^{c}$ such that
$P_{n}\frac{K_{h,\bm{v}_{0}}I(A=d^{\theta})}{G_{n}^{a}P_{0}K_{h,\bm{v}_{0}}}(Q_{n}-\Psi_{n})\left(\frac{\Delta^{c}-G_{n}^{c}}{G_{n}^{c}}\right)=o_{p}\left((nh^{r})^{-1/2}\right)$
and
$P_{n}\frac{K_{h,\bm{v}_{0}}(Q_{n}-\Psi_{n})}{P_{n}K_{h,\bm{v}_{0}}}\left(\frac{I(A=d^{\theta})-G_{n}^{a}}{G_{n}^{a}}\right)=o_{p}\left((nh^{r})^{-1/2}\right)$,
then

$\displaystyle
P_{0}\frac{K_{h,\bm{v}_{0}}(Q_{0}-\Psi_{n})}{P_{0}K_{h,\bm{v}_{0}}}\left(\frac{G_{0}-G_{n}}{G_{0}}\right)$
$\displaystyle=-(P_{n}-P_{0})\frac{K_{h,\bm{v}_{0}}I(A=d^{\theta})}{G_{0}^{a}P_{0}K_{h,\bm{v}_{0}}}(Q_{0}-\Psi_{0})\left(\frac{\Delta^{c}-G_{0}^{c}}{G_{0}^{c}}\right)$
$\displaystyle\hskip
14.45377pt-(P_{n}-P_{0})\frac{K_{h,\bm{v}_{0}}(Q_{0}-\Psi_{0})}{P_{0}K_{h,\bm{v}_{0}}}\left(\frac{I(A=d^{\theta})-G_{0}^{a}}{G_{0}^{a}}\right)+o_{p}\left((nh^{r})^{-1/2}\right).$

###  S4.5 The influence function.

Gathering all the terms, we have,

$\displaystyle\Psi_{nh}(\bm{v}_{0},\theta)-\Psi_{0}(\bm{v}_{0},\theta)=$
$\displaystyle(P_{n}-P_{0})\frac{K_{h,\bm{v}_{0}}\Delta^{c}I(A=d^{\theta})}{G_{0}P_{0}K_{h,\bm{v}_{0}}}I(Y\geq
t)$
$\displaystyle-(P_{n}-P_{0})\frac{K_{h,\bm{v}_{0}}I(A=d^{\theta})}{G_{0}^{a}P_{0}K_{h,\bm{v}_{0}}}Q_{0}\left(\frac{\Delta^{c}-G_{0}^{c}}{G_{0}^{c}}\right)$
$\displaystyle-(P_{n}-P_{0})\frac{K_{h,\bm{v}_{0}}Q_{0}}{P_{0}K_{h,\bm{v}_{0}}}\left(\frac{I(A=d^{\theta})-G_{0}^{a}}{G_{0}^{a}}\right)$
$\displaystyle-(P_{n}-P_{0})\frac{K_{h,\bm{v}_{0}}}{P_{0}K_{h,\bm{v}_{0}}}\Psi_{0h}(\bm{v}_{0},\theta)+h^{J}B_{0}(J)+o_{p}\left((nh^{r})^{-1/2}\right).$

Thus, the asymptotic normality follows from the choice of bandwidth such that
$(nh^{r})^{1/2}h^{J}=O(1)$ which implies $h_{n}=n^{-1/(2J+r)}$. Thus,
$(nh^{r})^{1/2}\\{\Psi_{nh}(\bm{v}_{0},\theta)-\Psi_{0}(\bm{v}_{0},\theta)\\}$
converges to a mean zero normal random variable when
$\frac{h_{n}}{n^{-1/(2J+r)}}\rightarrow 0$ (i.e., $h_{n}$ converges to zero
faster than $n^{-1/(2J+r)}$).

###  S4.6 The optimal bandwidth rate.

The rate constraints obtained in Lemma S3 (i.e., $h>n^{-\frac{2}{3r}}$) and
Section S4.5 (i.e., $h<n^{-\frac{1}{r+2J}}$) imply that the optimal bandwidth
rate is $h_{n}=n^{-\frac{1}{r+2J}}$. Note that for the constraint
$n^{-\frac{2}{3r}}<h<n^{-\frac{1}{r+2J}}$ to hold we must have $J>r/4$ which
is a very mild condition. For example, for $r<4$, we can still pick $J=1$.

##  Appendix S5 Proof of Theorem 4.3: Asymptotic linearity of the regimen-
response curve estimator for a fixed bandwidth $h$

$\displaystyle\Psi_{nh}(\bm{v}_{0},\theta)-\Psi_{0h}(\bm{v}_{0},\theta)=$
$\displaystyle
P_{n}K_{h,\bm{v}_{0}}\Psi_{n}\left(\frac{1}{P_{n}K_{h,\bm{v}_{0}}}-\frac{1}{P_{0}K_{h,\bm{v}_{0}}}\right)$
(S33)
$\displaystyle+(P_{n}-P_{0})\frac{K_{h,\bm{v}_{0}}}{P_{0}K_{h,\bm{v}_{0}}}(\Psi_{n}-\Psi_{0})+(P_{n}-P_{0})\frac{K_{h,\bm{v}_{0}}}{P_{0}K_{h,\bm{v}_{0}}}(\Psi_{0})$
$\displaystyle+P_{0}\frac{K_{h,\bm{v}_{0}}}{P_{0}K_{h,\bm{v}_{0}}}(\Psi_{n}-\Psi_{0}).$

###  S5.1 Asymptotic behaviour of the first term on the RHS of (S33)

Following the same steps used in the proof of Theorem 4.2 by consistency of
$\Psi_{n}$ and because, for fixed $h$,
$(P_{n}-P_{0})K_{h,\bm{v}_{0}}=O_{p}(n^{-1/2})$, we have

$\displaystyle