# Finding Optimal Dynamic Treatment Regimes Using Smooth Fisher Consistent
Surrogate Loss

Nilanjana<EMAIL_ADDRESS>[ Aaron Sonabend-
<EMAIL_ADDRESS>[ Rajarshi
<EMAIL_ADDRESS>[ Tianxi
<EMAIL_ADDRESS>[ Department of Biostatistics, Harvard
University, 677 Huntington Ave, Boston, MA 02115

###### Abstract

Large health care data repositories such as electronic health records (EHR)
opens new opportunities to derive individualized treatment strategies to
improve disease outcomes. We study the problem of estimating sequential
treatment rules tailored to patient’s individual characteristics, often
referred to as dynamic treatment regimes (DTRs). We seek to find the optimal
DTR which maximizes the discontinuous value function through direct
maximization of a fisher consistent surrogate loss function. We show that a
large class of concave surrogates fails to be Fisher consistent, which differs
from the classic setting for binary classification. We further characterize a
non-concave family of Fisher consistent smooth surrogate functions, which can
be optimized with gradient descent using off-the-shelf machine learning
algorithms. Compared to the existing direct search approach under the support
vector machine framework (Zhao et al., 2015), our proposed DTR estimation via
surrogate loss optimization (DTRESLO) method is more computationally scalable
to large sample size and allows for a broader functional class for the
predictor effects. We establish theoretical properties for our proposed DTR
estimator and obtain a sharp upper bound on the regret corresponding to our
DTRESLO method. Finite sample performance of our proposed estimator is
evaluated through extensive simulations and an application on deriving an
optimal DTR for treating sepsis using EHR data from patients admitted to
intensive care units.

, , and

00footnotetext: †: Equal Contributors

[]

##  1 Introduction

Due to the increasing adoption of electronic health records (EHR) and the
linkage of EHR with bio-repositories and other research registries, integrated
large datasets have become available for real world evidence based precision
medicine studies. These rich EHR data capture heterogeneity in response to
treatment over time and across patients, thereby offering unique opportunities
to optimize treatment strategies for individual patients over time. Sequential
treatment decisions tailored to patients’ individual characteristics at given
decision time points are often referred to as dynamic treatment regimes (DTRs)
in the statistical literature and reinforcement learning (RL) in the machine
learning literature. An optimal DTR can be defined as the sequential treatment
assignment rule that maximizes the expected counterfactual outcome, often
referred to as the value function in the DTR literature.

To estimate the optimal DTR, the most traditional approaches rely on modelling
the data-distribution or part of the data-distribution (Xu et al., 2016;
Zajonc, 2012). The most popular among the latter class are the regression-
based methods, including Q-learning, A-learning and marginal structural mean
models (Watkins, 1989; Murphy, 2003; Schulte et al., 2014; Orellana et al.,
2010; Robins, 2004). The regression-based methods, especially Q-learning,
offers the flexibility necessary for extension to a variety of settings
including, but not limited to, semi-supervised setting (Sonabend et al.,
2021), interactive model-building (Laber et al., 2014), discrete outcomes or
utilities (Moodie et al., 2014) etc. However, the underlying models in the
regression-based approaches are often high-dimensional, and susceptible to
mis-specification due to the sequential nature of the problem (Murphy et al.,
2001). Although A-learning and marginal structural mean models are more robust
to model mis-specification, they still require the contrast of Q-functions to
be correctly specified (cf. Schulte et al., 2014). These limitations of the
regression-based methods led the conception of the classification-based direct
search methods, which, in contrast, directly targets the counterfactual value
function.

The classification-based approaches essentially rely on the representation of
the counterfactual value function through importance sampling (Murphy et al.,
2001), whose maximization can be framed as a classification problem with
respect to the zero-one loss function (cf. Zhao et al., 2012, 2015; Chen et
al., 2016; Zhou and Kosorok, 2017; Song et al., 2015; Chen et al., 2017; Cui
and Tchetgen Tchetgen, 2020, and the references therein). The resulting
objective function is not amenable to efficient optimization owing to the
discontinuity of the zero-one loss. Therefore, following contemporary
classification literature (cf. Bartlett et al., 2006; Lin, 2004) the direct
search methods aim to replace the zero-one loss with alternative smoother
fisher consistent surrogate loss functions to facilitate efficient
classification methods. The paradigm shift of estimating DTRs by finding
classification rules is a powerful idea. Some authors indicate that existing
direct search methods outperform regression-based counterparts when the number
of stages is small (Kosorok and Laber, 2019; Luedtke and Van Der Laan, 2016).

Although initially developed for the one stage case, direct search method was
introduced to the multi-stage DTR by the novel work of Zhao et al. (2015).
Currently, it has two mainstream approaches. The first approach performs
binary classification stagewise in a backward fashion (cf. BOWL method of Zhao
et al., 2015; Jiang et al., 2019). However, at stage $t$, this approach can
only use those observations whose treatment assignment matches the optimal
treatment stage $t+1$ onward. As a result, the effective sample size of the
initial stages dwindles rapidly, which can be problematic during practical
implementation (Kosorok and Laber, 2019; Kallus, 2020). The other approach
builds on a simultaneous optimization method, which utilizes the whole data-
set for estimating each treatment assignment (SOWL, Zhao et al., 2015). While
it does not share the limitation of the BOWL-type approaches, this approach
hinges on a sequential weighted classification problem which is complicated by
the dependent nature of the DTR setting. Zhao et al. (2015) solves this
classification using a bivariate hinge-loss type surrogate. Although the idea
behind simultaneous optimization is powerful, the implementation via non-
smooth hinge-loss surrogate leads to a number of issues, scalability being one
of them. See Section 7 for more details. It is natural to ask whether the
hinge loss can be replaced by other surrogates. However, the answer is not
immediate because unlike BOWL, the simultaneous classification does not yield
to the binary classification theory on surrogate losses (Bartlett et al.,
2006). Although multicategory and multi-label classifications have apparent
resemblance with this classification problem, as we will see, they have
fundamental differences. This gives rise to the need for a unified study of
fisher consistent surrogate losses under the DTR setting. Our paper is the
first step towards that end.

For the ease of presentation, we focus on $k=2$ stage DTRs associated with two
time points in this paper. However, the main methodology easily extends to
general k-stage settings when $k>2$. Similar to most current works in direct
search methods, we consider only a binary treatment indicator, which is an
important practical case (Laber and Davidian, 2017). Direct search with multi-
level treatments would require substantially different techniques, and is out
of the scope of the present paper.

###  1.1 Main contributions:

In the sequel, we will refer to the classification problem resulting from the
simultaneous optimization approach as “the DTR classification problem” for
brevity. We will refer to our approach of achieving optimal DTR estimation via
surrogate loss optimization as DTRESLO.

##### Concave losses:

In Theorem 1, we establish that above-bounded smooth concave surrogates fail
to be Fisher consistent in the DTR context. The failure is not restricted to
only smooth concave surrogates since our Theorem 2 also shows that non-smooth
hinge loss also fails to be Fisher consistent. Moreover, we are unaware of any
concave Fisher consistent surrogate in the DTR context. Therefore, our
findings naturally raise the question of whether any concave loss function is
Fisher-consistent for this problem.

##### A class of Fisher consistent surrogates for DTR Estimation:

We introduce a class of Fisher consistent surrogate losses (see Theorem 3),
which are non-concave but smooth. Yielding to efficient gradient based
algorithms like stochastic gradient descent, these surrogates facilitate fast
and scalable optimization methods. Our approach offers flexibility for
learning the DTRs so that practitioners can tailor the method to the data and
problem at hand. In particular, the smoothness of our surrogate losses makes
the optimization problem suitable to a broad range of standard machine
learning algorithms including, but not limited to, neural networks, wavelet
series, and basis expansion. Interpretable treatment rules are also achievable
by coupling our DTRESLO method with interpretable classifiers, such as linear
or tree-based classifiers. Finally, since we optimize the primal objective
function, variable selection in our case is straight-forward via addition of
an $l_{1}$ penalty.

##### Theoretical guarantee for a class of DTR estimators

We provide sharp upper bound on the regret – the difference between the
optimal value function and the value attained by the estimated treatment
regime, with detailed analyses focused on searching for DTR within the neural
network classifiers. We perform a sharp analysis of our approximation error
(see Theorem 4) and estimation error under Tsybakov’s small noise condition
(Tsybakov et al., 2004). In Corollaries 1 and 3, we prove that provided the
optimization error is small, the regret of our DTRESLO method with neural
network and wavelet series classifiers decays at a fast rate. Here by fast, we
mean decay rate faster than $n^{-1/2}$ is achievable. It turns out that this
rate also matches the minimax rate of risk decay (upto a poly-logarithmic
factor) of binary classification under assumptions similar to ours (Audibert
et al., 2007). Since two stage DTR is unlikely to be simpler than one stage
DTR, we conjecture that that our rate is minimax-optimal (upto a poly-
logarithmic factor) in two stage DTR under our assumptions. In the special
case when treatment effects are bounded away from zero, we show that our
regret decays at the rate of $O_{p}(1/n)$ upto a poly-logarithmic order.

The rest of the article is organized as follows. In Section 2 we outline the
problem and discuss the mathematical formulation. In Section 3 we discuss
Fisher consistency in the DTR setting, show that a large class of concave
surrogates fail to be Fisher consistent, and establish the Fisher consistency
of a family of non-concave surrogates. In Section 4 we construct a method for
estimating the optimal DTRs using the Fisher consistent surrogates, and
discuss the potential sources of error that contribute to the regret. Section
5 and Section 6 are devoted towards obtaining theoretical upper bounds of the
regret of our DTRESLO method. Section 5 focuses on approximation error, which
is combined with the estimation error in Section 6 to yield the final regret
bound. Then in section 8 we illustrate our DTRESLO method’s empirical
performance with extensive simulations 111The code for our empirical study can
be found at the GitHub repository DTR-via-surrogate-loss. and an application
to a sepsis cohort. We continue with a discussion in Section 9. Additional
details and proofs of our theoretical results are deferred to the Supplement.

###  1.2 Notation

We let $\overline{\mathbb{R}}$ denote the extended real line
$\mathbb{R}\cup\\{\pm\infty\\}$ and write $\mathbb{R}_{+}$ for the positive
half line $\\{x\in\mathbb{R}:x>0\\}$. Denote by $\mathbb{N}=\\{1,2,\ldots\\}$
the set of all natural numbers and for any integer $t$, we let
$[t]=\\{1,2,...,t\\}$. We also let $\mathbb{Z}$ denote the set of all
integers. For $p\in\mathbb{N}$, we let $\|\cdot\|_{p}$ denote the $l_{p}$ norm
, i.e. for $v\in\mathbb{R}^{k}$,
$\|v\|_{p}=(\sum_{i=1}^{k}|v_{i}|^{p})^{1/p}$. If $v\in\mathbb{N}^{k}$, we
denote by $|v|_{1}$ the quantity $\sum_{i=1}^{k}v_{i}$.

For any probability measure $P$ and measurable function $f$, we denote by
$\|f\|_{P,k}$ the norm $\left(\int|f(x)|^{k}dP(x)\right)^{1/k}$. We will also
denote this norm kby $L_{k}(P)$. Also, $Pf$ will denote the integral $\int
fdP$. For a concave function $f:\mathbb{R}^{k}\mapsto\mathbb{R}$, the domain
$\text{dom}(f)$ will be defined as in (Hiriart-Urruty and Lemaréchal, 2004, p.
74), that is, $\text{dom}(f)=\\{x\in\mathbb{R}^{k}\ :\ f(x)>-\infty\\}$. For
$\psi:\mathbb{R}^{2}\mapsto\mathbb{R}$, we will denote

$\psi_{12}(x,y)=\partialderivative{\psi(x,y)}{x}{y}.$

For any $x\in\mathbb{R}$, we denote by $\sigma(x)$ the ReLU activation
function $x_{+}=\max(x,0)$. For any set $A$, use the notation $1[x\in A]$ to
denote the event $\\{x\in A\\}$. Also, we denote by $\text{int}(A)$ the
interior of the set $A$. The cardinality of $A$ will be denoted by $|A|$. In
this paper, we will use $C$ and $c$ to denote generic constants which may vary
from line to line.

Many results in this paper are asymptotic (in $n$) in nature and thus require
some standard asymptotic notations. If $a_{n}$ and $b_{n}$ are two sequences
of real numbers then $a_{n}\gg b_{n}$ (and $a_{n}\ll b_{n}$) implies that
${a_{n}}/{b_{n}}\rightarrow\infty$ (and ${a_{n}}/{b_{n}}\rightarrow 0$) as
$n\rightarrow\infty$, respectively. Similarly $a_{n}\gtrsim b_{n}$ (and
$a_{n}\lesssim b_{n}$) implies that
$\liminf_{n\rightarrow\infty}{{a_{n}}/{b_{n}}}=C$ for some $C\in(0,\infty]$
(and $\limsup_{n\rightarrow\infty}{{a_{n}}/{b_{n}}}=C$ for some
$C\in[0,\infty)$). Alternatively, $a_{n}=o(b_{n})$ will also imply $a_{n}\ll
b_{n}$ and $a_{n}=O(b_{n})$ will imply that $\limsup_{n\rightarrow\infty}\
a_{n}/b_{n}=C$ for some $C\in[0,\infty)$).

##  2 Mathematical formalism

We focus on the DTR estimation under a longitudinal setting where data are
collected over time periods indexed by $t\in\\{1,2\\}$. Let
$O_{t}\in\mathcal{O}_{t}\subset\mathbb{R}^{p_{t}}$ denote the $p_{t}$
dimensional vector of patient clinical variables collected at time $t$ and
$p=\max(p_{1},p_{2})$. At a given time $t$, a binary treatment decision
$A_{t}\in\\{\pm 1\\}$ is made for the patient and a response to such treatment
$Y_{t}\in\mathbb{R}$ is observed. Without loss of generality, we assume higher
values of response $Y_{t}$ are desirable. Let us denote the distribution
underlying the observed random vector
$\mathcal{D}=(O_{1},A_{1},Y_{1},O_{2},A_{2},Y_{2})$ by $\mathbb{P}$. Suppose
we sample $n$ i.i.d. observations from $\mathbb{P}$. The corresponding
empirical distribution function will be denoted by $\mathbb{P}_{n}$. Since
treatment decisions are often made based on all previous states including
prior treatments and responses, we define the patient history by

$H_{1}=O_{1},\mbox{ and }H_{2}=(O_{1},Y_{1},O_{2},A_{1}),$

where $H_{1}$ and $H_{2}$ take values in sets $\mathcal{H}_{1}$ and
$\mathcal{H}_{2}$, respectively. We denote by $\pi_{1}(a_{1}\mid H_{1})$ and
$\pi_{2}(a_{2}\mid H_{2})$ the propensity scores
$\mathbb{P}(A_{1}=a_{1}|H_{1})$ and $\mathbb{P}(A_{2}=a_{2}|H_{2})$,
respectively.

Our goal is to find the treatment regime
$d=(d_{1},d_{2}):\mathcal{H}_{1}\times\mathcal{H}_{2}\mapsto\\{\pm
1\\}\times\\{\pm 1\\}$ that maximizes the expected sum of rewards
$Y_{1}(d)+Y_{2}(d)$,

$V(d_{1},d_{2})=\mathbb{E}_{d}[Y_{1}(d)+Y_{2}(d)],$

where $Y_{t}(d)$ is the potential outcome associated with time $t\in\\{1,2\\}$
, and $\mathbb{E}_{d}$ is the expectation with respect to the data
distribution under regime $d$. To this end, first we make some assumptions on
the observed data distribution $\mathbb{P}$ so that $V(d_{1},d_{2})$ becomes
identifiable under $\mathbb{P}$.

Assumptions for identifiability:

  * I.

Positivity: There exists a constant $C_{\pi}\in(0,1)$ so that
$\pi_{t}(A_{t}\mid H_{t})>C_{\pi}$ for all $H_{t}$, $t=1,2$.

  * II.

Consistency: The observed outcomes $Y_{t}$ and covariates $O_{t}$ agree with
the potential outcomes and covariates under the treatments actually received.
See Schulte et al. (2014); Robins (1994) for more details.

  * III

Sequential ignorability: For each $t=1,2$, the treatment assignment $A_{t}$ is
conditionally independent of the future potential outcomes $Y_{t}$ and future
potential clinical profile $O_{t+1}$ given $H_{t}$. Here we take $O_{3}$ to be
the empty set.

Our version of sequential ignorability follows from Robins (1997); Murphy et
al. (2001). Assumptions I-III are standard in DTR literature (Schulte et al.,
2014; Sonabend et al., 2021; Murphy et al., 2001; Zhao et al., 2015, e.g.).

Under Assumptions I-III, $\mathbb{E}_{d}(Y_{1}+Y_{2})$ can be identified under
$\mathbb{P}$ as (Zhao et al., 2015)

$\displaystyle\mathbb{E}_{d}[Y_{1}(d)+Y_{2}(d)]=\mathbb{P}\bigg{[}\dfrac{(Y_{1}+Y_{2})1[A_{1}=d_{1}(H_{1})]1[A_{2}=d_{2}(H_{2})]}{\pi_{1}(A_{1}\mid
H_{1})\pi_{2}(A_{2}\mid H_{2})}\bigg{]}.$

Furthermore, the optimal DTR,
$d^{*}=(d_{1}^{*},d_{2}^{*})=\operatorname*{arg\,max}_{d_{1}\in\mathcal{H}_{1},d_{2}\in
H_{2}}\mathbb{E}_{d}[Y_{1}(d)+Y_{2}(d)]$ are identified as follows:

$\displaystyle d_{2}^{*}(H_{2})$
$\displaystyle=\operatorname*{arg\,max}_{a_{2}\in\\{\pm
1\\}}\mathbb{E}[Y_{2}\mid H_{2},A_{2}=a_{2}]$ $\displaystyle d_{1}^{*}(H_{1})$
$\displaystyle=\operatorname*{arg\,max}_{a_{1}\in\\{\pm
1\\}}\mathbb{E}\left[Y_{1}+U_{2}^{*}(H_{2})\mid H_{1},A_{1}=a_{1}\right],$ (1)

where

$U_{2}^{*}(H_{2})=\max_{a_{2}\in\\{\pm
1\\}}\mathbb{E}[Y_{2}|H_{2},A_{2}=a_{2}].$ (2)

Since the optimal decision rules remain unchanged if a constant $C$ is added
to both $Y_{1}$ and $Y_{2}$, in what follows, unless otherwise mentioned, we
assume that $Y_{1},Y_{2}>C$ for some $C>0$. This trick was also used in Zhao
et al. (2015).

######  Remark 1 (Uniqueness of $d_{1}^{*}$ and $d_{2}^{*}$).

It is worth noting that $d_{1}^{*}$ and $d_{2}^{*}$ defined in (2) may not be
unique because they are allowed to take any value in $\\{\pm 1\\}$ at the
boundary. To elaborate on this further, suppose some $H_{2}$ satisfies

$\mathbb{E}[Y_{2}|H_{2},A_{2}=1]=\mathbb{E}[Y_{2}|H_{2},A_{2}=-1].$

Such values of $H_{2}$ constitute the decision boundary for the second stage.
Then both versions $d_{2}(H_{2})=1$ and $d^{\prime}_{2}(H_{2})=-1$ qualify as
optimal rule for at $H_{2}$. Similarly, for $d_{1}^{*}$, we can show that if
$H_{1}$ belongs to the first stage decision boundary

$\bigg{\\{}h_{1}\in\mathcal{H}_{1}\ :\
\mathbb{E}\Big{[}Y_{1}+U_{2}^{*}(H_{2})\bigg{|}A_{1}=1,h_{1}\Big{]}=\mathbb{E}\Big{[}Y_{1}+U_{2}^{*}(H_{2})\bigg{|}A_{1}=-1,h_{1}\Big{]}\bigg{\\}},$

then $d_{1}^{*}(H_{1})$ can take either value $+1$ or $-1$. Thus, $d_{1}^{*}$
is not unique either. Consequently, to avoid confusion, we let $d_{1}^{*}=1$
and $d_{2}^{*}=1$ at both first and second stage decision boundaries. In what
follows, we shall also refer to this optimal rule as “the optimal rule”. ∎

There is an alternative way of formulating $d^{*}$. If $(f_{1}^{*},f_{2}^{*})$
is a maximizer of

$V(f_{1},f_{2})=\mathbb{P}\bigg{[}\dfrac{(Y_{1}+Y_{2})1[A_{1}f_{1}(H_{1})>0]1[A_{2}f_{2}(H_{2})>0]}{\pi_{1}(A_{1}\mid
H_{1})\pi_{2}(A_{2}\mid H_{2})}\bigg{]}$ (3)

over the class

$\mathcal{F}=\bigg{\\{}(f_{1},f_{2})\ \bigg{|}\
f_{1}:\mathcal{H}_{1}\mapsto\mathbb{R},\quad
f_{2}:\mathcal{H}_{2}\mapsto\mathbb{R}\text{ are measurable}\bigg{\\}},$ (4)

then $\text{sign}(f_{1}^{*})$ and $\text{sign}(f_{2}^{*})$ yield the optimal
rules $d_{1}^{*}$ and $d_{2}^{*}$, respectively (Zhao et al., 2015). If
$f_{1}^{*}$ and $f_{2}^{*}$ take the value zero, then $d_{1}^{*}$ and
$d_{2}^{*}$ can be either $+1$ or $-1$. Finally, even if $d_{1}^{*}$ and
$d_{2}^{*}$ are unique, $f_{1}^{*}$ and $f_{2}^{*}$ need not be unique.

At this stage, although it is intuitive to consider maximization of the sample
analogue of $V(f_{1},f_{2})$ to estimate the optimal decision rule, the non-
concavity and discontinuity of the zero-one loss function render the
maximization of $V_{\psi}(f_{1},f_{2})$ computationally hard. To deal with
issues of similar flavor, the classification literature (cf. Bartlett et al.,
2006) suggests using a suitable surrogate to the zero-one loss function. We
appeal to this very intuition and consider

$V_{\psi}(f_{1},f_{2})=\mathbb{P}\bigg{[}\dfrac{(Y_{1}+Y_{2})\psi\Big{(}A_{1}f_{1}(H_{1}),A_{2}f_{2}(H_{2})\Big{)}}{\pi_{1}(A_{1}\mid
H_{1})\pi_{2}(A_{2}\mid H_{2})}\bigg{]},$ (5)

where $\psi$ is some bivariate function. For example, Zhao et al. (2015) takes
$\psi(x,y)=\min(x-1,y-1,0)$, the bivariate concave version of the popular
hinge loss $\phi(x)=\max(1-x,0)$. Suppose

$(\tilde{f}_{1},\tilde{f}_{2})\in\operatorname*{arg\,max}_{(f_{1},f_{2})\in\mathcal{F}}V_{\psi}(f_{1},f_{2}),$
(6)

where $\mathcal{F}$ is as defined in (4). Of course, the set in the RHS of (6)
may be empty, or may not be a subset of $\mathcal{F}$. Even if the set of
maximizers is non-empty, $\tilde{f}_{1}$ and $\tilde{f}_{2}$ need not be
unique. Provided they exist, each $\tilde{f}_{1}$ and $\tilde{f}_{2}$
correspond to decision rules
$\tilde{d}_{1}(H_{1})=\text{sign}(\tilde{f}_{1}(H_{1}))$ and
$\tilde{d}_{2}(H_{2})=\text{sign}(\tilde{f}_{2}(H_{2}))$. If
$\tilde{f}(H_{t})=0$, then $\tilde{d}_{t}(H_{t})$ can be either $+1$ or $-1$.

Finally, we define excess risk in line with the excess risk in context of
classification. Letting $V^{*}=V(f_{1}^{*},f_{2}^{*})$ and
$V^{*}_{\psi}=V_{\psi}(\tilde{f}_{1},\tilde{f}_{2})$, we define the respective
regret and $\psi$-regret of using $(f_{1},f_{2})$ by

$V^{*}-V(f_{1},f_{2})\quad\text{and}\quad V^{*}_{\psi}-V_{\psi}(f_{1},f_{2}),$

respectively. Note that regret and the $\psi$-regret are always non-negative.

Throughout our paper, we will compare our DTR classification with binary
classification. Therefore, we will fix the notation for binary classification.
In the setting of binary classification, we have observations $X$ taking value
in an Euclidean space $\mathcal{X}$. Each $X$ is associated with a label $A$,
which plays the same role as our treatment assignments. The optimal rule or
the Bayes rule assigns label $+1$ if $\eta(X)=P(A=1\mid X)>1/2$ and label $-1$
otherwise (cf. Bartlett et al., 2006). If $\eta(X)=1/2$, both labels are
optimal. The Bayes rule minimizes the classification risk
$\mathcal{R}(f)=\mathbb{P}(Af(X)<0)$ over all measurable functions
$f:\mathcal{X}\to\mathbb{R}$. Also, we denote
$\mathcal{R}^{*}=\inf_{f}\mathcal{R}(f)$, where the infimum is taken over all
measurable functions. Replacing the zero-one loss with the surrogate $\phi$
results in the $\phi$-risk $\mathcal{R}_{\phi}(f)=\mathbb{E}[\psi(Af(X))]$. We
let $\mathcal{R}_{\phi}^{*}$ denote the optimized $\phi$-risk.

Some parallels with the DTR classification setting are immediate. For example,
$V(f_{1},f_{2})$, $V^{*}$, $V_{\psi}(f_{1},f_{2})$, and $V^{*}_{\psi}$
correspond to $R(f)$, $R^{*}$, $R_{\phi}(f)$, and $R^{*}_{\phi}$,
respectively. next, defining the maps
$\eta_{1}:\mathcal{H}_{1}\mapsto\mathbb{R}$ and
$\eta_{2}:\mathcal{H}_{2}\times\\{0,1\\}\mapsto\mathbb{R}$ by

$\eta_{1}(H_{1})=\frac{\mathbb{E}[Y_{1}+U_{2}^{*}(H_{2})|A_{1}=1,H_{1}]}{\mathbb{E}[Y_{1}+U_{2}^{*}(H_{2})|A_{1}=1,H_{1}]+\mathbb{E}[Y_{1}+U_{2}^{*}(H_{2})|A_{1}=-1,H_{1}]},$
(7)
$\eta_{2}(H_{2})=\frac{\mathbb{E}[Y_{1}+Y_{2}|A_{2}=1,H_{2}]}{\mathbb{E}[Y_{1}+Y_{2}|A_{2}=1,H_{2}]+\mathbb{E}[Y_{1}+Y_{2}|A_{2}=-1,H_{2}]},$
(8)

we observe that $\eta_{1}$ and $\eta_{2}$ play the same role in DTR setting as
the conditional probability $\eta$ in context of binary classification. To
elaborate, from the definitions of $d_{1}^{*}$ and $d_{2}^{*}$ in (2), it
follows that $d_{t}^{*}(H_{t})=+1$ if $\eta_{t}(H_{t})>1/2$, and $-1$
otherwise. Note also that the first stage and second stage decision boundaries
can be represented by the sets $\\{h_{1}:\eta_{1}(h_{1})=1/2\\}$ and
$\\{h_{2}:\eta_{2}(h_{2})=1/2\\}$.

##  3 Fisher consistency

A desirable $\psi$ should ensure that $\tilde{d}$ is consistent with $d^{*}$.
To concertize the idea, we need the concept of Fisher consistency.

######  Definition 1.

The surrogate $\psi$ is called Fisher consistent if for all $\mathbb{P}$
satisfying Assumption I-III, any $\\{f_{1n},f_{2n}\\}_{n\geq
1}\subset\mathcal{F}$ that satisfies

$V_{\psi}(f_{1n},f_{2n})\to V_{\psi}^{*},\quad\text{also satisfies}\quad
V(f_{1n},f_{2n})\to V^{*}.$

Our definition of Fisher consistency is in line with classification literature
(Bartlett et al., 2006). Note that Definition 1 does not require
$\tilde{f}_{1}$ and $\tilde{f}_{2}$ to exist or be measurable. However, if
$\tilde{f}_{1}$ and $\tilde{f}_{2}$ do exist, and they are in $\mathcal{F}$,
then Fisher consistency implies $V_{\psi}(\tilde{d})=V^{*}$, indicating
$\tilde{d}$ is a candidate for $d^{*}$. In context of binary classification,
the surrogate $\phi$ is Fisher consistent if and only if
$\mathcal{R}_{\phi}(f_{n})\to\mathcal{R}_{\phi}^{*}$ implies
$\mathcal{R}(f_{n})\to\mathcal{R}^{*}$, where $f_{n}$’s are measurable
functions mapping $\mathcal{X}$ to $\mathbb{R}$.

######  Remark 2 (Characterization of Fisher consistency).

In many classification problems, e.g. binary, multicategory, or multi-label
classification, Fisher consistency can be directly characterized by convex
hulls of points in the image space of $\psi$, and the related notion is known
as calibration (Bartlett et al., 2006; Zhang, 2010; Tewari and Bartlett, 2007;
Gao and Zhou, 2011). For example, Theorem 1 of Bartlett et al. (2006) shows
that a surrogate $\phi$ is Fisher-consistent for binary classification if and
only if the following condition holds.

######  Condition 1.

$\phi:\mathbb{R}\mapsto\mathbb{R}$ satisfies

$\sup_{x:x(2\eta-1)\leq
0}\Big{(}\eta\phi(x)+(1-\eta)\phi(-x)\Big{)}<\sup_{x\in\mathbb{R}}\Big{(}\eta\phi(x)+(1-\eta)\phi(-x)\Big{)}$

for all $\eta\in[0,1]$ such that $\eta\neq 1/2$.

However, due to the sequential nature of the DTR set-up, it is not easy to
represent Fisher consistency in terms of analytical properties of $\psi$. This
complicates the analysis of Fisher consistency in the DTR set-up. ∎

Traditionally, the first preference of surrogate losses have been the concave
(convex in context of minimization) surrogates because they ensure unique
optimum (Chen et al., 2017). In the binary setting, a univariate concave
surrogate $\phi$ is Fisher consistent if and only if it is differentiable at
$0$ with positive derivative (see Bartlett et al., 2006, Theorem 6). Many
commonly used univariate concave losses satisfy these conditions. We display
some of these in Figure 1. An important geometric property of these functions
is that they mimic the graph of the zero-one loss function. After proper
shifting and scaling, their image lies below that of the zero-one loss
function (see Figure 1). Of course, concavity is not necessary for
classification-calibration, and this geometric property is shared by non-
concave classification calibrated losses as well (see Lemma 9 of Bartlett et
al., 2006).

Figure 1: Plots of $\phi(x)$ vs $x$ for concave calibrated value function
$\phi$ for binary decision rules. Here are the functions, Zero-one:
$\phi(x)=1[x>0]$, Exponential: $\phi(x)=1-e^{-x}$, Hinge: $\phi(x)=\min(x,1)$,
Squared error: $\phi(x)=1-(1-x)^{2}$.

There are also classes of concave surrogates which are Fisher consistent for
multicategory classification with respect to the zero-one loss (Duchi et al.,
2018; Tewari and Bartlett, 2007; Neykov et al., 2016) or for multilabel
classification with respect to Hamming loss (Gao and Zhou, 2011). In that
light, it is not unnatural to expect concave surrogates will succeed in the
DTR classification setting as well. Unfortunately, as we will see in the next
section, this simple-minded extension of binary classification may not hold.

DTR classification bears resemblance with multilabel classification
(Dembczyński et al., 2012) but additional complication arises since $H_{2}$
contains $H_{1}$. Also, the Fisher consistency literature on multilabel
classification (Gao and Zhou, 2011) is based on Hamming loss and partial
ranking loss, which are substantially different from the zero-one loss. Our
problem also exhibits similarity with multiclass classification (Duchi et al.,
2018). However, a big difference arises because of the sequential structure.
Had $d_{1}$ been a map from $\mathcal{H}_{2}$ to $\\{\pm 1\\}$ similar to
$d_{2}$, existing theory on multiclass classification (Duchi et al., 2018)
could be readily used to provide conditions for a general function $\psi$ to
be Fisher consistent. However, during the treatment assignment $d_{1}$, one
has no knowledge of $A_{1}$, $Y_{1}$, and $O_{2}$. Tewari and Bartlett (2007)
and Zhang (2010) develop tools for general classification set-ups, but these
tools are too generalized for explosion of the specific sequential structure
of DTR classification. In fact, it is the binary classification, which seems
to have the most parallels with DTR classification.

(a) Bivariate zero-one loss: $1[x>0,y>0]$.

(b) Bivariate hinge loss: $\max(x-1,y-1,0)$.

(c) Bivariate squared error loss: $-(x^{2}+y^{2})$.

(d) Bivariate exponential loss: $-\exp(-(x+y))$.

(e) Bivariate logistic loss: $-\log(1+\exp(-x)+\exp(-y))$.

Figure 2: Bivariate zero-one loss function and some concave surrogate losses

###  3.1 Concave surrogates

In this section, we will establish that a large class of concave surrogates
fail to be Fisher consistent for DTR estimation. We first consider the case of
smooth concave losses because they amend to gradient based optimization
methods with good scalability properties. We will start our discussion with an
example. The smooth concave function $\phi(x)=-\exp(-x)$ is Fisher consistent
in the binary classification setting. Let us consider its bivariate extension
$\psi(x,y)=-\exp(-x-y)$. It tuns out that $\tilde{d}_{1}(H_{1})$ takes the
form

$\operatorname*{arg\,max}_{a_{1}\in\\{\pm
1\\}}\mathbb{E}\Big{[}h\Big{(}Y_{1}+\mathbb{E}[Y_{2}\mid
H_{2},A_{2}=1],Y_{1}+\mathbb{E}[Y_{2}\mid H_{2},A_{2}=-1\Big{)}\ \Big{|}\
H_{1},A_{1}=a_{1}\Big{]},$

where $h(x,y)=\sqrt{xy}$. However, $d_{1}^{*}(H_{1})$ takes the same form but
with $h(x,y)=\max(x,y)$. In general, therefore, $\tilde{d}_{1}(H_{1})$ and
$d_{1}^{*}(H_{1})$ do not agree. To see this, consider the toy example when
$Y_{1}=1$ and

$\displaystyle Y_{2}=$ $\displaystyle\ 4\times 1[A_{1},A_{2}=1]+3\times
1[A_{1}=1,A_{2}=-1]$ $\displaystyle\ +5\times
1[A_{1}=-1,A_{2}=1]+1[A_{1},A_{2}=-1].$

In this case, $d_{1}^{*}(H_{1})=-1$ but $\tilde{d}_{1}(H_{1})=1$ for all
$H_{1}$, and clearly, $\psi$ is not Fisher consistent. If we consider other
examples of smooth concave $\psi$, e.g. logistic or quadratic loss, we obtain
different $h$, but for these examples as well, $h$ is quite different from the
non-smooth $h(x,y)=\max(x,y)$.

The above heuristics indicate that the criteria of DTR Fisher consistency may
be incompatible with smooth concave losses. Theorem 1 below concretize the
above heuristics for an important class of concave smooth losses. Theorem 1
assumes that $\psi$ is closed and strictly concave. We say a function is
closed if it is upper semicontinuous everywhere, or equivalently, if its
superlevel sets are closed (pp. 78, Hiriart-Urruty and Lemaréchal, 2004). The
function $h$ is strictly concave if for any $\lambda\in(0,1)$, and
$x,y\in\text{dom}(h)$,

$h(\lambda x+(1-\lambda)y)>\lambda h(x)+(1-\lambda)h(y).$

######  Theorem 1.

Suppose $\psi$ is closed, strictly concave, and bounded above. In addition,
$\psi$ has continuous second order partial derivatives and $\psi_{12}$ has
continuous partial derivatives on $\text{int}(\text{dom}(\psi))$. Then $\psi$
can not be Fisher consistent for two stage DTR.

We list below some examples of $\psi$, also shown in Figure 2, which satisfy
the assumptions of Theorem 1.

Exponential: $\displaystyle\psi(x,y)=-\exp(-x-y);$ Logistic:
$\displaystyle\psi(x,y)=-\log(1+e^{-x}+e^{-y});$ Quadratic:
$\displaystyle\psi(x,y)=z^{T}Qz+b^{T}z+c,$ where $z=(x,y)^{T}$, $Q$ is
negative definite, and $b\in\mathbb{R}^{2}$.

The proof of Theorem 1 is given in Supplement 16.

######  Remark 3 (Main challenges in the proof of Theorem 1).

A main difficulty in proving Theorem 1 is that even under our pathological
case, $\tilde{f}_{1}(H_{1})$ and $\tilde{f}_{2}(H_{2})$ do not have closed
form expressions. They are implicitly defined as maximizers of complex
functionals of $\psi$. Therefore, if we consider a very large class of
$\psi$’s, characterization of $\tilde{f}_{1}(H_{1})$ and
$\tilde{f}_{2}(H_{2})$ becomes difficult. The assumptions on $\psi$ ensure
that the class of $\psi$’s under consideration is manageable, mitigating some
technical difficulties in the characterization of $\tilde{f}_{1}(H_{1})$ and
$\tilde{f}_{2}(H_{2})$. The latter is essential for learning the behaviour of
the signs of $\tilde{f}_{1}(H_{1})$ and $\tilde{f}_{2}(H_{2})$. Thus the
assumptions on $\psi$ are required for technical reasons in the proof. That is
to say that our conditions on $\psi$ are probably not necessary, and the
assertions of Theorem 1 may hold even without these assumptions. In fact, we
are not aware of any concave surrogates that are Fisher consistent in this
context. We defer further discussion on the assumptions in Theorem 1 to
Supplement 12.1. ∎

The failure of concave surrogates is not simply an artifact of the smoothness
of $\psi$ in Theorem 1. Theorem 2, which is proved in Supplement 17, shows
that a commonly used non-smooth loss $\psi$, a concave version of the
bivariate hinge loss $\max(x-1,y-1,0)$, is not Fisher consistent. See Figure
2(b) for a pictorial representation of this loss. In fact, we are not aware of
any concave surrogate, smooth or not, which is Fisher consistent for the DTR
classification problem.

######  Theorem 2 (Hinge loss not Fisher consistent).

The surrogate $\psi(x,y)=\min(x-1,y-1,0)$ is not Fisher consistent with
respect to two stage DTR.

######  Remark 4 (Pathological cases).

Our counter examples for Theorem 1 and Theorem 2 are based on a pathological
case where $O_{2}$ and $Y_{1}$ are deterministic functions of $H_{1}$. We
chose this case because it grants technical simplification. The realistic
cases are no more likely to yield under concave surrogates than this simple
pathological case. The calculations underlying the proof of Theorem 1 become
severely technically challenging when the second stage covariates are
potentially random given $H_{1}$. ∎

The Hinge loss has previously been shown as not Fisher consistent in some
classification problems (Liu, 2007). Smooth concave or convex surrogates can
also fail to be fisher consistent for complex machine learning problems. For
example, Gao and Zhou (2011) shows known convex surrogates are not Fisher
consistent for multilabel classification with ranking loss. Ranking is another
notable example, where convex losses fail for a number of losses including the
pairwise disjoint loss (Calauzenes et al., 2012). In fact, in the latter case,
the existence of a Fisher consistent concave surrogate would imply the
feedback arc-set problem to be polynomial-time solvable (Duchi et al., 2010),
which is conjectured to be NP complete (Karp, 1972).

The DTR classification problem shares one common feature with the above-stated
machine-learning problems where these surrogate losses fail. It does not
_organically_ reduce to a sequence of weighted binary classification problem,
which is apparently a common ingredient of all the classification problems
which are polynomial time solvable via convex surrogates, e.g. multicategory
loss with zero one loss function (Tewari and Bartlett, 2007), multilabel
classification with partial ranking and hamming loss (Gao and Zhou, 2011),
ranking with Hamming loss (Calauzenes et al., 2012), ordinal regression with
absolute error loss (Pedregosa et al., 2017) etc. Here we stress on the word
“organic” because DTR classification does reduce to sequences of binary
classification if it is framed as a sequential classification via exclusion of
data points at each stage; cf. BOWL (Zhao et al., 2015).

######  Remark 5.

Zhao et al. (2015)’s SOWL method bases on the bivariate hinge loss of Theorem
2, which was claimed to yield $\tilde{d}=d^{*}$. While there may be
distributions under which the latter holds, our analysis indicates that the
bivariate hinge loss is not Fisher consistent in the sense of Definition
1222Our counterexample shows that for this $\psi$, the criteria for Fisher
consistency fails for some $\mathbb{P}$’s for which $O_{2}$ and $Y_{1}$ are
deterministic functions of $H_{1}$. See Remark 4 for more details. . ∎

While we do not have an intuitive explanation for the apparent failure of such
a wide class of concave functions, we attempt at making one heuristic
reasoning. Even in the one stage case of binary classification, it was
observed that Fisher consistency requires the surrogates to mimic the shape of
the zero-one loss to some extent. It appears to us that for Fisher consistency
in two stage DTR, the function $\psi$ has to mimic the shape of the bivariate
zero-one loss function (see Figure 2(a)) more closely than that was necessary
in binary classification (see Figure 1). In other words, the non-concavity of
the zero-one loss function at the origin pushes the concave losses to failure,
thereby necessitating search for $\psi$ among non-concave losses.

###  3.2 Construction of Fisher consistent surrogates

In this section, we construct Fisher consistent loss functions for two stage
DTR classification. Noting the connection between binary classification and
DTR classification, we consider bivariate loss functions of form
$\psi(x,y)=\phi_{1}(x)\phi_{2}(y)$ where $\phi_{1}$ and $\phi_{2}$ themselves
are univariate loss functions. The most intuitive choice of $\phi_{i}$’s would
be the Fisher consistent losses for one stage DTR. However,
$\phi_{i}(x)=-\exp(-x)$ is Fisher consistent in one stage (Bartlett et al.,
2006; Chen et al., 2017) although the product $\phi_{1}(x)\phi_{2}(y)$ is
inconsistent for the two stage setting (see Section 3). The above indicates
$\phi_{i}$’s simply being Fisher consistent would not give the necessary
structure to $\psi$ so as to mimic the bivariate zero-one loss function
effectively.

In fact, our calculations hint that $\phi_{2}$ needs to share a particular
property of the zero-one loss function, that is for some constant $C>0$,

$\displaystyle\sup_{x\in\mathbb{R}}\Big{(}\eta\phi_{2}(x)+(1-\eta)\phi_{2}(-x)\Big{)}=C\max(\eta,1-\eta).$
(9)

The above property is satisfied by the sigmoid function, which is non-concave,
and Fisher consistent for binary classification (Bartlett et al., 2006).
Interestingly, (9) alone does not guarantee the fisher consistency of
$\psi=\phi_{1}\phi_{2}$. However, we can circumvent the challenges faced by
the concave surrogates if $\phi_{1}$ and $\phi_{2}$ are either the sigmoid or
have shape similar to the sigmoid function defined by the following condition

######  Condition 2.

$\phi$ is a strictly increasing function such that

  1. 1.

$\phi(x)>0$ for all $x\in\mathbb{R}$.

  2. 2.

For all $x\in\mathbb{R}$, $\phi(x)$ satisfies $\phi(x)+\phi(-x)=C_{\phi}$
where $C_{\phi}>0$ is a constant.

  3. 3.

$\lim_{x\to\infty}\phi(x)=C_{\phi}$ and $\lim_{x\to-\infty}\phi(x)=0$.

We show in Lemma 16 in Supplement 18.3.1 that a $\phi$ satisfying Condition 2
is Fisher consistent and satisfies (9). We provide some examples of functions
satisfying Condition 2 below.

######  Example 1.

The following odd functions are non-decreasing with range $[-1,1]$:

  1. 1.

$f_{a}(x)=\frac{x}{1+|x|}$.

  2. 2.

$f_{b}(x)=\frac{2}{\pi}\arctan\bigg{(}\frac{\pi x}{2}\bigg{)}$.

  3. 3.

$f_{c}(x)=\frac{x}{\sqrt{1+x^{2}}}$.

  4. 4.

$f_{d}(x)=\tanh(x)$, where $\tanh(x)=\frac{e^{x}-e^{-x}}{e^{x}+e^{-x}}$.

Then $\phi_{\jmath}(x)=1+f_{\jmath}(x)$ satisfies Condition 2 with
$C_{\phi}=2$, for $\jmath=a,b,c,d$. See Figure 3(a) for the pictorial
representation of these functions.

(a) Plot of $\phi$’s.

(b) Plot of the surrogate $\psi$’s obtained from the above $\phi$’s.

Figure 3: Plot of the $\phi$’s in Example 1 and the corresponding $\psi$’s.
Here $\psi(x,y)=\phi(x)\phi(y)$.

A property like (9) is not exclusive to sigmoid-type functions or those
satisfying Condition 2. Functions like $\phi(x)=\min\\{x-1,0\\}$ satisfy this
property. However, we do not pursue analysis for this $\phi$ because is
neither smooth nor does it lead to a concave $\psi$. Condition 2 is a
sufficient condition for fisher consistency but not a necessary condition.
However, our theoretical results suggest that the class of surrogates
satisfying Condition 2 produces a tight bound on the regret, which seems to
get reflected in its numerical performance. Finally, note that the shape of
the functions satisfying Condition 2 mimic the shape of the bivariate zero-one
loss function more closely than do the concave surrogates in Figure 2. To see
this, compare Figure 3(b) with Figure 2.

####  3.2.1 Fisher consistency of $\psi$

Instead of directly proving Fisher consistency, we will bound the true regret
$V^{*}-V(f_{1},f_{2})$ in terms of the $\psi$-regret
$V_{\psi}^{*}-V_{\psi}(f_{1},f_{2})$. The benefit of such bound is that the
rate of convergence of the true regret will be readily given by that of the
$\psi$-regret, which we actually minimize.

As mentioned earlier, the true regret and the $\psi$-regret parallel the
excess risk $\mathcal{R}(f)-\mathcal{R}^{*}$ and the $\phi$-excess risk
$\mathcal{R}_{\phi}(f)-\mathcal{R}_{\phi}^{*}$ in binary classification. The
relationship between the latters have been well-studied. For Fisher consistent
$\phi$, Bartlett et al. (2006) show that

$h_{\phi}\Big{(}\mathcal{R}(f)-\mathcal{R}^{*}\Big{)}\leq\mathcal{R}_{\phi}(f)-\mathcal{R}_{\phi}^{*},$

where $h_{\phi}$ is a convex function satisfying $h_{\phi}(0)=0$. In view of
the fact that the univariate sigmoid loss leads to a linear $h_{\phi}$
(Bartlett et al., 2006, Example 4), it is reasonable to expect that a similar
inequality holds when $\psi(x,y)=\phi(x)\phi(y)$ with $\phi$ as in Condition
2, as confirmed in following theorem.

######  Theorem 3.

Suppose $Y_{1},Y_{2}>0$ and Assumptions I-III hold. Let
$\psi(x,y)=\phi(x)\phi(y)$ with $\phi$ satisfying the conditions of Lemma 16
with some $C_{\phi}>0$. Then

$V^{*}-V(f_{1},f_{2})\leq\frac{\Big{(}V_{\psi}^{*}-V_{\psi}(f_{1},f_{2})\Big{)}}{(C_{\phi}/2)^{2}}.$
(10)

Theorem 3 immediately implies Fisher consistency because if
$V_{\psi}(f_{1n},f_{2n})$ converges to $V_{\psi}^{*}$ for some $(f_{1n}$,
$f_{2n})\in\mathcal{F}$, then $V(f_{1n},f_{2n})\to V^{*}$ as well. Theorem 3
is proved in Supplement 18.2.

######  Remark 6 (Scaling of the $\phi$’s).

The scaling factor $C_{\phi}/2$ appears in the regret of (10) because $\phi$
differ from the zero-one function in scale by a factor of $C_{\phi}/2$. To
understand the impact of the scaling factor in the regret bound, suppose
$\phi_{2}=a\phi_{2}$ for some $a>0$, and
$\psi_{t}(x,y)=\phi_{t}(x)\phi_{t}(y)$ for $t=1,2$. Then

$\frac{\Big{(}V_{\psi_{1}}^{*}-V_{\psi_{1}}(f_{1},f_{2})\Big{)}}{C_{\phi_{1}}^{2}/4}=\frac{\Big{(}V_{\psi_{2}}^{*}-V_{\psi_{2}}(f_{1},f_{2})\Big{)}}{C_{\phi_{2}}^{2}/4}.$

Thus, the regret bound in (10) does not depend on the scale of $\phi$.
Nevertheless, during our implementation, we take the scaling factor
$C_{\phi}/2$ to be one so that the surrogate loss is at the same scale with
the original zero-one loss.

##  4 Main methodology

In this section, we describe how we use the Fisher-consistent surrogate
derived in Section 3.2 to estimate the optimal treatment regimes. Suppose
$\psi(x,y)=\phi(x)\phi(y)$ where $\phi$ satisfies Condition 2. Define the
empirical $\psi$-value function

$\widehat{V}_{\psi}(f_{1},f_{2})=\mathbb{P}_{n}\bigg{[}\dfrac{(Y_{1}+Y_{2})\psi\Big{(}A_{1}f_{1}(H_{1}),A_{2}f_{2}(H_{2})\Big{)}}{\pi_{1}(A_{1}\mid
H_{1})\pi_{2}(A_{2}\mid H_{2})}\bigg{]}$ (11)

Because $\mathbb{P}$ is unknown, we maximize $\widehat{V}_{\psi}(f_{1},f_{2})$
instead of $V_{\psi}(f_{1},f_{2})$. Ideally, one should maximize
$\widehat{V}_{\psi}(f_{1},f_{2})$ over $\mathcal{F}$ but brute force search
over $\mathcal{F}$ is impossible unless $\mathcal{H}_{1}$ and
$\mathcal{H}_{2}$ are discrete spaces with finite cardinality. Therefore, in
practice, one may optimize $\widehat{V}_{\psi}(f_{1},f_{2})$ over a nested
class

$\mathcal{U}_{1}\subset\ldots\subset\mathcal{\mathcal{U}}_{n}\subset\mathcal{F},$

where $\mathcal{U}_{n}$ is some rich class of classifiers, preferably a
universal class (see Zhang et al., 2018a). We will discuss them in more detail
later in Section 5.

Whatever is the choice of $\mathcal{U}_{n}$, maximization of
$\widehat{V}_{\psi}(f_{1},f_{2})$ over $(f_{1},f_{2})\in\mathcal{U}_{n}$
generally leads to a non-convex optimization problem. Our classifiers may be
local maximizers of $\widehat{V}_{\psi}(f_{1},f_{2})$ over
$(f_{1},f_{2})\in\mathcal{U}_{n}$, but it may also be a saddle point, or it
may not even be a critical point at all if early stopping is used. Therefore,
our classifiers may be far from the global optimizers of
$\widehat{V}_{\psi}(f_{1},f_{2})$ over $(f_{1},f_{2})\in\mathcal{U}_{n}$. That
is not a huge concern for us, however, as long as the empirical surrogate
value function at our classifiers and the said global optimizers are close.
That is because if the surrogate value functions of two function-pairs are
close, then their true value functions are also close by Theorem 3. The regret
bounds in Section 6 concertizes this idea. From a clinical point of view,
finding a treatment regime with high value function is of more significance
than finding a treatment regime close to the optimal regime.

The surrogate loss based DTR optimization allows flexibility in the choice of
$\mathcal{U}_{n}$ and the modification of the empirical loss
$\widehat{V}_{\psi}(f_{1},f_{2})$ to accommodate high dimensional covariates,
non-linear effects, and variable selection. One can maximize
$\widehat{V}_{\psi}(f_{1},f_{2})+\mathcal{P}(f_{1},f_{2})$ instead of
$\widehat{V}_{\psi}(f_{1},f_{2})$ to enable variable selection and attain
stable estimation, where $\mathcal{P}(f_{1},f_{2})$ is a penalty term. One can
include complex basis functions in $\mathcal{U}_{n}$ to incorporate non-linear
effects. For example, tree and list based methods (Zhang et al., 2018b; Sun
and Wang, 2021; Laber and Zhao, 2015) as well as neural networks [REF] can be
potentially adapted to construct $\mathcal{U}_{n}$. Our method can be extended
to $K$ stages by taking $\psi(x_{1},\ldots,x_{k})=\prod_{i=1}^{k}\phi(x_{i})$.

###  4.1 Decomposition of errors

In this section, we will discuss the decomposition of $\psi$-regret of DTRESLO
into three sources of errors. To that end, let us denote our classifiers by
$(\widehat{f}_{n,1},\widehat{f}_{n,2})$. We will provide upper bounds for the
$\psi$-regret $V_{\psi}(\widehat{f}_{n,1},\widehat{f}_{n,2})$, which readily
produces an upper bound for the true regret
$V(\widehat{f}_{n,1},\widehat{f}_{n,2})$ by Theorem 3. Before going into
further detail, we point out that

$V(\widehat{f}_{n,1},\widehat{f}_{n,2})=\mathbb{P}\bigg{[}(Y_{1}+Y_{2})\frac{1[A_{1}\phi(\widehat{f}_{n,1}(H_{1}))>0]1[A_{2}\phi(\widehat{f}_{n,2}(H_{2}))>0]}{\pi_{1}(A_{1}\mid
H_{1})\pi_{2}(A_{2}\mid H_{2})}\bigg{]}$

is a random quantity because here we assume that $H_{t},A_{t},Y_{t}$ $(t=1,2)$
are drawn from $\mathbb{P}$ independent of
$(\widehat{f}_{n,1},\widehat{f}_{n,2})$. The same holds regarding the
$\psi$-regret $V_{\psi}(\widehat{f}_{n,1},\widehat{f}_{n,2})$. We decompose
the $\psi$-regret according to three sources of errors: (i) approximation
error due to the approximation of $\mathcal{F}$ by $\mathcal{U}_{n}$; (ii)
estimation error due to the use of finite sample; and (iii) optimization error
due to the possibility of not achieving global maximization for
$\widehat{V}_{\psi}(f_{1},f_{2})$ since $\psi$ is non-concave. We define the
optimization error as

$\texttt{Opt}_{n}=\sup_{(f_{1},f_{2})\in\mathcal{U}_{n}}\widehat{V}_{\psi}(f_{1},f_{2})-\widehat{V}_{\psi}(\widehat{f}_{n,1},\widehat{f}_{n,2}).$

We first provide some heuristics and intuitions for the error decomposition.
For the time being, let us assume that
$\operatorname*{arg\,max}_{(f_{1},f_{2})\in\mathcal{U}_{n}}V_{\psi}^{*}(f_{1},f_{2})$
is attained at some $(\tilde{f}_{n,1},\tilde{f}_{n,2})\in\mathcal{U}_{n}$. The
existence of $(\tilde{f}_{n,1},\tilde{f}_{n,2})$ is not guaranteed in general,
and even if they exist, $(\tilde{f}_{n,1},\tilde{f}_{n,2})$ will be hard to
characterize for an arbitrary $\mathcal{U}_{n}$. We thus do not assume the
existence of $(\tilde{f}_{n,1},\tilde{f}_{n,2})$ in our proof. We define the
map
$\xi_{f_{1},f_{2},g_{1},g_{2}}:\mathcal{H}_{1}\times\mathcal{O}_{2}\times\mathbb{R}^{2}\times\\{\pm
1\\}^{2}\mapsto\mathbb{R}$ by

$\displaystyle\xi_{f_{1},f_{2},g_{1},g_{2}}(\mathcal{D})$ $\displaystyle:=$
$\displaystyle\
\frac{(Y_{1}+Y_{2})\left\\{\psi(A_{1}g_{1}(H_{1}),A_{2}g_{2}(H_{2}))-\psi(A_{1}f_{1}(H_{1}),A_{2}f_{2}(H_{2}))\right\\}}{\pi_{1}(A_{1}\mid
H_{1})\pi_{2}(A_{2}\mid H_{2})}.$ (12)

Elementary algebra shows that the $\psi$ regret can be decomposed as follows:

$\displaystyle V_{\psi}^{*}-V_{\psi}(\widehat{f}_{n,1},\widehat{f}_{n,2})$
$\displaystyle=$ $\displaystyle\
\underbrace{V_{\psi}^{*}-V_{\psi}(\tilde{f}_{n,1},\tilde{f}_{n,2})}_{\text{Approximation
error
}}+(V_{\psi}-\widehat{V}_{\psi})(\tilde{f}_{n,1},\tilde{f}_{n,2})-(V_{\psi}-\widehat{V}_{\psi})(\widehat{f}_{n,1},\widehat{f}_{n,2})$
$\displaystyle\
+\widehat{V}_{\psi}(\tilde{f}_{n,1},\tilde{f}_{n,2})-\sup_{(f_{1},f_{2})\in\mathcal{U}_{n}}\widehat{V}_{\psi}(f_{1},f_{2})+\underbrace{\sup_{(f_{1},f_{2})\in\mathcal{U}_{n}}\widehat{V}_{\psi}(f_{1},f_{2})-\widehat{V}_{\psi}(\widehat{f}_{n,1},\widehat{f}_{n,2})}_{\text{Optimization
error: }\texttt{Opt}_{n}}$ $\displaystyle\leq$ $\displaystyle\
\text{Approximation error
}+\underbrace{|(\mathbb{P}_{n}-\mathbb{P})[\xi_{f_{1},f_{2},\tilde{f}_{n,1},\tilde{f}_{n,2}}]|}_{\text{Estimation
error }}+\texttt{Opt}_{n}$ (13)

Clearly, $\texttt{Opt}_{n}$ depends on the optimization method used to
maximize $\widehat{V}_{\psi}(f_{1},f_{2})$ over $\mathcal{U}_{n}$. In this
paper, we do not analyze $\texttt{Opt}_{n}$, and focus on the approximation
error and the generalization error instead. However, our simulations in
Section 8 show that our DTRESLO methods are competitive to other popular
methods, indicating that $\texttt{Opt}_{n}$ does not weigh down the
performance of our DTRESLO methods relative to others. In our sharp analysis
of the $\psi$-regret, the estimation error bound depends on the approximation
error in an intricate manner; see Supplement 12.2 for more details. To keep
our presentation short and focused, therefore, we discuss the approximation
error in Section 5, and present the final regret-bound in Section 6. Details
on the explicit analysis of the estimation error can be found in Supplement
22. Finally, owing the potential non-concavity, the sharp analysis of the
$\psi$-regret is significantly more subtle compared to existing results and
approaches in the literature. We elaborate on this more in a detailed
discussion presented in Supplement 12.2.

In what follows, similar to Zhao et al. (2015), we assume that the propensity
scores, i.e. $\pi_{1}$ and $\pi_{2}$ are known. This will hold in particular
under a clinical trial like SMART (Kosorok and Laber, 2019), but not for
observational data. When $\pi_{1}$ and $\pi_{2}$ are unknown, they can be
estimated using a logistic regression model. This additional estimation step
will not change the approximation error but the estimation error will likely
change.

##  5 Approximation error

###  5.1 Assumptions

To establish the convergence rate of the approximation error, we require two
assumptions. First, we require the standard assumption that the outcomes are
bounded.

######  Assumption A.

Outcomes $Y_{1},Y_{2}$ satisfy $\max(Y_{1},Y_{2})\leq C_{y}$.

The second assumption is the DTR version of Tsybakov’s small noise assumption
(Tsybakov et al., 2004; Audibert et al., 2007). Before stating the assumption,
we first define the treatment effect contrasts

$\mathcal{T}_{1}(H_{1})=\mathbb{E}[Y_{1}+U_{2}^{*}(H_{2})|A_{1}=1,H_{1}]-\mathbb{E}[Y_{1}+U_{2}^{*}(H_{2})|A_{1}=-1,H_{1}],$
(14)

and

$\mathcal{T}_{2}(H_{2})=\mathbb{E}[Y_{1}+Y_{2}|A_{1},A_{2}=1,H_{2}]-\mathbb{E}[Y_{1}+Y_{2}|A_{1},A_{2}=-1,H_{2}],$
(15)

where $U_{2}^{*}$ is as defined in (2). The above quantities are called the
optimal-blip to-zero function, or sometimes simply the blip function, in the
literature (Robins, 2004; Schulte et al., 2014; Luedtke and Van Der Laan,
2016). Because we assume $Y_{1}$ and $Y_{2}$ are bounded away from zero,
$\eta_{t}(H_{t})=1/2$ if and only if $\mathcal{T}_{t}(H_{t})=0$ for $t=1,2$.
Therefore, the treatment boundary $\\{h_{t}:\eta_{t}(h_{t})=1/2\\}$ can also
be formulated as $\\{h_{t}:\mathcal{T}_{t}(h_{t})=0\\}$.

In classification literature, it is well-noted that obtaining a fast rate of
convergence (faster than $n^{-1/2}$) requires control on the distribution of
the random variable $\eta(X)-1/2$ near the decision boundary
$\\{x:\eta(x)=1/2\\}$ to some degree, which gives rise to the so-called margin
conditions (Tsybakov et al., 2004; Audibert et al., 2007). Similarly, in DTR
context, even with regression based methods, regulation near the treatment
effect boundary $\\{H_{t}:\mathcal{T}_{t}(H_{t})=0\\}$ are generally required;
see Appendix 12.3 for a detailed discussion. Thus it is expected that we too
would require control on the rate of decay of $\eta_{1}-1/2$ and
$\eta_{2}-1/2$ near the treatment boundary $\\{h_{1}:\eta_{1}(h_{1})=1/2\\}$
and $\\{h_{2}:\eta_{2}(h_{2})=1/2\\}$, respectively, to obtain sharp bound on
the $\psi$-regret. Among many variants of margin condition, we consider the
Tsybakov small noise condition (Assumption MA of Audibert et al. (2007); see
also Proposition 1 of Tsybakov et al. (2004)), which has seen wide use in the
literature (Audibert et al., 2007; Steinwart et al., 2007; Blanchard et al.,
2008). The DTR formulation of Tsybakov small noise condition takes the
following form:

######  Assumption B (Tsybakov small noise assumption).

There exist a constant $C>0$, a small number $t_{0}\in(0,1)$, and positive
reals $\alpha_{1},\alpha_{2}$ such that

$P(0<|\eta_{1}(H_{1})-1/2|\leq t)\leq Ct^{\alpha_{1}},\quad
P(0<|\eta_{2}(H_{2})-1/2|\leq t)\leq Ct^{\alpha_{2}}$

for all $t<t_{0}$.

The parameters $\alpha_{1}$ and $\alpha_{2}$ are the _Tsybakov noise
exponents_. We already noted that the $Y_{i}$’s are bounded below. Since the
outcomes are also bounded above by Assumption A, Assumption B is equivalent to
saying

$P(0<\mathcal{T}_{1}(H_{1})<t)+P(0<\mathcal{T}_{2}(H_{2})<t)\leq
Ct^{\alpha},\quad\text{ for all }t\leq t_{0}.$ (16)

This alternative version is more common in precision medicine literature (Qian
and Murphy, 2011; Luedtke and Van Der Laan, 2016). See Supplement 12.3 for
more details on small noise assumption or similar assumptions in precision
medicine literature. Finally, observe that if a stage has Tsybakov noise
exponent $\alpha$, then it also has noise exponent $\alpha^{\prime}$ for all
$\alpha^{\prime}<\alpha$. Thus, to keep our calculations short, we will assume
that both stages have noise exponent $\alpha$ where
$\alpha=\min(\alpha_{1},\alpha_{2})$. We postpone further discussion on
Assumption B till Supplement 12.3.

Under Assumption B, it turns out that, our surrogates satisfying Condition 2
do not exhibit identical approximation error. The difference in the rate stems
from the difference in their respective derivatives. Thus, it will be
convenient to split the above-mentioned surrogates into two types.

######  Definition 2.

We say a surrogate $\phi$ satisfying Condition 2 is of type A if there exists
a constant $B_{\phi}>0$ and $\kappa\geq 2$ such that
$|\phi^{\prime}(x)|<B_{\phi}(1+|x|)^{-\kappa}$ for all $x\neq 0$. We say a
surrogate $\phi$ satisfying Condition 2 is of type B if there exists a
constant $B_{\phi}>0$ and $\kappa>0$ such that
$|\phi^{\prime}(x)|<B_{\phi}\exp(-\kappa|x|)$ for all $x\neq 0$.

$\phi(x)$ | Type | $B_{\phi}$ | $\kappa$
---|---|---|---
(a) $x/(1+|x|)+1$ | A | 1 | 2
(b) $\frac{2}{\pi}\arctan(\pi x/2)+1$ | A | 2 | 2
(c) $x/\sqrt{1+x^{2}}+1$ | A | $2^{3/2}$ | $3$
(d) $2(1+e^{-x})^{-1}$ | B | 2 | $1$
Table 1: Caption

First, Definition 2 assumes $\phi$ to be smooth everywhere except perhaps at
the origin. This restriction rules out non-smooth $\phi$’s, but they are
uninteresting from our implementation perspective anyways. All the $\phi$’s we
have considered in Example 1 are differentiable at $\mathbb{R}/\\{0\\}$ (see
Table 1; more details can be found in Supplement 24). Second, type A merely
means $\phi^{\prime}$ decays polynomially in $|x|^{-\kappa}$, where type B
$\phi$’s enjoy exponential decay of the derivative.

###  5.2 Approximation Error Rate

Theorem 4 summarizes the approximation error rate.

######  Theorem 4.

Suppose $\mathbb{P}$ satisfies Assumptions I-III, Assumption A and Assumption
B with small noise coefficient $\alpha$. Let $0<a_{n}\to\infty$ be any
sequence of positive reals. Further suppose there exist a small number
$\delta_{n}\in(0,1)$ and maps
$\tilde{h}_{n,1}:\mathcal{H}_{1}\mapsto\mathbb{R}$ and
$\tilde{h}_{n,2}:\mathcal{H}_{2}\times\\{0,1\\}\mapsto\mathbb{R}$ so that

$\|\tilde{h}_{n,1}-(\eta_{1}-1/2)\|_{\infty}+\|\tilde{h}_{n,2}-(\eta_{2}-1/2)\|_{\infty}\leq\delta_{n}$

where $\eta_{1}$ and $\eta_{2}$ are defined in (7) and (8). Then for any
$\phi$ of type A, the following holds for any $\alpha^{\prime}\in(0,\alpha)$
satisfying $\alpha-\alpha^{\prime}<1$:

$\displaystyle
V_{\psi}^{*}-V_{\psi}(a_{n}\tilde{h}_{n,1},a_{n}\tilde{h}_{n,2})$
$\displaystyle\lesssim$
$\displaystyle\begin{cases}a_{n}^{1-\kappa}+\min(\delta_{n}^{2+\alpha}a_{n},\delta_{n}^{1+\alpha})&\text{if
}\kappa<2+\alpha\\\
\frac{a_{n}^{-\frac{1+\alpha}{1+(\alpha-\alpha^{\prime})/(\kappa-1)}}}{\alpha-\alpha^{\prime}}+\min(\delta_{n}^{2+\alpha}a_{n},\delta_{n}^{1+\alpha})+\frac{\delta_{n}^{\alpha^{\prime}+2-\kappa}}{(\alpha-\alpha^{\prime})a_{n}^{\kappa-1}}&\text{if
}\kappa\geq 2+\alpha.\end{cases}$

Suppose $\phi$ is of type B. Then

$V_{\psi}^{*}-V_{\psi}(a_{n}\tilde{h}_{n,1},a_{n}\tilde{h}_{n,2})\lesssim\frac{(\log
a_{n})^{1+\alpha}}{a_{n}^{1+\alpha}}+\min(a_{n}\delta_{n}^{2+\alpha},\delta_{n}^{1+\alpha})+a_{n}\delta_{n}\exp(-\kappa
a_{n}\delta_{n}/2).$

Theorem 4 entails that if $\\{\tilde{h}_{n,1},\tilde{h}_{n,2}\\}$ approximates
$\\{\eta_{1}-1/2,\eta_{2}-1/2\\}$ well in the sup-norm, then their scaled
versions $\tilde{f}_{n,1}=a_{n}\tilde{h}_{n,1}$ and
$\tilde{f}_{n,2}=a_{n}\tilde{h}_{n,2}$ incur small regret. It may appear a bit
unusual in that we require $\tilde{f}_{n,i}$’s to be close to the functions
$\eta_{i}-1/2$’s, where $V_{\psi}$ is actually maximized at
$(\tilde{f}_{1},\tilde{f}_{2})$ (see Lemma 15). To that end, note that the
extended real valued functions $\tilde{f}_{i}$’s can not be approximated by
any real valued $f_{i}$’s because $\|f_{i}-\tilde{f}_{i}\|_{\infty}$ is
infinity for all such $f_{i}$’s. However, the proof of Theorem 4 ensures that
$a_{n}(\eta_{i}-1/2)$’s are good proxy for the $\tilde{f}_{i}$’s because

$V_{\psi}(\tilde{f}_{1},\tilde{f}_{2})-V_{\psi}\Big{(}a_{n}(\eta_{1}-1/2),a_{n}(\eta_{2}-1/2)\Big{)}$

is small. The bounds in Theorem 4 holds for any small $\delta_{n}$, whose
optimal rate will be found during from from the the estimation error
calculation.

###  5.3 Special case: strong separation with $\alpha=\infty$

We next describe the convergence rates under a special case of
$\alpha=\infty$, which will be referred as the _strong separation_ :

######  Assumption C (Strong separation).

$\eta_{1}$ and $\eta_{2}$ are bounded away from zero on their respective
domains.

Under this setting, the treatment effects are uniformly bounded away from
zero. See Section 12.3 for related discussion. This special case is of great
interest in the current literature, cf. Zhao et al. (2012, 2015); Qian and
Murphy (2011). When the strong separation holds, the approximation error can
be made much smaller than that of Theorem 4 as detailed in Proposition 1.

######  Proposition 1.

Suppose the conditions of Theorem 4 hold except $\mathbb{P}$ satisfies
Assumption C instead of Assumption B. Suppose
$\tilde{h}_{n,1}:\mathcal{H}_{1}\mapsto\mathbb{R}$ and
$\tilde{h}_{n,2}:\mathcal{H}_{2}\times\\{0,1\\}\mapsto\mathbb{R}$ satisfy

$\|\tilde{h}_{n,1}-(\eta_{1}-1/2)\|_{\infty}+\|\tilde{h}_{n,2}-(\eta_{2}-1/2)\|_{\infty}\leq
c$

for some $c>0$. Then there exists $C>0$ so that the following assertion holds
for any large positive number $a_{n}$:

$\displaystyle
V_{\psi}^{*}-V_{\psi}(a_{n}\tilde{h}_{n,1},a_{n}\tilde{h}_{n,2})\lesssim\begin{cases}a_{n}^{-(\kappa-1)}&\text{if
}\phi\text{ is of type A,}\\\ \exp(-\kappa Ca_{n})&\text{if }\phi\text{ is of
type B}.\end{cases}$

######  Remark 7.

Proposition 1 implies that under the strong separation assumption,
$\tilde{h}_{n,1}$ and $\tilde{h}_{n,2}$ do not even need to approximate
$\eta_{1}-1/2$ and $\eta_{2}-1/2$ very precisely and the approximation error
decays substantially faster. To get a sense of how fast the regret diminishes,
we consider the case when $a_{n}=n$. In this case, the regret is $O_{p}(1/n)$
for a $\phi$ of type A because $\kappa\geq 2$ for a type A $\phi$. The regret
decays exponentially fast for a type B $\phi$, which can be attributed to the
exponential decay of its derivatives. Recall that the derivative of the zero-
one loss function is exactly zero at any $x\neq 0$. Thus the type B surrogates
bear closer resemblance to the original zero-one loss function.

Theorem 4 or Proposition 1 bound the approximation error because if
$\mathcal{U}_{n}=\mathcal{U}_{1n}\times\mathcal{U}_{2n}$ is such that

$\displaystyle\inf_{f_{t}\in\mathcal{U}_{tn}}\|f_{t}-(\eta_{t}-1/2)\|_{\infty}<\delta_{n},\quad
t=1,2,$ (17)

then Theorem 4 or Proposition 1 upper bound
$V_{\psi}^{*}-\sup_{(f_{1},f_{2})\in\mathcal{U}_{n}}V_{\psi}(f_{1},f_{2})$.

##  6 Estimation error

In this section, we focus on the estimation error in (4.1), and provide sharp
regret-bound for a selected set of classifiers by combining all sources of
error. We assume that
$(\widehat{f}_{n,1},\widehat{f}_{n,2})\in\mathcal{U}_{n}=\mathcal{U}_{1n}\times\mathcal{U}_{2n}$,
where $\mathcal{U}_{1n},\mathcal{U}_{2n}$ are classes of functions. Our
analysis in this section is fully nonparametric because our $\mathcal{U}_{n}$
is agnostic of the underlying data generating mechanism. We first present some
theorems (Theorems 5, 6 and 7) for general $\mathcal{U}_{n}$’s. Then we will
move to study the particular examples of neural networks and wavelets.

####  6.0.1 Estimation error when $\mathcal{U}_{n}$ is a general function-
class

For general function-classes, we need some assumptions to control the
complexity of $\mathcal{U}_{n}$. Such assumptions are widely used for bounding
the expectation of the estimation error (Koltchinskii, 2011; Bartlett et al.,
2006; Bartlett and Mendelson, 2002). To define complexity in context of
function-classes, we need to introduce the concept of the bracketing entropy.
Given two functions $f_{l}$ and $f_{u}$, the bracket $[f_{l},f_{u}]$ is the
set of all function $f$ satisfying $f_{l}\leq f\leq f_{u}$. Suppose
$\|\cdot\|$ is a norm on the function-space and $\epsilon>0$. Then
$[f_{l},f_{u}]$ is called an $\epsilon$-bracket if $\|f_{u}-f_{l}\|<\epsilon$.
For a function-class $\mathcal{G}$, we define the bracketing entropy $N_{[\
]}(\epsilon,\mathcal{G},\|\cdot\|)$ to be the minimum number of
$\epsilon$-brackets needed to cover $\mathcal{G}$. This is a measure of the
complexity of $\mathcal{G}$. We will see that the estimation error directly
depends on the bracketing entropy of $\mathcal{U}_{n}$. We will first consider
the case when just the small noise assumption (Assumption B) holds, and then
we will move to the special case when the strong separation (Assumption C)
holds.

We derive the estimation error of DTRESLO under the small noise assumption
(Assumption B) when

$N_{[\
]}(\epsilon,\mathcal{U}_{tn},\|\cdot\|_{\infty})\lesssim\bigg{(}\frac{A_{n}}{\epsilon}\bigg{)}^{\rho_{n}},\quad
t=1,2$ (18)

where $A_{n},\rho_{n}>0$. This leads to the regret bound of Theorem 5 that
depends on $A_{n}$ and $\rho_{n}$. The $\mathcal{U}_{n}$’s that satisfy (18)
are called VC-type classes (p. 41 Koltchinskii, 2011). In all our examples,
$\mathcal{U}_{n}$ will satisfy (18) for appropriate $A_{n}$ and $\rho_{n}$.

######  Theorem 5.

Suppose $\mathcal{U}_{n}$ is such that there exists $A_{n}>0$ and
$\rho_{n}\in\mathbb{R}$ so that (18) holds with $\liminf_{n}\rho_{n}>0$,
$\rho_{n}\log A_{n}=o(n)$, and $\liminf_{n}\rho_{n}\log A_{n}>0$. Further
suppose there exist $(\tilde{f}_{n,1},\tilde{f}_{n,2})\in\mathcal{U}_{n}$ so
that

$\|\tilde{f}_{n,1}/a_{n}-(\eta_{1}-1/2)\|_{\infty}+\|\tilde{f}_{n,2}/a_{n}-(\eta_{2}-1/2)\|_{\infty}\leq\bigg{(}\frac{\rho_{n}\log
A_{n}}{n}\bigg{)}^{1/(2+\alpha)}$ (19)

for some $a_{n}=n^{a}$ where $a>1$. We also assume that $\mathbb{P}$ satisfies
Assumptions I-III, Assumption A, and Assumption B with coefficient $\alpha>0$.
Then there exist $C>0$ and $N_{0}\geq 1$ such that for all $n\geq N_{0}$ and
all $x>0$,

$V_{\psi}^{*}-V_{\psi}(\widehat{f}_{n,1},\widehat{f}_{n,2})\leq
C\max\bigg{\\{}(1+x)^{2}(\log n)^{2}\bigg{(}\frac{\rho_{n}\log
A_{n}}{n}\bigg{)}^{\frac{1+\alpha}{2+\alpha}},\texttt{Opt}_{n}\bigg{\\}}$

with probability at least $1-\exp(-x)$.

Theorem 5 is proved in Section 20. In our examples with neural networks and
wavelets, we will see that the regret bound in Theorem 5 leads to sharp rates
provided $\eta_{1}$ and $\eta_{2}$ satisfy some smoothness conditions.

Next we will derive regret bounds under the special case of strong separation
(Assumption C). We begin with a theorem that consider $\mathcal{U}_{n}$’s
satisfying the following entropy bound:

$\displaystyle\log N_{[\
]}(\epsilon,\mathcal{U}_{tn},L_{2}(\mathbb{P}_{n}))\lesssim\bigg{(}\frac{A_{n}}{\epsilon}\bigg{)}^{2\rho_{n}},\quad
t=1,2.$ (20)

General Hölder classes satisfy (20) (cf. p.154, Van Der Vaart et al., 1996).

######  Theorem 6.

Suppose the function-class $\mathcal{U}_{n}$ is such that (20) holds with
$\rho_{n}\in(0,1)$, $A_{n}>1$, where it also holds that
${A_{n}^{2\rho_{n}}}/{n}\to 0$. Further suppose the approximation error

$V_{\psi}^{*}-\sup_{(f_{1},f_{2})\in\mathcal{U}_{n}}V_{\psi}(f_{1},f_{2})=O((\log
n)^{k}/{n})$

for some $k\in\mathbb{N}$, the optimization error $\texttt{Opt}_{n}<1/2$, and
Assumptions I-III, A, and C hold. Then there exist $C>0$ and $N_{0}>0$ such
that for all $n\geq N_{0}$ and for any $x>0$,

$V_{\psi}^{*}-V_{\psi}(\widehat{f}_{n,1},\widehat{f}_{n,2})\leq
C\max\bigg{\\{}(1+x)^{2/(1+\rho_{n})}{A_{n}^{2\rho_{n}/(1+\rho_{n})}}{n^{-1/(1+\rho_{n})}},\
\texttt{Opt}_{n}\bigg{\\}}$

with probability at least $1-\exp(-x)$.

Theorem 6 is proved in Supplement 22.1. The proof of Theorem 6 shows that
under (20), the estimation error is larger than the approximation error. Thus
in this case, the biggest contribution in the regret comes from the estimation
error term, which is of the order
$O_{p}\bigg{(}A_{n}^{2\rho_{n}/(1+\rho_{n})}n^{-1/(1+\rho_{n})}\bigg{)}$. This
type of rate is also observed for the $L_{2}$ risk in regression problems
where the regression function is a member of a class with similar entropy
bounds (p. 75, Example 4, Koltchinskii, 2011). When $\rho_{n}\to 0$, the
complexity of the class decreases, and as a result, the estimation error
decreases as well.

Our next theorem, which is proved in Section 22.2, considers VC-type
$\mathcal{U}_{n}$’s similar to Theorem 5. The corresponding entropy bound is
smaller than that in (20), which results in smaller regret bound compared to
Theorem 6.

######  Theorem 7.

Suppose $\mathcal{U}_{n}$ is a function-class such that there exist $A_{n}>0$
and $\rho_{n}>0$ so that

$N_{[\
]}(\epsilon,\mathcal{U}_{n},L_{2}(\mathbb{P}_{n}))\lesssim\bigg{(}\frac{A_{n}}{\epsilon}\bigg{)}^{\rho_{n}},$
(21)

and $\liminf_{n}(\rho_{n}\log A_{n})>0$. Further suppose the approximation
error

$V_{\psi}^{*}-\sup_{(f_{1},f_{2})\in\mathcal{U}_{n}}V_{\psi}(f_{1},f_{2})=O((\log
n)^{k}/{n})$ (22)

for some $k\in\mathbb{N}$. Then under Assumptions I-III, A, and C, there exist
$C>0$ and $N_{0}>0$ such that for all $n>N_{0}$ and any $x>0$,

$V_{\psi}^{*}-V_{\psi}(\widehat{f}_{n,1},\widehat{f}_{n,2})\leq
C\max\bigg{\\{}\frac{(1+x)^{2}(\log n)^{2}(\rho_{n}\log A_{n})^{2}+(\log
n)^{k}}{n},\ \texttt{Opt}_{n}\bigg{\\}}$

with $\mathbb{P}$-probability at least $1-\exp(-x)$.

We will see that in our neural network example, $A_{n}$ is a polynomial in
$n$, and $\rho_{n}$ is a poly-log term, which leads to a regret of order
$O_{p}(1/n)$ upto a polylog term. In binary classification context, many
methods are known to attain this sharp rate (Massart et al., 2006; Blanchard
et al., 2008). However, to the best of our knowledge, our DTRESLO method is
the first DTR method for which such a sharp regret bound is established under
our nonparametric set-up. Theorem S.1 of Zhao et al. (2015) implies that if
$\eta_{1}$ and $\eta_{2}$ are uniformly bounded away from both zero and one,
then the regret of BOWL and SOWL estimators can be pushed to the order of
$O_{p}(n^{-\omega/(1+\omega)})$ under Assumption C, where $\omega>0$ is a
fixed quantity, called geometric noise exponent, that depends only on
$\mathbb{P}$. This rate is slower than our rate when the optimization error
$\texttt{Opt}_{n}$ is of the order $O_{p}(1/n)$. Theorem 3.1 of Qian and
Murphy (2011) can be used to obtain regret bounds for Q-learning type methods
under Assumption C and D. However, the resulting regret bound has the same
order as the $L_{2}(\mathbb{P})$ estimation error of the Q-functions (cf.
equation 3.6 of Qian and Murphy, 2011), which can not be faster than
$n^{-2\beta/(2\beta+p)}$ under nonparametric set-up (Yang, 1999).

###  6.1 Examples of regret bounds with specific $\mathcal{U}_{n}$’s

We will first state some assumptions on $\eta_{1}$ and $\eta_{2}$ that ensure
the approximability of $G_{1}$ and $G_{2}$ by $\mathcal{U}_{n}$ when
$\mathcal{U}_{n}$ corresponds to basis-expansion type classes. Next we will
elaborate on the special cases when $\mathcal{U}_{n}$ corresponds to neural
networks and wavelets classes.

####  6.1.1 Smoothness assumption

First we explain why smoothness conditions on $\eta_{1}$ and $\eta_{2}$ are
required. To control the estimation error, some structures on
$\mathcal{U}_{n}$ are desirable because smaller search spaces for the
$f_{t}$’s results in smaller estimation error. Restricting the search space is
equivalent to restricting the complexity of the class $\mathcal{U}_{n}$
(Audibert et al., 2007). We require structural assumptions on $\eta_{1}$ and
$\eta_{2}$ to ensure that $\eta_{1}$ and $\eta_{2}$ are well approximable by
such $\mathcal{U}_{n}$’s – giving rise to the so-called complexity
assumptions. Thus, the complexity assumption enables the attainment of small
estimation error without necessarily blowing up the approximation error. See
Audibert et al. (2007), Koltchinskii (2011) and Tsybakov et al. (2004), among
others, for more detailed account on the necessity of complexity assumptions.
In classification context, $\mathcal{U}_{n}$ is taken to be some smoothness
class or VC class because most popular classifiers, e.g. neural network,
basis-expansion type classifiers, wavelets etc. belong to such classes. Such
classes can approximate $\eta_{1}$ and $\eta_{2}$ well if the latter are
smooth. Therefore, we will assume our $\eta_{1}$ and $\eta_{2}$ belong to
smoothness classes. To that end, we define the Hölder classes with smoothness
index $\beta>0$ below.

Let $p\in\mathbb{N}$. A function
$f:\mathcal{\mathcal{X}}\subset\mathbb{R}^{p}\mapsto\mathbb{R}$ is said to
have Hölder smoothness index $\beta>0$ if for all
$u=(u_{1},\ldots,u_{p})\in\mathbb{N}^{p}$ satisfying $|u|_{1}<\beta$,
$\partial^{u}f=\partial^{u_{1}}\partial^{u_{2}}\ldots\partial^{u_{p}}f$ exists
and there exists a constant $C>0$ so that

$\frac{|\partial^{u}f(x)-\partial^{u}f(y)|}{|x-y|^{\beta-\left\lfloor\beta\right\rfloor}}<C\quad\text{
for all }x,y\in\mathcal{X}.$

For some $\mathcal{Y}>0$, we denote by
$\mathcal{C}^{\beta}_{d}(\mathcal{X},\mathcal{Y})$ the Hölder class of
functions given by

$\displaystyle\bigg{\\{}f:\mathcal{X}\subset\mathbb{R}^{d}\mapsto\mathbb{R}\
\bigg{|}\
\sum_{u:|u|_{1}<\beta}\|\partial^{u}f\|_{\infty}+\sum_{u:|u|=\left\lfloor\beta\right\rfloor}\sup_{\begin{subarray}{c}x,y\in\mathcal{\mathcal{X}}\\\
x\neq
y\end{subarray}}\frac{|\partial^{u}f(x)-\partial^{u}f(y)|}{|x-y|^{\beta-\left\lfloor\beta\right\rfloor}}\leq\mathcal{Y}\bigg{\\}}.$
(23)

Since $H_{t}$ may include categorical variables such as smoking status, we
separate the continuous and categorical parts of $H_{t}$ as
$H_{t}=(H_{ts},H_{tc})\in\mathcal{H}_{t}=\mathcal{H}_{ts}\otimes\mathcal{H}_{tc}$,
where $H_{ts}\in\mathcal{H}_{ts}\subset\mathbb{R}^{p_{ts}}$ and
$H_{tc}\in\mathcal{H}_{tc}\subset\mathbb{R}^{p_{tc}}$ correspond to the
continuous and categorical part of $H_{t}$, for $t=1,2$.

######  Assumption D (Smoothness assumption).

$\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ are compact and $\mathcal{H}_{1c}$ and
$\mathcal{H}_{2c}$ are finite sets. Also, there exist $\beta>0$ and
$\mathcal{Y}>0$ so that the following hold:

  * 1.

Let $\mathcal{X}=\mathcal{H}_{1s}$. For each $h\in\mathcal{H}_{1c}$, the map
$\eta_{1}(\cdot,h):\mathcal{X}\mapsto\mathbb{R}$ is in
$\mathcal{C}^{\beta}_{p_{1s}}(\mathcal{X},\mathcal{Y})$.

  * 2.

Let $\mathcal{X}=\mathcal{H}_{2s}$. For each
$(h,a)\in\mathcal{H}_{2c}\times\\{\pm 1\\}$, the function
$\eta_{2}(\cdot,h,a):\mathcal{X}\mapsto\mathbb{R}$ is in
$\mathcal{C}^{\beta}_{p_{2s}}(\mathcal{X},\mathcal{Y})$.

We formulated the smoothness assumption in terms of $\eta_{1}$ and $\eta_{2}$
so that our results are consistent with contemporary classification
literature. However, our proofs show that one could formulate the assumptions
in terms of the smoothness of the blip functions in (14) and (15) as well. The
compact support assumption for $H_{t}$, which is typically satisfied in real
applications, is also commonly required in the DTR literature (Zhao et al.,
2012, 2015; Sonabend et al., 2021; Zhang et al., 2018b). Under the compactness
assumption, $\mathcal{H}_{1c}$ and $\mathcal{H}_{2c}$ are finite sets.
Smoothness conditions as Assumption D have appeared in DTR literature in
context of nonparametric estimation (Sun and Wang, 2021). Compared to the
parametric assumptions often imposed on the blip functions in Q-learning or
A-learning (Schulte et al., 2014), our smoothness assumptions are much weaker.
Our smoothness assumption includes non-differentiable functions as well. Next,
we will establish regret bounds for neural network and wavelets classes.

####  6.1.2 Neural networks as an example of $\mathcal{U}_{n}$

We consider the neural network space in line with Schmidt-Hieber (2020)’s
construction. Let ${\mathcal{F}}(L,W,s,\mathcal{Y})$ be the class of ReLU
networks uniformly bounded by $\mathcal{Y}>0$, with depth $L\in\mathbb{N}$,
width vector $W$, sparsity $s\in\mathbb{N}$, and weights bounded by one. The
output layer of the networks in ${\mathcal{F}}(L,W,s,\mathcal{Y})$ uses a
linear gate. In this example, we consider that for $t=1,2$, the class
$\mathcal{U}_{tn}$ corresponds to
${\mathcal{F}}(L_{n},W_{n},s_{n},\mathcal{Y}_{n})$ where $L_{n}$, $W_{n}$,
$s_{n}$, and $\mathcal{Y}_{n}$ may depend on $n$. To avoid cumbersome
notation, we drop $n$ from $L_{n}$, $W_{n}$, $s_{n}$, and $\mathcal{Y}_{n}$,
and simply denote them by $L$, $p$, $s$, and $\mathcal{Y}$, respectively. One
can control the complexity of this class via pre-specifying upper bounds on
the depth, width, and sparsity of the network. We will first consider regret
bound under Assumption B, and then we move to the special case of strong
separation (Assumption C).

Corollary 1 establishes the regret bound of DTRESLO with neural network
classifier under Assumption B.

######  Corollary 1.

Suppose $\mathbb{P}$ satisfies Assumptions I-III, Assumption A, Assumption B
with parameter $\alpha>0$, and Assumption D with parameter $\beta>0$. Let
$\mathcal{U}_{n,1}$ and $\mathcal{U}_{n,2}$ be of the form
$\mathcal{F}(L,W,s,\infty)$ with appropriate $W_{1}$, where
$\mathcal{F}(L,W,s,\infty)$ is as defined in Section 6.1.2. Suppose
$L=c_{1}\log n$, $s=c_{2}n^{p/((2+\alpha)\beta+p)}$, and the maximal width
$\max W\leq c_{3}s/L$ where $c_{1},c_{2},c_{3}>0$. Then there exist $N_{0}>0$
and $C>0$ depending on $\mathbb{P}$ and $\psi$ such that if
$c_{1},c_{2},c_{3}>C$, then for $n\geq N_{0}$ and any $x>0$, the following
holds with probability at least $1-\exp(-x)$:

$V^{*}_{\psi}-V_{\psi}(\widehat{f}_{n,1},\widehat{f}_{n,2})\leq
C\max\bigg{\\{}(1+x)^{2}(\log
n)^{\frac{6+4\alpha}{2+\alpha}}n^{-\frac{1+\alpha}{2+\alpha+p/\beta}},\
\texttt{Opt}_{n}\bigg{\\}}.$

The proof of Corollary 1 can be found in Section 21.1. The proof of Corollary
1 assumes that $p$ is fixed, i.e. it does not grow with $n$. The generic
constant $C$ in Corollary 1 may depend on $p$ as well. Under Assumptions
similar to A, B, and D, the rate $n^{-\frac{1+\alpha}{2+\alpha+p/\beta}}$ is
minimax in context of binary classification (Audibert et al., 2007). Since
two-stage weighted classification problem is not easier than binary
classification, this rate is expected to be the minimax rate under our set-up
as well. To the best of our knowledge, no other nonparametric DTR method has
better guarantees for the regret under set-up similar to ours. See Section 7
for a comparison of the regret bound of our DTRESLO method with some other
existing methods.

######  Corollary 2.

Suppose $\mathbb{P}$ satisfies Assumptions I-III, Assumption A with
$\mathcal{B}>0$, Assumption C, and Assumption D with parameter $\beta>0$. Let
$\mathcal{U}_{n,1}$ and $\mathcal{U}_{n,2}$ be of the form
$\mathcal{F}(L,W,s,\infty)$ with appropriate $W_{1}$, where
$\mathcal{F}(L,W,s,\infty)$ is as defined in Section 6.1.2. Suppose
$L=c_{1}\log n$, $s=c_{2}(\log n)^{p/\beta}$, and the maximal width $\max
W\leq c_{3}\log n$ where $c_{1},c_{2},c_{3}>0$. Then there exist
$C,C^{\prime}>0$ depending on $p$, $\beta$, and $\psi$ such that if
$c_{1},c_{2},c_{3}>C$, then for all $n\geq N_{0}$ and any $x>0$,

$V_{\psi}^{*}-V_{\psi}(\widehat{f}_{n,1},\widehat{f}_{n,2})\leq
C^{\prime}\max\bigg{\\{}\frac{(1+x)^{2}(\log n)^{3+p/\beta}}{n},\
\texttt{Opt}_{n}\bigg{\\}}$

with probability at least $1-\exp(-x)$.

The proof of Corollary 2 can be found in Section 21.2. Corollary 2 implies
that DTRESLO method with neural network classifier can attain regret
$O_{p}(1/n)$ upto a poly-log term, as announced earlier.

####  6.1.3 Wavelets as an example of $\mathcal{U}_{n}$

In this section, we will define the wavelet estimators that we will use. Let
us consider father and mother wavelets $\nu$ and $\zeta$. We will assume that
these wavelets are S-regular wavelets (cf. Definition 4.2.14, p.326, Giné and
Nickl, 2015) with $S>\beta$ and they have compact support, e.g., Daubechies
wavelets (p. 318, Giné and Nickl, 2015). The corresponding $p$-dimensional
wavelets are constructed using tensor products

$\Gamma({x})=\nu(x_{1})\cdots\nu(x_{p}),\quad{x}=(x_{1},\dots,x_{p})\in K,$

where $K\subset\mathbb{R}^{p}$ is a compact set. Define
$\Gamma_{k}({x})=\Gamma({x}-{k})$ for $k\in\mathbb{Z}$, where $\mathbb{Z}$ was
denoted to be the set of all integers. Observe that since $\nu$ has compact
support, $\nu(\cdot-k)$ restricted to any compact set is non-zero only for
finitely many $k$’s. Therefore, $\Gamma_{k}$ restricted to compact set $K$ is
non-zero only for finitely many $k$’s. We denote the set of such $k$’s by
$\mathcal{K}_{0}$. For each $i\in\\{0,1\\}^{p}\setminus(0,\ldots 0)$, we also
consider the $p$-dimensional tensor product

$\mathcal{Z}^{i}({x})=\zeta^{i_{1}}(x_{1})\cdots\zeta^{i_{p}}(x_{p}),$

where $\zeta$ is the mother wavelet function. For $x\in K$, we define the
wavelet function
${\mathcal{Z}}^{i}_{lk}({x})=2^{lr/2}{\mathcal{Z}}^{i}(2^{l}{x}-{k})$, where
$l\in\mathbb{N}\cup\\{0\\}$, and $k\in\mathbb{Z}$. Notice that since $\zeta$
has compact support and $K$ is also compact, the function
$x\mapsto{\mathcal{Z}}(2^{l}x-{k})$ is non-zero only for some specific $k$’s
belonging to a finite set. We will call this set $\mathcal{K}(l)$. It is easy
to show that there exist $C_{1}$ and $C_{2}>0$ so that the cardinality of
$\mathcal{K}(l)$ satisfies $C_{1}2^{lp}\leq|\mathcal{K}(l)|<C_{2}2^{lp}$ for
all $l\geq 1$. We let $l\in\\{0,1,\ldots,b_{n}\\}$, where $b_{n}$ is a
sequence of integers diverging to $\infty$. This number $b_{n}$ will be called
the level of the wavelet class. Finally, we define the wavelet function-class
by

$\displaystyle\mathcal{H}^{M}_{n}(K)=\bigg{\\{}f_{n}:K\mapsto\mathbb{R}\
\bigg{|}\ $ $\displaystyle\
f_{n}=\sum_{{k}\in\mathcal{K}_{0}}c_{k}\Gamma_{k}+\sum_{l=0}^{b_{n}}\sum_{{k}\in\mathcal{K}(l),{i}\in\mathcal{I}}c_{l{k}}{\mathcal{Z}}_{l{k}}^{i}\text{
where }c_{k}\in\mathbb{R}$ $\displaystyle\ \text{ for all
}k\in\mathcal{K}_{0}\text{ and }c_{lk}\in\mathbb{R}\text{ for all
}l\in\mathbb{Z}\text{ and }k\in\mathcal{K}(l),$ $\displaystyle\
\sup_{k\in\mathcal{K}_{0}}c_{k}+\sup_{l\geq
0}\Big{\\{}2^{l(\beta+p/2)}\sup_{k\in\mathcal{K}(l)}|c_{lk}|\Big{\\}}<M,$ (24)
$\displaystyle\ \mathcal{I}=\\{0,1\\}^{p}\setminus(0,\ldots,0)\bigg{\\}}.$
(25)

The following two corollaries give the regret bound for the wavelet-based
classifiers. As usual, we assume that $p\ll n$. The following corollary, which
is proved in Section 21.3, indicates that under the small noise condition, the
regret decay rate is similar to that of the neural network example upto a
poly-log term. Therefore, similar to the neural network example, the regret
decay rate matches the minimax rate of risk decay in binary classification
under conditions similar to Assumptions A, B, and D (Audibert et al., 2007).

######  Corollary 3.

Suppose $\\{{\Gamma}_{k}\\}_{k\in\mathcal{K}_{0}}$ and
$\\{\mathcal{Z}_{lk}:{k\in\mathcal{K}(l),l\in\mathbb{N}\cup\\{0\\}}\\}$ form
an compactly supported S-regular wavelet basis where $S>\beta>0$. For $t=1,2$,
let $\mathcal{U}_{tn}$ denote the wavelet class $\mathcal{H}_{n}^{M}(K_{t})$
defined in (6.1.3) with $M=2n$ and level $b_{n}$, where $K_{t}$’s are compact
sets such that $K_{1}\supset\text{dom}(\eta_{1})$ and
$K_{2}\supset\text{dom}(\eta_{2})$. If in addition Assumptions I-III,
Assumption A hold, Assumption B holds with $\alpha>0$, and Assumption D holds
with $\beta>0$, then there exist $C,N_{0}>0$ such that for all $n\geq N_{0}$
and any $x>0$,

$V_{\psi}^{*}-V_{\psi}(\widehat{f}_{n,1},\widehat{f}_{n,2})\leq
C\max\bigg{\\{}(1+x)^{2}(p2^{p})^{\frac{1+\alpha}{2+\alpha}}(\log
n)^{\frac{5+3\alpha}{2+\alpha}}n^{-\frac{1+\alpha}{2+\alpha+p/\beta}},\
\texttt{Opt}_{n}\bigg{\\}}$

with probability at least $1-\exp(-x)$.

The following corollary gives the regret decay rate under Assumption C, the
strong separation assumption. This rate is also similar to the analogous rate
for the neural network example upto a poly-log term when $p\ll n$. Corollary 4
is proved in Section 21.4.

######  Corollary 4.

Suppose $\mathbb{P}$ and $\mathcal{U}_{n}$ are as in Corollary 3. Further
suppose Assumptions I-III, A, C, and D hold with smoothness parameter
$\beta>0$. Then there exist constants $C,C^{\prime},N_{0}>0$ such that if
$n\geq N_{0}$ and $b_{n}>C$, then for any $x>0$,

$V_{\psi}^{*}-V_{\psi}(\widehat{f}_{n,1},\widehat{f}_{n,2})\leq
C^{\prime}\max\bigg{\\{}\frac{(1+x)^{2}2^{2Cp}(\log
n)^{4}}{n},\texttt{Opt}_{n}\bigg{\\}}$

with probability at least $1-\exp(-x)$.

##  7 Connection to related literature

### Comparison with SOWL (Zhao et al., 2015)

As mentioned earlier, our work is inspired by the SOWL method of Zhao et al.
(2015), who use the hinge loss $\psi(x,y)=\min(x-1,y-1,0)$ as the surrogate.
Since hinge loss is non-smooth, the resulting surrogate value function is not
amenable to gradient based optimization methods. The surrogate value function
is therefore optimized via the dual formulation, which leads to SVM type
methods. However, SVM type approaches have several limitations. First, the
number of linear constraints of the SVM formulation scales linearly with the
sample size, which leads to computational complexity of $O(n^{3})$ (Williams
and Seeger, 2001). This makes SVM based methods not scalable for large scale
EHR studies where $n$ is large. Second, the dual objective does not facilitate
straightforward variable selection via the addition of an $l_{1}$ penalty.
Finally, SVM type methods are not flexible enough to accommodate tree-based or
neural-network type classifiers with rigorous statistical guarantees. Our
surrogate loss based approach is flexible, and it has the advantage of both
computational efficiency and the ability to incorporate variable selection.

### Comparison of regret upper bound

The direct search methods are often advocated as advantageous when compared
with regression-based methods because it does not require modeling of nuisance
parameters (Zhao et al., 2012, 2015). We will present some heuristics to show
how our theoretical findings strengthen this claim. Notably, our numerical
finding does not refute this claim either; see Section 8 for more details.
Table 2 compares our regret bound with neural network and wavelets classifiers
with the available regret bound for the following nonparametric DTR methods:
BOWL/SOWL (Zhao et al., 2015), nonparametric Q-learning, the BART-based method
of Zhang et al. (2018b), and the classification list-based method of Sun and
Wang (2021). An important fact about the last two methods is that they
specially target enhanced interpretability. We also remark that the upper
bound on Q-learning regret under small noise assumption is derived using Qian
and Murphy (2011)’s results, which were for stage one. We refer to Supplement
13 for more details on these methods, and the construction of Table 2 using
the results in the associated papers. In Table 2, by “noise assumption”, we
indicate assumptions on the data-distribution at the treatment boundary;
Assumption B is an example but see Section 5.1 for more details. By
“smoothness assumption”, we refer to assumptions on the smoothness of
functions such as $\eta_{1}$, $\eta_{2}$, the blip functions in (14) and (15),
or their derivatives.

It is important to note that the regret bounds of Zhao et al. (2015) are
originally produced under the geometric noise condition of Steinwart et al.
(2007). We use the correspondence between the latter condition and our
Assumptions B and D, which follows from Steinwart et al. (2007) itself. See
Supplement 13 for more details on this. Also, the interpretable methods and
Q-learning do not consider the approximation error, and thus their regrets are
with respect to the best treatment regime within the classes under
consideration.

Method |  Regret bound exponent |  Noise assumption type |  Smoothness assumption type
---|---|---|---
DTRESLO with NN and wavelets | $\min(\frac{1+\alpha}{2+\alpha+p/\beta},-\log_{n}\texttt{Opt}_{n})$ | Small noise | Hölder ($\beta$)
BOWL/ SOWL | $\begin{matrix}\frac{1}{2}\frac{1+\alpha}{1+\alpha+p/\beta}-\epsilon^{*}&\text{if }\alpha\geq p/\beta-1\\\ \frac{1+\alpha}{3+3\alpha+p/\beta}-\epsilon^{*}&\text{o.w.}\end{matrix}$ | Small noise | Hölder ($\beta$)
Nonparametric Q-learning† | $\frac{(1+\alpha)}{2+\alpha}\frac{2}{2+p/\beta}$ | Small noise | Hölder ($\beta$)
ST-RL | $\frac{1}{3}\frac{2}{2+p/\beta}$ | $\text{Different}^{**}$ | Hölder ($\beta$)
List- based | $\Big{(}\frac{2}{3}\Big{)}^{l^{{\dagger}{\dagger}}}\frac{1}{2+p/\beta}$ | $\text{Different}^{**}$ |  Modulus of smoothness $(\beta)^{**}$
Table 2: Best possible regret-decay rate given by available upper bounds on
different DTRs: here the decay rate is $O(n^{-\text{exponent}})$ upto a
polylog term. Also, $\texttt{Opt}_{n}$ is the optimization error mentioned in
Section 4.1. Here NN: neural network.
$*$Here $\epsilon$ is any positive number.
${**}$ See Supplement 13 for more details.
${\dagger}$ This regret bound holds for one stage case.
${\dagger}{\dagger}$ $l$ is the length of the longest classification list.

First note that, since $\alpha>0$, we always have $1/3<(1+\alpha)/(2+\alpha)$.
Thus Table 2 implies ST-RL has a larger regret bound than nonparametric
Q-learning. However, ST-RL is not far from a Bayesian version of nonparametric
Q-learning. The list-based method of Zhang et al. (2018b) also faces such an
apparent loss in efficiency. The less efficiency for these two methods is best
perceived as the cost of enhanced interpretability (Sun and Wang, 2021).
Anyway, since no lower bound on these regrets are available, nothing can be
said for certain.

By elementary algebra, it can be shown that when $\texttt{Opt}_{n}$ is
negligible, our rate is sharper than the rate garuantees of all other methods
in Table 2. For the sake of simplicity, in the rest of the discussion, we
assume that $\texttt{Opt}_{n}$ is negligible. We have already mentioned that
that Zhao et al. (2015)’s original analysis uses the geometric noise
assumption. In Supplement 13, we show that the geometric noise assumption is
satisfied when MA-type small noise assumption and Assumption D both hold. From
Steinwart et al. (2007) it follows that under our assumptions, binary
classification with hinge loss almost attains the regret decay rate
$n^{-(1+\alpha)/(2+\alpha+p/\beta)}$. (We use the word “almost” here because
the term $\epsilon$ in Table 2 still can not be removed.) Had the small noise
condition been used, exploiting the direct correspondence between BOWL and
binary classification, it may be possible to obtain a similar upper bound for
BOWL’s regret. Therefore it is possible that BOWL also attains the regret
decay rate of $n^{-(1+\alpha)/(2+\alpha+p/\beta)-\epsilon}$ under our
assumptions.

We speculate that improving the Q-learning regret bound may not be easy
without imposing stricter assumptions on the Q-function. The reason is that we
anticipate the term $2/(2+p/\beta)$ probably can not be improved. This term
stems from the lower bound on $L_{2}(\mathbb{P})$ error rates of nonparametric
methods; see Supplement 13 for more details. Therefore, we speculate that, at
least under our setting, there may be a gap between the regret of NP
Q-learning type methods and direct search methods such as ours or BOWL. To
rigorously prove this claim, one has to calculate the lower bound on the
regret of two stage NP Q-learning, which is beyond the limits of the current
paper. Nevertheless, our comparisons at least does not refute the claim that
direct search methods are more efficient in principle, and partially explains
Kosorok and Laber (2019)’s comments on direct search methods being more
efficient implementation-wise in settings with smaller stages.

However, it must be kept in mind that the extra efficiency may come at the
cost of computational hardness as discussed in Section 3. Of course BOWL,
being a convex method, does not face theoretical computational barriers, but
it excludes data points in the initial stages and does not scale well to large
size of data. There is also a possibility that if we opt for more
interpretable classifiers, the theoretical efficiency will decline as in the
case of ST-RL and the classification list based method of Zhang et al.
(2018b). Hence, there is a possibility of a trade-off among theoretical
efficiency, computational hardness, and interpretability, at least in
principle. Future research in this direction may be able to shed more light
onto this trade-off.

##  8 Empirical Analysis

We have performed extensive simulations to evaluate the finite sample
properties of our DTRESLO methods and to compare with existing methods under
different generating data mechanisms. We additionally evaluate the performance
of the proposed method using EHR data to identify optimal treatment rules for
patients in intensive care units (ICU) with sepsis.

###  8.1 Simulations

We compare the performance of our DTRESLO method with the regression-based
method Q-learning and the direct search methods BOWL and SOWL (Zhao et al.,
2015). For our DTRESLO method, we take $\phi(x)=1+2/\pi\cdot\arctan(\pi x/2)$
because simulation shows that it has slightly better performance than the
other smooth surrogates considered in Example 1. We consider several choices
for the class of classifiers $\mathcal{U}_{1n}$ and $\mathcal{U}_{2n}$. When
we consider the linear treatment policies, $\mathcal{U}_{1n}$ and
$\mathcal{U}_{2n}$ are the class of all linear functions on $\mathcal{H}_{1}$
and $\mathcal{H}_{2}$, respectively. We consider cubic spline, wavelets, and
neural network (NN) as the non-linear treatment policies, with
$\mathcal{U}_{1n}$ and $\mathcal{U}_{2n}$ being the respective function-
classes in these cases. For the comparators, i.e. Q-learning, BOWL, and SOWL,
we consider both linear and non-linear policies as well. Following Zhao et al.
(2015), we incorporate non-linear policies for BOWL and SOWL using a
reproducing kernel Hilbert space (RKHS) with RBF kernel; see Zhao et al.
(2015) for more details. The non-linear treatment policy for Q-learning is
achieved by letting the Q-functions to be in neural network classes. See
Supplement 14 for more details on the implementation of these methods.

We considered five broad simulation settings as detailed in Section 14 of the
Supplement:

  1. 1.

All covariates are discrete. Hence, exhaustive search over $\mathcal{F}$ is
possible using saturated models.

  2. 2.

This is a setting with non-linear decision boundaries in both stages. However,
$Y_{2}$ does not depend on $A_{1}$.

  3. 3.

This setting is inspired by Setting 2 of Zhao et al. (2015), where the outcome
models i.e. $\mathbb{E}[Y_{t}\mid H_{t}]$’s are linear function of $H_{t}$ for
$t=1,2$. We will call this setting the linear setting.

  4. 4.

This has highly non-linear and even non-smooth decision boundaries.

  5. 5.

This setting has higher number of covariates. In this case, the first stage
outcome model is linear but the second stage outcome model is non-linear.

Setting 1 is a simple toy setting. The motivation behind including this
setting is verifying the consistency of DTRESLO. We will use the linear
setting 3 to check if linear treatment policies perform well when the outcome
models are linear. On the other hand, we include settings 2 and 4 to examine
if the methods with non-linear treatment policies have an edge over the
methods with linear treatment policies when the decision boundaries are non-
linear. Finally, setting 5 is included to compare the performance of different
methods when the dimension of $\mathcal{O}_{1}\cup\mathcal{O}_{2}$ is
comparatively larger.

Under each listed setting, we estimate the DTRs based on samples of size
$n=250,2500,5000$. For each estimated DTR $\widehat{d}$, we estimate the value
function $V(\widehat{d}_{1},\widehat{d}_{2})$ by the empirical value function
based on an independent sample of size 10,000. We estimate the expectation and
the standard deviation of these value function estimates using 500 Monte Carlo
replications. We also estimate the optimal value function
$V^{*}=V(d_{1}^{*}(H_{1}),d_{2}^{*}(H_{2}))$ for each setting using these 500
Monte Carlo replications. Figures 4 and 5 compare the estimated expected value
functions of the different methods under consideration. In these figures, we
use the neural network DTRESLO as the non-linear DTRESLO because this method
is comparable to neural network Q-learning. The average value functions
corresponding to the other non-linear DTRESLO methods can be found in Table 5
in Supplement 14. The overall performance of all the non-linear DTRESLO
methods are quite similar although NN DTRESLO is slightly better than the
rest.

(a) Setting 1.

(b) Setting 2.

(c) Setting 4.

Figure 4: Plot of the estimated average value functions for the settings 1, 2
and 4. Here the black horizontal line corresponds to the true value function.
The left and right panels correspond to the linear and non-linear treatment
policies, respectively. Here the non-linear DTRESLO corresponds to the neural
network classifier. The error bars are given by $\pm 2$ SD.

(a) Setting 3.

(b) Setting 5.

Figure 5: Plot of the estimated average value functions for the settings 3 and
5. Here the black horizontal line corresponds to the true value function. The
left and right panels correspond to the linear and non-linear treatment
policies, respectively. Here the non-linear DTRESLO corresponds to the neural
network classifier. The error bars are given by $\pm 2$ SD.

First of all, Figures 4 and 5 entail that DTRESLO consistently performs better
or at least as good as the other methods under all our settings and all sample
sizes. No other method has reliable performance across all settings. First, we
will investigate the five settings in more details. Then we will look more
closely into the comparison between DTRESLO and the other methods. Finally, we
will compare the run-time of different methods.

Figure 4(a) underscores that in the simple setting 1, DTRESLO outperforms all
other methods under both linear and non-linear treatment policies. Figures
4(b) and 4(c) show that in the non-linear settings 2 and 4, as expected, the
non-linear versions of DTRESLO, BOWL, and Q-learning perform better than the
linear counterparts. The only exception is the case of SOWL, which we will
discuss later in more detail. We also observe that setting 4 is quite hard in
that the expected value function of all methods are noticeably lower than the
optimal value function. In settings 2 and 4, non-linear DTRESLO performs
noticeably better than non-linear Q-learning in small sample ($n=$ 250). As
sample size increases, the difference decreases. SOWL has poor performance
under both settings. Although BOWL has better performance than SOWL, its
performance improves rather slowly with $n$ when compared to DTRESLO. This
difference is most noticeable for the non-linear treatment policies under
large samples.

Figure 5(a) implies that under the linear setting 3, value function estimates
of the linear treatment policies are as large as the non-linear policies for
all methods except SOWL. Setting 5, which has larger number of variables, is a
relatively more complicated setting. Although the second stage outcome models
are non-linear in this setting, Figure 5(b) underscores that linear DTRESLO
performs quite comparably to non-linear DTRESLO in large samples under this
setting. Table 5 in Supplement 14 implies that the situation with the other
non-linear DTRESLO methods is similar. Similar to settings 2, in this case,
non-linear DTRESLO has a noticeable edge over all other methods when the
sample size is 250.

Under all settings, non-linear DTRESLO and Q-learning exhibit one particular
pattern, which merits some discussion. Non-linear DTRESLO performs better than
non-linear Q learning in small samples, but their performance becomes almost
similar when the sample size increases to 5000. The relative underperformance
of nonparametric Q-learning in small samples may be due to its heavy reliance
on the correct estimation of Q-functions. Nonparametric estimation of
functions is harder unless the sample size is sufficiently large. In contrast,
DTRESLO only needs to estimate the sign of the blip functions, which is easier
than estimation of the whole function. Finally, this difference may be the
manifestation of the speculated faster regret decay of neural network DTRESLO
(see Section 7). Thus our simulation study complements the theoretical
comparison of the regrets between nonparametric Q-learning and DTRESLO.

DTRESLO outperforms the other direct search method BOWL and SOWL under all
settings except the linear setting, i.e. setting 3. The difference is most
pronounced for non-linear treatment policies in large samples. DTRESLO’s
advantage over BOWL may be attributed to DTRESLO’s simultaneous optimization
approach as opposed to BOWL’s stagewise approach. The latter reduces the
effective sample size in the first stage. In general, SOWL’s average value
function stays quite below the optimal value function. Its performance is
comparable to other methods only in setting 3, where classification is
comparatively easy. What is perhaps more surprising is that the estimated
value function of non-linear SOWL, which is a nonparametric method, does not
even increase with the sample size. The last observation raises the question
whether the approximation error of SOWL decays to zero as sample size
increases. This observation does not refute Theorem 2, and aligns with Remark
5, which speculates that hinge loss as a surrogate may not work for all data
distributions.

Setting | $n$ | DTRESLO | BOWL | SOWL | Q-Learning
---|---|---|---|---|---
Linear | Wavelet | Spline | NeuNet | Linear | RBF | Linear | RBF | Linear | NeuNet
1 | 250 | 0.04 | 0.05 | 0.04 | 0.1 | 1.24 | 21.01 | 0.1 | 0.16 | 0.07 | 0.18
2500 | 0.42 | 0.49 | 0.43 | 0.94 | 13.11 | 655.43 | 54.19 | 80.82 | 0.69 | 2.12
5000 | 0.89 | 1.02 | 0.8 | 2.15 | 77.45 | 3913.85 | 400.32 | 534.36 | 1.4 | 3.94
2 | 250 | 0.04 | 0.06 | 0.04 | 0.16 | 1.36 | 3.48 | 1.3 | 1.33 | 0.09 | 0.19
2500 | 0.54 | 0.6 | 0.42 | 1.01 | 27.75 | 271.08 | 773.79 | 822.42 | 0.73 | 1.83
5000 | 0.88 | 1.25 | 0.91 | 3.7 | 136.88 | 5139.53 | 5901.54 | 5755.75 | 1.49 | 4.15
3 | 250 | 0.05 | 0.05 | 0.04 | 0.15 | 12.59 | 24.39 | 0.08 | 0.13 | 0.08 | 0.2
2500 | 0.41 | 0.71 | 0.42 | 1.03 | 25.75 | 704.16 | 46.55 | 106.67 | 0.73 | 2.04
5000 | 1.22 | 1.04 | 0.84 | 2.19 | 107.99 | 4063.34 | 345.83 | 859.37 | 1.46 | 5.36
4 | 250 | 0.04 | 0.06 | 0.04 | 0.1 | 1.66 | 3.21 | 1.28 | 1.33 | 0.09 | 0.19
2500 | 0.59 | 0.49 | 0.42 | 1.04 | 20.12 | 424.81 | 806.54 | 833.64 | 0.73 | 1.88
5000 | 0.86 | 1.32 | 0.91 | 1.5 | 70.86 | 3317.64 | 5674.96 | 5778.79 | 1.47 | 3.61
5 | 250 | 0.04 | 0.05 | 0.04 | 0.09 | 10.97 | 18.05 | 1.26 | 1.32 | 0.07 | 0.18
2500 | 0.6 | 0.48 | 0.42 | 1.53 | 33.46 | 222.42 | 810.59 | 833.1 | 0.72 | 1.92
5000 | 1.19 | 1.26 | 1.21 | 2.1 | 169.31 | 1012.86 | 5645.25 | 6010.88 | 1.38 | 3.79

Table 3: Run-time for estimating DTR for our smooth surrogates (DTRESLO), Zhao
et al. (2015)’s BOWL & SOWL, and $Q$-learning under settings 1–5.

Table 3 tabulates the run-time of the DTR estimation methods. Run times for
DTRESLO with linear, wavelets, and spline-based treatment policies are
relatively similar. The run-time doubles for neural network treatment
policies. Nonetheless, they are all less than three seconds. Both linear and
neural network $Q$-learning methods are slightly slower than their DTRESLO
counterparts but the difference in run-time is negligible. This is not
surprising because DTRESLO and $Q$-learning methods are trained in a similar
way. They all use stochastic gradient descent with RMSprop for optimization of
the respective loss functions. All these methods are trained for 20 epochs and
use a batch size of 128. As expected, BOWL and SOWL have much larger run-time,
which also increases sharply with $n$. This larger run-time is expected
because SVMs utilize the dual space for optimization. The time cost is
especially high in settings 2 and 4 which have highly non-linear decision
boundaries, and setting 5, which has over 32 features.

To summarize, DTRESLO improves the scalability of existing direct search
methods, achieving run-time as small as Q-learning. We also observe that
within the same class of treatment regimes, i.e. linear or neural network,
DTRESLO outperforms regression-based Q-learning in small sample. This
observation aligns with the existing claim in the literature that
classification is easier than regression in the context of DTR especially in
small samples (Zhao et al., 2015; Kosorok and Laber, 2019). This may happen
because regression-based methods focus on minimizing the $L_{2}(\mathbb{P})$
loss, where the estimation of optimal rules only requires minimization of the
zero-one loss. This mismatch of loss has previously been discussed in
literature (Murphy, 2005; Qian and Murphy, 2011). Our observation thus hints
that bypassing regression may result in better-quality treatment regimes at
least in small samples.

###  8.2 EHR Data Application: DTR for ICU Patients with sepsis

We evaluate our DTRESLO methods and benchmarks on a cohort of $n=9,872$ ICU
patients with sepsis from the Medical Information Mart for Intensive Care
version IV (MIMIC-IV) data (Johnson et al., 2020). Sepsis is a situation when
body’s response to an infection overwhelms the body’s immune system,
potentially causing damage to tissue, multi-organ failure, and in some cases,
death. Individualization of treatment strategies is highly important for
managing sepsis due to the high dissimilarity among sepsis patients (Lat et
al., 2021; Komorowski et al., 2018; Sonabend et al., 2020). Physicians usually
treat it with a high, constant dose of antibiotics. They also use vasopressors
to control blood pressure. Not all patients benefit from vasopressors,
however, and it is also unclear when to administer it (Lat et al., 2021). Here
$d_{t}$ corresponds to the decision regarding whether vasopressors should be
administered at a given time $t$, where we consider $t=1$ at baseline and
$t=2$ at $4$ hour after diagnosis of sepsis. We code our actions as $A_{t}=-1$
if no dose is necessary and $A_{t}=1$ otherwise. The state space is comprised
of 46 covariates $O_{t}$ at each time-step. Measured variables include age,
body mass index, diastolic and systolic blood pressure, etc. We use an inverse
transformation of lactate acid level as the outcome $Y_{t}$. In particular
$Y_{t}=(LA_{t}+5)^{-1}+2$ where $LA_{t}$ stands for lactic acid level at time
$t$. Adding the offset $2$ ensures that $Y_{t}>0$, which is required by
Assumption A.

We evaluate DTRESLO using linear, wavelets, splines and neural network
functions. The comparators are as in Section 8.1. For DTRESLO, we choose from
four choices of $\phi$ presented in Example 1. To this end, we use a five fold
cross validation. To elaborate, for each split, we estimate the DTRs using 4
folds of the sample. Then we estimate the value function of each treatment
regime using the remaining fold. Thus DTRESLO with above-mentioned four
choices of $\phi$ yield four potentially different treatment regimes in the
first step. However, we pick only one regime among these four regimes. We
consider that treatment regime, which has the highest value function estimate.
We also report the mean and SE of the value function estimates for both
DTRESLO and the comparators.

We estimate the value functions of the derived DTRs using a doubly robust
approach (Jiang and Li, 2016; Thomas and Brunskill, 2016; Sonabend et al.,
2021). The corresponding value function estimator is given by

$\displaystyle\begin{split}\widehat{V}_{DR}(\widehat{d}_{1},\widehat{d}_{2})=&\mathbb{P}_{n}\bigg{[}\widehat{Q}_{1}(H_{1},\hat{d}_{1}(H_{1}))\\\
+&\frac{1[A_{1}=\hat{d}_{1}(H_{1})]}{\widehat{\pi}_{1}(A_{1}|H_{1})}\left[Y_{1}-\left\\{\widehat{Q}_{1}(H_{1},\hat{d}_{1}(H_{1}))-\widehat{Q}_{2}(H_{2},\hat{d}_{2}(H_{2}))\right\\}\right]\\\
+&\frac{1[A_{1}=\hat{d}_{1}(H_{1}),A_{2}=\hat{d}_{2}(H_{2})]}{\widehat{\pi}_{1}(A_{1}|H_{1})\widehat{\pi}_{2}(A_{2}|H_{2})}\left\\{Y_{2}-\widehat{Q}_{2}(H_{2},\hat{d}_{2}(H_{2}))\right\\}\bigg{]},\end{split}$

where $\widehat{\pi}_{1}$ and $\widehat{\pi}_{2}$ are the propensity score
estimators, and $\widehat{Q}_{1}$ and $\widehat{Q}_{2}$ are estimators of the
Q-functions corresponding to the first and second stage, respectively. Here
the $Q$ functions are estimated using a neural network with the same
specifications as our nonparametric Q-learning DTR estimator.

|  | DTRESLO
---
(linear)
| DTRESLO
---
(splines)
| DTRESLO
---
(wavelets)
| DTRESLO
---
(NN)
| BOWL
---
(linear)
| BOWL
---
(RBF)
| SOWL
---
(linear)
| SOWL
---
(RBF)
| Q-learning
---
(linear)
| Q-learning
---
(NN)
$\hat{V}_{DR}(\hat{d}_{1},\hat{d}_{2})$ | 1.352 | 1.415 | 1.467 | 1.458 | 1.248 | 1.450 | 0.845 | 1.137 | 1.062 | 1.144
SD | 0.043 | 0.084 | 0.024 | 0.023 | 0.039 | 0.071 | 0.148 | 0.082 | 0.095 | 0.102

Table 4: Estimated Value function and Standard deviation for the DTRs derived
with different methods. Here NN stands for neural network.

Table 4 reports the estimated values using the doubly robust estimator for
each method and their corresponding standard errors. We observe that DTRESLO
with wavelet treatment policies has the highest value function estimate. Also,
the estimated standard deviations of non-linear (nonparametric) DTRESLO
estimators are smaller compared to the nonparametric counterparts of BOWL,
SOWL, and Q-learning. It must be noted that the linear DTRESLO outperforms
linear $Q$-learning, which was observed in all our simulation settings as
well. Also, as in our simulations, SOWL has lower value function estimate
compared to all other methods. This aligns with our observations with SOWL in
Section 8.1.

##  9 Discussion

Our work _is the first step_ towards a unified understanding of general
surrogate losses in the simultaneous optimization context. Our work leaves
ample room for modification and generalization to complex real-world
scenarios. We list below some important open questions:

The theoretical results in this paper consider the propensity scores to be
known. They may be available in SMART studies, but they need to be estimated
for observational studies. At best, we may be able to estimate the propensity
scores at $n^{-1/2}$-rate. Therefore it is possible that in this situation,
our regret-decay rate will slow down. We also do not know if it is at all
possible to push the regret decay to $O_{p}(1/n)$ in this situation because we
do not know the minimax rate of regret-decay in this context. However, there
is a more pressing issue with the use of inverse propensity score weighting.
The weight will grow smaller as the number of stages increase, leading to a
highly volatile method (Kosorok and Laber, 2019). However, there are
strategies (Kallus, 2018) that can be incorporated to ensure robustness.
Research in this direction is needed to increase the stability of our DTRESLO
method.

Also, there are many choices of $\phi$’s that satisfy Condition 2, and hence
can be used for DTRESLO. In this paper, we have not considered the problem of
selecting a $\phi$. We fixed a particular $\phi$ in our empirical study but
the performance may be improved by a more careful tuning of $\phi$.

The DTRESLO method easily extends to $K>2$ by using a surrogate
$\psi(x_{1},\ldots,x_{K})=\phi(x_{1})\ldots\phi(x_{K})$. Although we do not
yet know whether Fisher consistency still holds, our proof techniques are
readily extendable to the higher stages via mathematical induction. If our
DTRESLO method is Fisher consistent for general $K$ stages, the pattern of
error accumulation over stages will be an immediate interest. For Q-learning,
the regret grows exponentially with the number of stages (Murphy, 2005). In
view of Wang et al. (2020), exponential error accumulation may sometimes be
inevitable under very general conditions. However, we wonder whether our
simultaneous maximization procedure escapes the exponential error accumulation
in presence of noise conditions.

Despite being of immense practical interest, this area greatly lacks direct
search method with rigorous guarantees at multi-stage setting. Direct search
method with more than two levels of treatment requires integration of
multicategory classification with the sequential setting of DTR, and hence is
conceptually more challenging than the regression-based counterparts. However,
we expect that DTRESLO can be extended to identify optimal DTRs under this
more complex setting. Detailed strategies for identifying the surrogate loss
and implementing algorithms to estimate DTRs in practice warrant future
research.

##  10 Acknowledgments

Rajarshi Mukherjee and Nilanjana Laha’s research was supported by National
Institutes of Health grant P42ES030990. Tianxi Cai’s research was supported by
National Institutes of Health grant R01LM013614. Aaron Sonabend’s research was
supported by the Boehringer-Ingelheim Fellowship at Harvard.

##  11 Supplement

The supplement contains the following:

  1. 1.

Further details on some assumptions made in this paper, further details on the
regret decomposition in Section 4.1, and further details on Table 2.

  2. 2.

More details on the simulation settings in Section 8.1.

  3. 3.

Proof of all theorems, lemmas, and corollaries.

## References

  * Andersen et al. (2021) Andersen, M. S., Dahl, J., and Vandenberghe., L. (2021).  Cvxopt: A python package for convex optimization, version 1.2.  Available at https://cvxopt.org, 2021.
  * Audibert et al. (2007) Audibert, J.-Y., Tsybakov, A. B., et al. (2007).  Fast learning rates for plug-in classifiers.  The Annals of statistics, 35(2), 608–633.
  * Bartlett and Mendelson (2002) Bartlett, P. L. and Mendelson, S. (2002).  Rademacher and gaussian complexities: Risk bounds and structural results.  Journal of Machine Learning Research, 3(Nov), 463–482.
  * Bartlett et al. (2006) Bartlett, P. L., Jordan, M. I., and McAuliffe, J. D. (2006).  Convexity, classification, and risk bounds.  Journal of the American Statistical Association, 101(473), 138–156.
  * Blanchard et al. (2008) Blanchard, G., Bousquet, O., Massart, P., et al. (2008).  Statistical performance of support vector machines.  The Annals of Statistics, 36(2), 489–531.
  * Calauzenes et al. (2012) Calauzenes, C., Usunier, N., and Gallinari, P. (2012).  On the (non-) existence of convex, calibrated surrogate losses for ranking.  Advances in Neural Information Processing Systems 25 (NIPS 2012), pages 197–205.
  * Chen et al. (2016) Chen, G., Zeng, D., and Kosorok, M. R. (2016).  Personalized dose finding using outcome weighted learning.  Journal of the American Statistical Association, 111(516), 1509–1521.
  * Chen et al. (2017) Chen, S., Tian, L., Cai, T., and Yu, M. (2017).  A general statistical framework for subgroup identification and comparative treatment scoring.  Biometrics, 73(4), 1199–1209.
  * Cui and Tchetgen Tchetgen (2020) Cui, Y. and Tchetgen Tchetgen, E. (2020).  A semiparametric instrumental variable approach to optimal treatment regimes under endogeneity.  Journal of the American Statistical Association, 116(533), 162–173.
  * Dembczyński et al. (2012) Dembczyński, K., Waegeman, W., Cheng, W., and Hüllermeier, E. (2012).  On label dependence and loss minimization in multi-label classification.  Machine Learning, 88(1-2), 5–45.
  * Duchi et al. (2018) Duchi, J., Khosravi, K., Ruan, F., et al. (2018).  Multiclass classification, information, divergence and surrogate risk.  The Annals of Statistics, 46(6B), 3246–3275.
  * Duchi et al. (2010) Duchi, J. C., Mackey, L. W., and Jordan, M. I. (2010).  On the consistency of ranking algorithms.  In ICML.
  * Eberts and Steinwart (2011) Eberts, M. and Steinwart, I. (2011).  Optimal learning rates for least squares svms using gaussian kernels.  Advances in neural information processing systems, 24, 1539–1547.
  * Gao and Zhou (2011) Gao, W. and Zhou, Z.-H. (2011).  On the consistency of multi-label learning.  In Proceedings of the 24th annual conference on learning theory, pages 341–358.
  * Giné and Zinn (1984) Giné, E. and Zinn, J. (1984).  Some limit theorems for empirical processes.  The Annals of Probability, pages 929–989.
  * Giné and Nickl (2015) Giné, E. and Nickl, R. (2015).  Mathematical Foundations of Infinite-Dimensional Statistical Models, volume 40 of Cambridge series in statistical and probabilistic mathematics.  Cambridge University Press, Cambridge.
  * Hinton et al. (2012) Hinton, G., Srivastava, N., and Swersky, K. (2012).  Neural networks for machine learning, overview of mini-batch gradient descent.  Available at https://www.cs.toronto.edu/~tijmen/csc321/slides/lecture_slides_lec6.pdf.
  * Hiriart-Urruty and Lemaréchal (2004) Hiriart-Urruty, J.-B. and Lemaréchal, C. (2004).  Fundamentals of convex analysis.  Springer Science & Business Media.
  * Jiang et al. (2019) Jiang, B., Song, R., Li, J., and Zeng, D. (2019).  Entropy learning for dynamic treatment regimes.  Statistica Sinica, 29(4), 1633–1655.
  * Jiang and Li (2016) Jiang, N. and Li, L. (2016).  Doubly robust off-policy value evaluation for reinforcement learning.  arXiv.org.
  * Johnson et al. (2020) Johnson, A., Bulgarelli, L., Pollard, T., Horng, S., Celi, L. A., and Mark., R. (2020).  MIMIC-IV (version 0.4). PhysioNet. .  https://doi.org/10.13026/a3wn-hq05.
  * Kallus (2018) Kallus, N. (2018).  Balanced policy evaluation and learning.  Advances in neural information processing systems, 31.
  * Kallus (2020) Kallus, N. (2020).  Comment: Entropy learning for dynamic treatment regimes.  arXiv preprint arXiv:2004.02778.
  * Karp (1972) Karp, R. M. (1972).  Reducibility among combinatorial problems.  In Complexity of computer computations, pages 85–103. Springer.
  * Koltchinskii (2011) Koltchinskii, V. (2011).  Oracle inequalities in empirical risk minimization and sparse recovery problems: École d’Été de Probabilités de Saint-Flour XXXVIII-2008, volume 2033.  Springer Science & Business Media.
  * Komorowski et al. (2018) Komorowski, M., Celi, L. A., Badawi, O., Gordon, A. C., and Faisal, A. A. (2018).  The artificial intelligence clinician learns optimal treatment strategies for sepsis in intensive care.  Nature medicine, 24(11).
  * Kosorok and Laber (2019) Kosorok, M. R. and Laber, E. B. (2019).  Precision medicine.  Annual review of statistics and its application, 6, 263–286.
  * Laber and Davidian (2017) Laber, E. B. and Davidian, M. (2017).  Dynamic treatment regimes, past, present, and future: A conversation with experts.  Statistical methods in medical research, 26(4), 1605–1610.
  * Laber and Zhao (2015) Laber, E. B. and Zhao, Y.-Q. (2015).  Tree-based methods for individualized treatment regimes.  Biometrika, 102(3), 501–514.
  * Laber et al. (2014) Laber, E. B., Lizotte, D. J., Qian, M., Pelham, W. E., and Murphy, S. A. (2014).  Dynamic treatment regimes: Technical challenges and applications.  Electronic journal of statistics, 8(1), 1225.
  * Lat et al. (2021) Lat, I., Coopersmith, C. M., and De Backer, D. (2021).  The surviving sepsis campaign: fluid resuscitation and vasopressor therapy research priorities in adult patients.  Intensive Care Medicine Experimental, 9(1), 1–16.
  * Lederer (2020) Lederer, J. (2020).  Risk bounds for robust deep learning.  arXiv preprint arXiv:2009.06202.
  * Lin (2004) Lin, Y. (2004).  A note on margin-based loss functions in classification.  Statistics & probability letters, 68(1), 73–82.
  * Liu (2007) Liu, Y. (2007).  Fisher consistency of multicategory support vector machines.  In Artificial intelligence and statistics, pages 291–298.
  * Luedtke and Van Der Laan (2016) Luedtke, A. R. and Van Der Laan, M. J. (2016).  Statistical inference for the mean outcome under a possibly non-unique optimal treatment strategy.  Annals of statistics, 44(2), 713.
  * Mammen et al. (1999) Mammen, E., Tsybakov, A. B., et al. (1999).  Smooth discrimination analysis.  Annals of Statistics, 27(6), 1808–1829.
  * Massart et al. (2006) Massart, P., Nédélec, É., et al. (2006).  Risk bounds for statistical learning.  The Annals of Statistics, 34(5), 2326–2366.
  * Moodie et al. (2014) Moodie, E. E., Dean, N., and Sun, Y. R. (2014).  Q-learning: Flexible learning about useful utilities.  Statistics in Biosciences, 6(2), 223–243.
  * Murphy (2005) Murphy, S. (2005).  A generalization error for q-learning.  Journal Of Machine Learning Research, 6, 1073–1097.
  * Murphy (2003) Murphy, S. A. (2003).  Optimal dynamic treatment regimes.  J. R. Stat. Soc. Series. B Stat. Methodol., 65, 331–355.
  * Murphy et al. (2001) Murphy, S. A., van der Laan, M. J., and Robins, J. M. (2001).  Marginal mean models for dynamic regimes.  Journal of the American Statistical Association, 96(456), 1410–1423.
  * Neykov et al. (2016) Neykov, M., Liu, J. S., and Cai, T. (2016).  On the characterization of a class of fisher-consistent loss functions and its application to boosting.  The Journal of Machine Learning Research, 17(1), 2498–2529.
  * Orellana et al. (2010) Orellana, L., Rotnitzky, A., and Robins, J. 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(2).
  * Pedregosa et al. (2017) Pedregosa, F., Bach, F., and Gramfort, A. (2017).  On the consistency of ordinal regression methods.  Journal of Machine Learning Research, 18, 1–35.
  * Qian and Murphy (2011) Qian, M. and Murphy, S. A. (2011).  Performance guarantees for individualized treatment rules.  Annals of statistics, 39(2), 1180.
  * Robins (2004) Robins, J. M. . (2004).  Optimal structural nested models for optimal sequential decisions.  Proceedings of the second Seattle symposium on biostatistics In D. Y. Lin and P. Heagerty (Eds.) (pp. 189–326). New York: Springer.
  * Robins (1994) Robins, J. M. (1994).  Correcting for non-compliance in randomized trials using structural nested mean models.  Communications in Statistics-Theory and methods, 23(8), 2379–2412.
  * Robins (1997) Robins, J. M. (1997).  Causal inference from complex longitudinal data.  In Latent variable modeling and applications to causality, pages 69–117. Springer.
  * Schmidt-Hieber (2020) Schmidt-Hieber (2020).  Nonparametric regression using deep neural networks with relu activation function.  Annals of Statistics, 48(4), 1875–1897.
  * Schulte et al. (2014) Schulte, P. J., Tsiatis, A. A., Laber, E. B., and Davidian, M. (2014).  $\mathbf{Q}$\- and $\mathbf{A}$-learning methods for estimating optimal dynamic treatment regimes.  Statist. Sci., 29(4), 640–661.
  * Sonabend et al. (2020) Sonabend, A., Lu, J., Celi, L. A., Cai, T., and Szolovits, P. (2020).  Expert-supervised reinforcement learning for offline policy learning and evaluation.  In Advances in Neural Information Processing Systems, volume 33, pages 18967–18977.
  * Sonabend et al. (2021) Sonabend, A., Laha, N., Ananthakrishnan, A. N., Cai, T., and Mukherjee, R. (2021).  Semi-supervised off policy reinforcement learning.
  * Song et al. (2015) Song, R., Kosorok, M., Zeng, D., Zhao, Y., Laber, E., and Yuan, M. (2015).  On sparse representation for optimal individualized treatment selection with penalized outcome weighted learning.  Stat, 4(1), 59–68.
  * Steinwart et al. (2007) Steinwart, I., Scovel, C., et al. (2007).  Fast rates for support vector machines using gaussian kernels.  The Annals of Statistics, 35(2), 575–607.
  * Sun and Wang (2021) Sun, Y. and Wang, L. (2021).  Stochastic tree search for estimating optimal dynamic treatment regimes.  Journal of the American Statistical Association, 116(533), 421–432.
  * Tewari and Bartlett (2007) Tewari, A. and Bartlett, P. L. (2007).  On the consistency of multiclass classification methods.  Journal of Machine Learning Research, 8(May), 1007–1025.
  * Thomas and Brunskill (2016) Thomas, P. and Brunskill, E. (2016).  Data-efficient off-policy policy evaluation for reinforcement learning.  In International Conference on Machine Learning, pages 2139–2148. PMLR.
  * Tsybakov et al. (2004) Tsybakov, A. B. et al. (2004).  Optimal aggregation of classifiers in statistical learning.  The Annals of Statistics, 32(1), 135–166.
  * Van Der Vaart et al. (1996) Van Der Vaart, A. W., van der Vaart, A. W., van der Vaart, A., and Wellner, J. (1996).  Weak convergence and empirical processes: with applications to statistics.  Springer Science & Business Media.
  * Wang et al. (2020) Wang, R., Foster, D. P., and Kakade, S. M. (2020).  What are the statistical limits of offline RL with linear function approximation?  arXiv preprint arXiv:2010.11895.
  * Watkins (1989) Watkins, C. J. C. H. (1989).  Learning from delayed rewards.
  * Wellner (2005) Wellner, J. A. (2005).  Empirical processes: Theory and applications.  Lecture notes.
  * Williams and Seeger (2001) Williams, C. and Seeger, M. (2001).  Using the nyström method to speed up kernel machines.  In T. Leen, T. Dietterich, and V. Tresp, editors, Advances in Neural Information Processing Systems, volume 13. MIT Press.
  * Xu et al. (2016) Xu, Y., Müller, P., Wahed, A. S., and Thall, P. F. (2016).  Bayesian nonparametric estimation for dynamic treatment regimes with sequential transition times.  Journal of the American Statistical Association, 111(515), 921–950.
  * Yang (1999) Yang, Y. (1999).  Minimax nonparametric classification. i. rates of convergence.  IEEE Transactions on Information Theory, 45(7), 2271–2284.
  * Zajonc (2012) Zajonc, T. (2012).  Bayesian inference for dynamic treatment regimes: Mobility, equity, and efficiency in student tracking.  Journal of the American Statistical Association, 107(497), 80–92.
  * Zhang et al. (2018a) Zhang, J., Liu, T., and Tao, D. (2018a).  On the rates of convergence from surrogate risk minimizers to the bayes optimal classifier.  arXiv preprint arXiv:1802.03688.
  * Zhang (2010) Zhang, T. (2010).  Analysis of multi-stage convex relaxation for sparse regularization.  Journal of Machine Learning Research, 11(3).
  * Zhang et al. (2018b) Zhang, Y., Laber, E. B., Davidian, M., and Tsiatis, A. A. (2018b).  Interpretable dynamic treatment regimes.  Journal of the American Statistical Association, 113(524), 1541–1549.
  * Zhao et al. (2012) Zhao, Y., Zeng, D., Rush, A. J., and Kosorok, M. R. (2012).  Estimating individualized treatment rules using outcome weighted learning.  Journal of the American Statistical Association, 107(499), 1106–1118.
  * Zhao et al. (2015) Zhao, Y.-Q., Zeng, D., Laber, E. B., and Kosorok, M. R. (2015).  New statistical learning methods for estimating optimal dynamic treatment regimes.  Journal of the American Statistical Association, 110, 583–598.
  * Zhou and Kosorok (2017) Zhou, X. and Kosorok, M. R. (2017).  Augmented outcome-weighted learning for optimal treatment regimes.  arXiv preprint arXiv:1711.10654.

Supplement to “Finding the Optimal Dynamic Treatment Regime Using Smooth
Fisher Consistent Surrogate Loss”

The supplement is organized as follows. In Supplement 12, we provide further
details on some assumptions made in our paper and the regret decomposition in
Section 4.1. Supplement 13 contains more details on table 2. Supplement 14
provides more details on the simulation settings in Section 8.1. Rest of the
supplement collects the proofs.

##  12 Additional discussion on assumptions and regret decomposition

###  12.1 On the assumptions in Theorem 1

Closure: Closure ensures that the concave function is well behaved at the
boundary of its domain. In particular, if $h(x)>-\infty$ for all
$x\in\mathbb{R}^{k}$, then $h$ is automatically closed. Thus all loss
functions considered in this paper are closed. However, if $\text{dom}(\psi)$
is a strict subset of $\mathbb{R}^{2}$, then whether $\psi$ is closed or not
depends on the value of $\psi$ at the boundary of $\text{ri}(\text{dom}(h))$.
Precisely, $h$ is closed if for all $x\in\mathbb{R}^{k}$,
$h(x)=\lim_{t\downarrow 0}h(x+t(x^{\prime}-x))$ for some
$x^{\prime}\in\text{ri}(\text{dom}(h))$ (Hiriart-Urruty and Lemaréchal, 2004,
Proposition 1.2.5). we assume $\psi$ to be closed in order to avoid the