Policy Learning for Optimal Individualized Dose Intervals

Guanhua Chen Xiaomao Li Menggang Yu

University of Wisconsin, Madison<EMAIL_ADDRESS>University of Wisconsin,
Madison<EMAIL_ADDRESS>University of Wisconsin, Madison
<EMAIL_ADDRESS>
###### Abstract

We study the problem of learning individualized dose intervals using
observational data. There are very few previous works for policy learning with
continuous treatment, and all of them focused on recommending an optimal dose
rather than an optimal dose interval. In this paper, we propose a new method
to estimate such an optimal dose interval, named probability dose interval
(PDI). The potential outcomes for doses in the PDI are guaranteed better than
a pre-specified threshold with a given probability (e.g., $50\%$). The
associated nonconvex optimization problem can be efficiently solved by the
Difference-of-Convex functions (DC) algorithm. We prove that our estimated
policy is consistent, and its risk converges to that of the best-in-class
policy at a root-n rate. Numerical simulations show the advantage of the
proposed method over outcome modeling based benchmarks. We further demonstrate
the performance of our method in determining individualized Hemoglobin A1c
(HbA1c) control intervals for elderly patients with diabetes.

##  1 INTRODUCTION

Increasing evidence has shown that personalized prescriptions based on
patients’ health conditions and medical history will result in more effective
treatment (Mancinelli et al.,, 2000; Collins and Varmus,, 2015). Hence, it is
desirable to learn personalized treatment policy based on collected data and
recommend treatment for future patients based on such a policy. We focus on
off-policy learning, which is about evaluating and optimizing new personalized
treatment policies based on observational data.

Despite the abundance of research on policy learning for categorical
treatments (Liang and Yu,, 2020; Lou et al.,, 2018; Chen et al.,, 2017;
Kosorok and Laber,, 2019), policy learning for continuous treatments has not
been studied sufficiently. Meanwhile, many applications in the medical domain
involve determining optimal continuous treatment: the optimal dosage of a
medical drug, e.g., Warfarin (International Warfarin Pharmacogenetics
Consortium,, 2009; Rich et al.,, 2014), and the optimal length of stay in the
hospital for hospitalized patients (Borghans et al.,, 2012). Existing
approaches for continuous treatments focused on recommending a fixed optimal
dose for each patient. However, for many medical applications, recommending
optimal dose interval can be more practical. For example, for chronic disease
management metrics (e.g. HbA1c level, blood pressure, and cholesterol),
targeted intervals, instead of a particular level, are frequently recommended
from major medical organizations such as the American Diabetes Association
(American Association of Clinical Endocrinologists,, 2019) and the American
Heart Association (Flack and Adekola,, 2020). Another example is cancer
radiation therapy, where radiation received by targeted tumors, not delivered
by physicians, can be controlled only at an interval due to the irregular
shape of solid tumors and the need to radiate surrounding regions (Scott et
al.,, 2017). Therefore we focus on policy learning for optimal recommendation
of an interval instead of a single dose. In addition, this shifted focus will
help relax the local approximation needed for optimal dose finding, therefore
improving both computational stability and theoretical properties, including a
faster convergence rate for the target function over an existing method (Chen
et al.,, 2016).

##  2 PROBLEM SETUP

###  2.1 Settings and Assumptions

Denote the dataset with $n$ observations by $(X_{i},A_{i},Y_{i})$,
$i=1,\ldots,n$. Here, $X_{i}\in\mathcal{X}$ is a $d$-dimensional vector
representing patient-level covariates.
$A_{i}\in\mathcal{A}\equiv[a_{L},a_{U}]$ is the observed treatment dose whose
range has the lower and upper bounds ($a_{L}$ and $a_{U}$). $Y_{i}$ is the
observed outcome from receiving $A_{i}$. We adopt the well-known potential
outcome notation in causal inference (Imbens and Rubin,, 2015) to articulate
our research problem. For any dose $a\in\mathcal{A}$, the potential outcome
$Y_{i}(a)$ is the outcome that would have been observed if patient $i$
received treatment $a$. Note that when $A_{i}=a$, naturally $Y_{i}=Y_{i}(a)$
by the causal consistency assumption (Imbens and Rubin,, 2015). An
individualized policy or dose assignment rule is a functional
$f(X):\mathcal{X}\rightarrow\mathcal{A}$. For simplicity of presentation, we
assume bigger $Y(a)$ is better. The optimal policy is the one that maximizes
the policy value $V(f)\equiv E[Y(f(X))]$. To estimate this causal quantity
using the observed data, we adopt the following three commonly assumed
conditions (Imbens and Rubin,, 2015): 1) the stable unit treatment value
assumption (SUTVA); 2) the no unmeasured confounding assumption: $Y(a)\perp
A|X$; and 3) the positivity assumption which assumes $P(A=a|X)>0$ for all
$a\in\mathcal{A}$. $P(A=a|X)$ is known as the (generalized) propensity score
function that models the treatment assignment mechanism (Imbens and Rubin,,
2015).

An Individualized Dose Interval (IDI) is desirable for a patient if the doses
inside the interval are more likely associated with better outcomes than doses
outside the interval. Each IDI is therefore connected with a threshold $S$
that separates good outcomes from bad ones for the patient. In addition, a
confidence probability $\alpha$ reflects how strongly the desirable interval
is associated with the good outcomes. Specifically, let $f_{L}(x)$ and
$f_{U}(x)$ be the lower and upper bounds of the IDI for a patient with
covariates $x$ such that $a_{L}\leq f_{L}(x)\leq f_{U}(x)\leq a_{U},\;\forall
x\in\mathcal{X}$. Then, we define a level-$\alpha$ Probability Dose Interval
(PDI) as follows:

$\displaystyle\text{PDI}_{\alpha}(x)\equiv\big{[}f_{L}(x),f_{U}(x)\big{]}:\
\forall a\in[f_{L}(x),f_{U}(x)],$ (1) $\displaystyle
P\\{Y(a)>S|x\\}\geq\alpha.$

If a clinician follows the PDI recommendation when treating a patient with
covariates $x$, the patient is guaranteed to have the outcome above $S$ with a
probability of at least $\alpha$.

###  2.2 Related work

Indirect methods: Indirect methods for estimating treatment rules are based on
outcome modeling (Zhao et al.,, 2009; Moodie et al.,, 2012; Chakraborty and
Moodie,, 2013). This framework has two steps. The first step models the
conditional expectation of the outcome given the treatment and covariates,
$E(Y|A,X)$. The models can be parametric such as regression with Lasso-type
regularization (Tibshirani,, 1996), or nonparametric such as regression tree
(Breiman et al.,, 1984; Breiman,, 2001) and Support Vector Machine (Steinwart
and Christmann,, 2008; Vapnik,, 2013). In the second step, the treatment
option which maximizes the conditional expection in the first step is set as
the optimal treatment. For most cases, the second step involves numerical
optimization such as grid search due to lack of close-form solutions. Since
there are no existing approaches for PDI learning, we propose the following
indirect PDI learning approach to serve as a benchmark. In particular, the
first step for indirect PDI learning would construct an outcome model to
estimate the conditional probability $P\big{(}Y>S\,|\,A=a,X=x\big{)}$, for
example, based on a classifier for predicting the binary label $I\\{Y>S\\}$
using $A$ and $X$. Subsequently, a grid search for doses is done to identify
the dose interval where the predicted conditional probability is greater than
the confidence level $\alpha$. It has been noted in binary treatment settings
that solutions from the indirect methods can be unstable due to the need to
model $P\big{(}Y>S\,|\,A=a,X=x\big{)}$ for all $a$ and $x$ (Kosorok and
Laber,, 2019). Furthermore, the grid search step could be computationally
challenging.

Direct methods: On the other hand, for optimal policy learning with
categorical treatments, direct methods have been developed to learn the
optimal policy directly, without using the outcome modeling. These approaches
(Moodie et al.,, 2009; Xu et al.,, 2015; Lou et al.,, 2018; Chen et al.,,
2017; Liang and Yu,, 2020) are believed to have better performance when
signal-to-noise ratio is low. For learning optimal policy for continuous
treatments, Chen et al., (2016) proposed an outcome weighted learning
(O-learning) method which is based on a weighted regression. In particular,
O-learning maximizes the value function associated with a policy $f$: i.e.
$f_{opt}=\text{argmax}_{f}V(f)$ with $V(f)=E[Y(f(X))]$. Furthermore, Chen et
al., (2016) and other direct methods (Kallus and Zhou,, 2018; Bica et al.,,
2020; Sondhi et al.,, 2020; Zhu et al.,, 2020; Schulz and Moodie,, 2020; Zhou
et al.,, 2021) use inverse probability weighted (IPW) or augmented IPW
estimators of the value function to address confounding issues by reweighting
the observational data to mimic randomized trials. However, none of these
approaches are directly applicable for the PDI problem.

##  3 DIRECT POLICY LEARNING FOR PDI

Because the recommended dose interval needs to satisfy probability constraint
$P\big{(}Y(a)>S\,|\,X=x\big{)}\geq\alpha$ for $a\in\text{PDI}_{\alpha}(x)$,
maximizing the associated value function for a dose interval becomes less
intuitive. Instead, we propose a loss function such that the estimated optimal
policy that minimizes the loss function would (asymptotically) satisfy the
constraint.

There are two types of errors that an PDI can incur. When a patient receives a
dose inside the PDI, $a\in\text{PDI}_{\alpha}(x)$, but has an outcome below
the threshold, $Y(a)<S$. It presents a false-positive error. On the other
hand, when the patient receives a dose outside the interval,
$a\not\in\text{PDI}_{\alpha}(x)$, but has an outcome above the threshold,
$Y(a)>S$. It presents a false-negative error. The goal is to minimize these
two errors. Instead of working with a multi-objective optimization problem, we
use $\alpha$ to balance the magnitudes of the two errors and work with the
weighted sum, that is, $\alpha I\\{Y(a)\leq
S,a\in\text{PDI}_{\alpha}(x)\\}+(1-\alpha)I\\{Y(a)>S,a\not\in\text{PDI}_{\alpha}(x)\\}$.
This weighted sum leads to nice theoretical guarantee of the resulting
interval for our proposed method as we demonstrate in Section 4 below.

The indirect methods rely on modeling the outcome $Y$ to predict $I(Y(a)>S)$.
Instead, we circumvent modeling of the outcome and directly minimize the above
error based on observed data. Specifically, we consider the following risk
function for minimization to obtain an estimate for $\text{PDI}_{\alpha}(x)$
with any given $\alpha$, and $S$.

$\displaystyle E\bigg{[}\Big{\\{}\alpha I(Y\leq
S)I(A\in\text{PDI}_{\alpha}(X))$ (2)
$\displaystyle+(1-\alpha)I(Y>S)I(A\notin\text{PDI}_{\alpha}(X))\Big{\\}}\frac{1}{P(A\,|\,X)}\bigg{]}$

Here $P(A\,|\,X)$ is the generalized propensity score account for possible
nonrandom treatment assignment in observational studies.

Under mild conditions, the minimizer of the risk function (LABEL:risk-PDI)
becomes an interval $[f_{L}(x),f_{U}(x)]$ where $f_{L}(x)$ and $f_{U}(x)$ are
boundary functions that minimize the following risk

$\displaystyle R(f_{L},f_{U})=E\bigg{[}\Big{\\{}\alpha I(Y\leq
S)I(A\in[f_{L}(X),f_{U}(X)])$ (3)
$\displaystyle+(1-\alpha)I(Y>S)I(A\notin[f_{L}(X),f_{U}(X)])\Big{\\}}\frac{1}{P(A\,|\,X)}\bigg{]}$

In cases when Conditions C1 and C2 defined in Section 4 are satisfied, then
only the lower bound is needed to separate the desirable doses from the
undesirable doses, i.e., $\text{PDI}_{\alpha}(x)=[f_{L}(X),a_{U}]$.
Accordingly, we just need to minimize the following to estimate $f_{L}(x)$.

$\displaystyle R(f_{L})=E\bigg{[}\Big{\\{}\alpha I(Y\leq
S)I\big{(}A\in[f_{L}(X),a_{U}]\big{)}$ (4)
$\displaystyle+(1-\alpha)I(Y>S)I\big{(}A\in[a_{L},f_{L}(X)]\big{)}\Big{\\}}\frac{1}{P(A\,|\,X)}\bigg{]}$

The functional forms for the boundary functions $f_{L}(x)$ and $f_{U}(x)$ can
be chosen quite flexibly in principle, as in many machine learning approaches.
In practice, one may put some some smoothness conditions on them for given $x$
to ensure that a small perturbation of $x$ will not lead to a dramatic change
in the resulting PDI.

###  3.1 Non-Convex Surrogate Loss Relaxation

As in the classical classification tasks, it is difficult to directly optimize
the empirical risk due to the discontinuity of indicator functions in (3) and
(4). Surrogate loss functions have been proposed to replace the indicator
functions as long as the surrogate loss functions satisfy certain conditions
such as Fisher consistency. Note that unbounded loss functions such as the
hinge loss are not good options here as demonstrated in the optimal dose
problem (Chen et al.,, 2016). The reason is due to the sensitivity of the
solution of the unbounded surrogate loss to outliers, that is, to observations
whose received doses are far from the optimal bounds. Therefore we propose a
truncated hinge loss as a surrogate similar to Chen et al., (2016). Such non-
convex surrogate losses have also been used in robust SVM (Wu and Liu,, 2007).

The risk function (4) of the one-sided $\text{PDI}_{\alpha}(X)$ after
relaxation is as follows. See Figure 1 for a graphical representation.

$\displaystyle R_{\epsilon}(f_{L})=E\bigg{[}\Big{\\{}\alpha I\big{(}Y\leq
S\big{)}\Psi_{\epsilon}(A,{f}_{L})$ (5)
$\displaystyle+\big{(}1-\alpha\big{)}I\big{(}Y>S\big{)}\Psi_{\epsilon}({f}_{L},A)\Big{\\}}\frac{1}{P(A\,|\,X)}\bigg{]}$

where $\Psi_{\epsilon}(a,b)=min\\{{\epsilon}^{-1}{(a-b)_{+}},1\\}.$ Note that
$\Psi_{\epsilon}(a,b)=\Psi_{\epsilon,1}(a,b)-\Psi_{\epsilon,2}(a,b)$, where
$\Psi_{\epsilon,1}(a,b)=\\{{\epsilon}^{-1}{(a-b)}\\}_{+}$ and
$\Psi_{\epsilon,2}(a,b)=\\{{\epsilon}^{-1}{(a-b)}-1\\}_{+}$ are convex.

The risk function (3) of the two-sided $\text{PDI}_{\alpha}(X)$ after
relaxation is as follows.

$\displaystyle R_{\epsilon}(f_{L},f_{U})=E\bigg{[}\Big{\\{}\alpha
I\big{(}Y\leq S\big{)}\Psi_{\epsilon}^{in}\big{(}{f}_{L},A,f_{U}\big{)}$ (6)
$\displaystyle+\big{(}1-\alpha\big{)}I\big{(}Y>S\big{)}\Psi_{\epsilon}^{out}\big{(}{f}_{L},A,f_{U}\big{)}\Big{\\}}\frac{1}{P(A\,|\,X)}\bigg{]}$

where

$\Psi_{\epsilon}^{in}\big{(}a,b,c\big{)}=\begin{cases}0,&b<a<c\\\
(b-a)/\epsilon,&a<b<a+\epsilon<c\\\ 1,&a+\epsilon<b<c-\epsilon\\\
(c-b)/\epsilon,&a<c-\epsilon<b<c\\\ 0,&a<c<b\\\ \end{cases}\\\ \text{, and}$

$\Psi_{\epsilon}^{out}\big{(}a,b,c\big{)}=1-\Psi_{\epsilon}^{in}\big{(}a,b,c\big{)}.$

Figure 1: Loss Functions $\Psi_{\epsilon}$,$\Psi_{\epsilon,1}$ and
$\Psi_{\epsilon,2}$

###  3.2 DC Algorithm for One-sided PDI

Assume that $f_{L}$ belongs to a reproducing kernel Hilbert space (RKHS)
$\mathcal{H}$ (Shawe-Taylor et al.,, 2004). We propose minimizing the
following regularized empirical risk of $R_{\epsilon}(f_{L})$ corresponding to
(5):

$\displaystyle\hat{f}_{L,n}$
$\displaystyle=\arg\min_{f_{L}\in\mathcal{H}}\sum_{i=1}^{n}\bigg{[}\Big{\\{}\alpha
I(Y_{i}\leq S)\Psi_{\epsilon}(A_{i},f_{L}(X_{i}))$
$\displaystyle+\big{(}1-\alpha\big{)}I(Y_{i}>S)\Psi_{\epsilon}(f_{L}(X_{i}),A_{i})\Big{\\}}\frac{1}{P(A_{i}|X_{i})}\bigg{]}$
$\displaystyle+\frac{\lambda_{n}}{2}\|f_{L}\|_{\mathcal{H}},$ (7)

where the regularization term
$\frac{\lambda_{n}}{2}\|f_{L}\|_{\mathcal{H}}^{2}$ penalizes the complexity of
$f_{L}$ and $\|\cdot\|_{\mathcal{H}}^{2}$ is a norm defined in $\mathcal{H}$.

Optimizing such an objective function is challenging due to the nonconvexity
of the loss function. We adopt the Difference-of-Convex-functions (DC)
algorithm (Le Thi Hoai and Tao,, 1997) to iteratively solve a series of convex
minimization problems. For simplicity, assume
$f_{L}(X_{i};\mathbf{v})=\sum_{j=1}^{n}v_{j}k(X_{i},X_{j})+v_{0}$ where $k$ is
the reproducing kernel of $\mathcal{H}$ and denote
$\mathbf{v}=(v_{1},.....,v_{n},v_{0})$ for the unknown parameters. We can
write the regularized empirical risk in (7) as

$\displaystyle\mathcal{L}(\mathbf{v})$
$\displaystyle\equiv\sum_{i=1}^{n}\bigg{[}\Big{\\{}\alpha I(Y_{i}\leq
S)\Psi_{\epsilon}\big{(}A_{i},f_{L}(X_{i};\mathbf{v})\big{)}$
$\displaystyle+\big{(}1-\alpha\big{)}I(Y_{i}>S)\Psi_{\epsilon}\big{(}f_{L}(X_{i};\mathbf{v}),A_{i}\big{)}\Big{\\}}\frac{1}{P(A_{i}|X_{i})}\bigg{]}$
$\displaystyle+\frac{\lambda_{n}}{2}\mathbf{v}^{\top}\mathbf{K}\mathbf{v}$ (8)

where $\mathbf{K}$ is a $n\times n$ matrix with $k(X_{i},X_{j})$ ($i=1\ldots
n;j=1,\ldots,n$) as its elements. Let $\mathcal{L}_{1}(\mathbf{v})$ and
$\mathcal{L}_{2}(\mathbf{v})$ be similarly defined as
$\mathcal{L}(\mathbf{v})$ except with $\Psi_{\epsilon}$ replaced with
$\Psi_{\epsilon,1}$ and $\Psi_{\epsilon,2}$ (without the regularization term)
respectively in (8). Then
$\mathcal{L}(\mathbf{v})=\mathcal{L}_{1}(\mathbf{v})-\mathcal{L}_{2}(\mathbf{v})$
and the non-convex objective function $\mathcal{L}(\mathbf{v})$ is the
difference of two convex functions $\mathcal{L}_{1}(\mathbf{v})$ and
$\mathcal{L}_{2}(\mathbf{v})$. By replacing the second function
$\mathcal{L}_{2}(\mathbf{v})$ by its first order Taylor expansion expanded
around the solution from previous iteration and the fact that a linear
function is convex and concave, the DC algorithm solves a convex optimization
problem at each iteration.

In summary, the DC algorithm estimates $\mathbf{v}$ by repeatedly solving the
following optimization problem

$\displaystyle\mathbf{v}^{t+1}$ $\displaystyle=\arg\min_{\mathbf{v}}\hskip
4.87553pt\mathcal{L}_{1}(\mathbf{v})-\\{\nabla\mathcal{L}_{2}(\mathbf{v}^{t})\\}^{\top}(\mathbf{v}-\mathbf{v}^{t})$
$\displaystyle=\arg\min_{\mathbf{v}}\hskip
4.87553pt\mathcal{L}_{1}(\mathbf{v})-\\{\nabla\mathcal{L}_{2}(\mathbf{v}^{t})\\}^{\top}\mathbf{v},$

where $\mathbf{v}^{t+1}$ and $\mathbf{v}^{t}$ are solutions of the DC-
algorithm at $t+1$-th and $t$-th iterations, respectively. The optimization
problem above is a quadratic programming problem which can be solved by
standard packages. The details of the derivations are documented in the
supplementary. The convergence property of the DC algorithm has been shown in
Lanckriet and Sriperumbudur, (2009) such that the solution is guaranteed to
converge to a local minimum.

So far, we only discussed the case of one-sided PDIs with a lower bound
function. In the case where an upper bound is desired, one could also utilize
the algorithm described above with a minor modification. By taking the
negative value of the dose, the upper bound of the original dose will become
the lower bound for the negative dose. To find two-sided PDIs, one can use the
following strategy to convert it into two separate one-sided PDI problems.
Recall that the optimal dose function $f_{opt}=\arg\max_{f}V(f)$ defined in
Section $2.2$ and assume that $f_{L,opt}(x)<f_{opt}(x)<f_{U,opt}(x)$ for all
$x\in\mathcal{X}$. Then, we can use samples with observed $A<f_{opt}(x)$ to
learn the lower bound and samples with observed $A>f_{opt}(x)$ to learn the
upper bound. To estimate the unknown function $f_{opt}(x)$, we can use the
existing methods for optimal dose finding (Chen et al.,, 2016; Kallus and
Zhou,, 2018; Zhu et al.,, 2020; Zhou et al.,, 2021). The choice of which PDI
(one-sided or two-sided) to use usually depends on practical considerations.
When recommending HbA1c level for diabetes patients, because over-control
HbA1c is the main concern, learning a one-sided PDI is sufficient (see Section
5.2 for details). In other settings, e.g. optimal personalized dose intervals
for the Warfarin drug (an anticoagulant drug), a two-sided PDI may be
preferred as overdosing on Warfarin will increase risk of blooding, and
underdosing will increase risk of stroke.

Note that PDIs depend on the threshold $S$ and the probability $\alpha$. The
combinations of these quantities with different values could lead to different
interpretations of the resulting PDIs. The pre-determined threshold $S$ that
defines the favorable outcomes can also be chosen as dependent on some
covariates $z\subset x$, e.g., age or age groups.

###  3.3 Generalized Propensity Score for Continuous Treatment

To conduct off-policy learning for a continuous treatment, we need to estimate
$P(a\,|\,x)$ which is the conditional density of dose $a$ given the covariates
$x$ (Austin,, 2018) and weigh each observation using the inverse of the
estimated $P(a\,|\,x)$. Multiple approaches are available for estimating
$P(a\,|\,x)$. One natural approach is regression. In particular, one could
assume the density is from parametric or nonparametric distributions (Austin,,
2018). This approach usually performs satisfactorily when the distributional
assumptions are reasonable. For our simulation settings as well as the data
application, the approach led to quite unstable inverse weights. To remedy the
problem, we adopt recent methods that directly estimate the weights (instead
of the inverse of estimated $P(a\,|\,x)$) that could balance the covariates
across the range of a continuous treatment (Fong et al.,, 2018; Kallus and
Santacatterina,, 2019; Tübbicke,, 2020; Huling et al.,, 2021). In particular,
we use a modified version of Kernel Optimal Orthogonality Weighting (Kallus
and Santacatterina,, 2019) and the independence-inducing weights by Huling et
al., (2021) for our numerical analysis due to their stable finite sample
performance.

##  4 THEORETICAL RESULTS

We include results for one-sided PDIs here and leave those for two-sided PDIs
to the supplement. If $[f_{L}(x),a_{U}]$ is a one-sided level-$\alpha$ PDI,
then it is reasonable to assume $[f_{L}(x)+\eta(x),a_{U}]$ is also a one-sided
level-$\alpha$ PDI for any $\eta(x)\geq 0$ as it is a narrower interval. We
ensure this by making the following assumptions throughout the section:

(C1) For any $x$, the boundaries $a_{L}$ and $a_{U}$ for the dose range
satisfy $P\\{Y(a_{L})>S\big{|}x\\}\leq\alpha$ and
$P\big{\\{}Y(a_{U})>S|x\big{\\}}\geq\alpha$.
(C2) For any $x$, there is a unique solution $a_{x}\in[a_{L},a_{U}]$ such that
$P\\{Y(a_{x})>S|x\\}=\alpha$.

The first condition (C1) ensures that our choices of $S$ and $\alpha$ are
reasonable in the sense that the level-$\alpha$ is attainable. The second
condition (C2) ensures that the limit of our learning algorithm is well
defined. Details of the proof can be found in the supplementary file. Given
(C1) and (C2), our solution $f_{L}(x)$ can attain the lowest possible value
(i.e., the PDI is as wide as possible).

Note that if $P(Y>S|A=a,X=x)$ is a monotonic function in $a$ for all $x$, i.e.
$P(Y>S|A=a_{1},X=X)\leq P(Y>S|A=a_{2},X=X)$, $\forall
a_{1},a_{2}\in\text{PDI}_{\alpha}(X)$ and $a_{1}<a_{2}$, then Conditions (C1)
and (C2) will hold with properly chosen $\alpha$ and $S(X)$. However, (C1) and
(C2) do not require $P(Y(a)>S|X=x)$ to be monotonic. In other words, the
monotonicity is sufficient but not necessary for (C1) and (C2) to hold. We
make assumptions that (C1) and (C2) hold to guarantee the existence of
$\text{PDI}_{\alpha}(X)$. Hence, our method can work with multi-modal outcomes
as long as (C1) and (C2) are satisfied. On the other hand, for many medical
applications, one can safely expect $P(Y>S|a,X)$ to be monotonic or unimodal
in $a$. A unimodal $P(Y>S|a,X)$ will require two-sided PDIs.

###  4.1 Fisher Consistency

Theorem 4.1 below shows that minimizing $R(\cdot)$ defined in (4) will result
in a one-sided PDI lower bound, above which the doses will produce outcomes
larger than $S$ with a probability of at least $\alpha$.

######  Theorem 4.1.

Let $f_{L,opt}=\arg\min_{f}R(f)$. Then, under the assumptions (C1) and (C2),
$f_{L,opt}(x)$ has the following two properties:
(1) $P\big{\\{}Y\big{(}f_{L,opt}(x)\big{)}>S\,|\,x\big{\\}}=\alpha,\forall x$.
(2) For any measurable function $a(\cdot)$ s.t. $f_{L,opt}(x)\leq a(x)\leq
a_{U}$, $\forall x\in\mathcal{X}$, the potential outcome $Y(a(x))$ satisfies
$P(Y(a(x))>S|x)\geq\alpha$.

Theorem 4.1 shows that for any given value of $X=x$, the length of
$[f_{L,opt}(x),a_{U}]$ is the widest among all valid $\text{PDI}_{\alpha}(x)$.

###  4.2 Convergence Rate

Theorem 4.2 below shows that, for one-sided PDI lower bounds, the difference
between the risk function $R(\cdot)$ and its non-convex relaxation
$R_{\epsilon_{n}}(\cdot)$ converges to zero with the same rate as
$\epsilon_{n}\rightarrow 0$. That is, the approximation bias goes to $0$. Then
in Theorem 4.3 we derive the convergence rates of the risk function based on
the estimated $\hat{f}_{L,n}$ and the general results of empirical risk
minimization and approximation. This theorem tells us how the empirically
optimal policy performs out-of-sample. Here, we assume the Gaussian kernel
with bandwidth $\gamma_{n}$ (i.e.
$\mathbf{K}(x,x^{\prime})=\exp(-\gamma_{n}^{2}||x-x^{\prime}||^{2}))$ was used
in estimating $\hat{f}_{L,n}$. Since ${f}_{L,opt}$ may not belong to the RKHS
we optimize in, we choose the Gaussian kernel (which is dense in $L_{2}$) to
make sure $\hat{f}_{L,n}$ can well approximate ${f}_{L,opt}$.

######  Theorem 4.2.

For any measurable function $f:\mathcal{X}\mapsto\mathbb{R}$, we have
$|R(f)-R_{\epsilon_{n}}(f)|\leq C\epsilon_{n}$, where $C$ is a constant.

######  Theorem 4.3.

Let $B^{\delta}_{1,\infty}(\mathbb{R}^{d})$ be a Besov space defined as
$B^{\delta}_{1,\infty}(\mathbb{R}^{d})=\\{f\in
L_{\infty}(\mathbb{R}^{d}):\sup_{t>0}(t^{-\delta}w_{(r,L_{1})}(f,t))<\infty\\}$,
where $w$ is the modulus of continuity. Assume that $f_{L,opt}\in
B^{\delta}_{1,\infty}(\mathbb{R}^{d})$. Then for any $\zeta>0$,
$d/(d+\tau)<p<1$, and $\tau>0$,

$\displaystyle R(\hat{f}_{L,n})-R(f_{L,opt})\leq
c_{1}\frac{\lambda_{n}}{\gamma_{n}^{d}}+c_{2}\gamma_{n}^{\delta}+$
$\displaystyle
c_{3}{\gamma_{n}^{-\frac{(1-p)(1+\zeta)d}{2-p}}\lambda_{n}^{-\frac{p}{2-p}}n^{-\frac{1}{2-p}}}+c_{4}\frac{\tau^{\frac{1}{2}}}{n^{\frac{1}{2}}}+c_{5}\frac{\tau}{n}+c_{6}\epsilon_{n}$

with probability $1-3e^{-\tau}$. Furthermore, if we choose these constants to
be

$\displaystyle\gamma_{n}=\mathcal{O}\bigg{(}\frac{1}{n}\bigg{)}^{\frac{1}{2\delta+d}}\hskip
6.68129pt\lambda_{n}=\mathcal{O}\bigg{(}\frac{1}{n}\bigg{)}^{\frac{\delta+d}{2\delta+d}}\hskip
6.68129pt\epsilon_{n}=\mathcal{O}\bigg{(}n^{-\frac{\delta}{2\delta+d}}\bigg{)}.$

We have the following optimal convergence rate:
$R(\hat{f}_{L,n})-R(f_{L,opt})=\mathcal{O}\bigg{(}n^{-\frac{\delta}{2\delta+d}}\bigg{)}.$

Theorem 4.3 shows that the O-learning for PDI can potentially achieve a rate
of nearly $\mathcal{O}(n^{-\frac{1}{2}})$, when the smoothness parameter
$\delta$ of the optimal bound function $f_{L,opt}$ is relatively large
compared to the covariate dimension $d$. This rate is faster than the
convergence rate of about $\mathcal{O}(n^{-\frac{1}{4}})$ in Chen et al.,
(2016), which focused on finding the optimal dose. This improvement is due to
the fact that more sample points are contained in the neighborhood of the
entire optimal interval than in the neighborhood of the optimal dose. Note
that, we have assumed the generalized propensity score $P(a|x)$ can be
estimated at the rate of $\mathcal{O}(n^{-1/2})$ or faster. Consequently, the
estimation error of the generalized propensity score is negligible. However,
if the generalized propensity score can not be estimated at this rate, then
convergence rate in Theorem 4.3 will be dominated by the possibly slower
convergence rate of the generalized propensity score estimator.

##  5 EXPERIMENTS

###  5.1 Simulation

We compare our methods with three indirect methods for learning one-sided
$\text{PDI}_{\alpha}(X)$. Note that existing optimal dose learning methods can
not yield a dose interval, hence they are not included in the comparison. For
our methods, we implemented linear and Gaussian kernels and denote them as
L-O-Learning and G-O-Learning. As the DC-algorithm is an iterative algorithm
which relies on initial values, we used covariates-independent lower bounds as
such initial values. All indirect methods involve classification of outcomes
to find the PDIs. In particular, the indicator $I{\\{Y>S\\}}$ was used as the
binary outcome for classification so that we could build a model for
$\hat{P}(Y>S|a,X)\equiv h(A,X)$. Then we found the PDIs as a range of doses
$\mathcal{A}_{1}(X)$ such that for $\forall a\in\mathcal{A}_{1}(X)$,
$h(a,X)>\alpha$. To find $\mathcal{A}_{1}(X)$, we evaluated different $a$
value in estimated $h(a,X)$. All $a$ that satisfied $h(a,X)>\alpha$ were
included in $\mathcal{A}_{1}(X)$. In our simulations, we used $200$ evenly
spaced grids to construct $\mathcal{A}_{1}(X)$.

We used logistic regression with $L_{1}$ penalty, Support Vector Machine
(SVM), and Random Forests (RF) as indirect methods. For the penalized logistic
regression, we included linear and quadratic terms of the dose, covariates,
and pairwise interactions between dose and covariates as predictors. For SVM,
we used the Gaussian kernel with the kernel bandwidth fixed using the “median
heuristic”. For RF, we set the number of trees to be $1000$. All methods
required parameter tuning and we tuned the $L_{2}$ regularization parameters
for O-learnings, the $L_{1}$ regularization parameter for the logistic
regression, the margin trade-off parameter for the SVM, and the percentage of
variables that can be used for splitting at each node for the RF. Tuning
parameters for all methods were selected using 5-fold cross-validation.

We considered multiple simulation settings with various data generating
processes, sample sizes, numbers of covariates as well as signal to noise
levels. We set the confidence probability $\alpha$ to be $0.5$. Scenario $1$
is for linear lower bounds and Scenario $2$ for nonlinear lower bounds. For
both scenarios, there is a unique value of $a$ such that
$P(Y(a)>S(x)|x)=\alpha$ for a given $x$.

We assume that the training data contain $n$ i.i.d. observations
$(A_{i},Y_{i},X_{i}),i=1\ldots n$. The dose $A_{i}$ is generated from a
truncated normal distribution $\mathcal{N}(\mu_{A},0.5)$ with the truncation
limits $[-2,2]$. $X_{i}$ is generated from a uniform distribution
$\text{Uniform}(-1,1)^{d}$. The outcome $Y$ is generated from a normal
distribution $\mathcal{N}(\mu_{Y},\sigma^{2})$. We assume there exist
confounders (i.e., $X$ can impact both treatment assignment and the outcome)
to mimic the observational study settings. In the supplement, we also include
simulations for randomized trial settings.

Scenario 1: Linear bounds
$\mu_{A}=0.3X_{1}+0.3X_{2}+0.3X_{3};$
$\mu_{Y}={5e^{10A}}/(e^{10\mu_{A}}+e^{10A})+D(X);$
$D(X)=0.6X_{2}+0.6X_{3}+0.6X_{4}.$

Scenario 2: Non-linear bounds
$\mu_{A}=0.75\log(|X_{1}|+1)-0.2\cos(\pi X_{2})+0.2I(X_{3}>0)-0.4;$
$\mu_{Y}={5e^{10A}}/(e^{10C}+e^{10A})+D(X);$
$D(X)=0.4\sin(\pi X_{2})+0.4I(X_{3}>0)+0.4|X_{4}|.$

We repeat each setting for $100$ times and report the results of all methods
by evaluating the empirical risk or the empirical value of $R(\hat{f}_{L})$
defined in Equation (4) on corresponding independent testing data sets with
sample sizes of $10000$. For O-learning methods, we estimate weights that
could approximate $1/P(a|x)$ using the approaches mentioned in the Section
3.3. Meanwhile when evaluating the performances of all methods, we plug in the
true values of $1/P(a|x)$ for calculating the empirical risks of estimated
policies on the testing data sets.

Table 1: Empirical risk ($R(\hat{f}_{L})$) with SD in (.) from comparison methods under different settings | n | d | L-O-Learning | G-O-Learning | Logistic | SVM | RF
---|---|---|---|---|---|---|---
Scenario 1 | 200 | 10 | 0.048 (0.004) | 0.049 (0.004) | 0.054 (0.023) | 0.053 (0.005) | 0.053 (0.005)
$\sigma^{2}=2.25$ | 200 | 50 | 0.055 (0.005) | 0.052(0.002) | 0.176(0.074) | 0.077 (0.041) | 0.058 (0.004)
| 400 | 10 | 0.045 (0.002) | 0.047 (0.002) | 0.044(0.002) | 0.048 (0.002) | 0.049 (0.002)
| 400 | 50 | 0.056 (0.004) | 0.051 (0.002) | 0.071(0.025) | 0.055 (0.003) | 0.053 (0.002)
Scenario 1 | 200 | 10 | 0.122 (0.004) | 0.121 (0.003) | 0.182(0.049) | 0.139 (0.026) | 0.134 (0.010)
$\sigma^{2}=9$ | 200 | 50 | 0.132 (0.005) | 0.122(0.002) | 0.215(0.049) | 0.238 (0.032) | 0.142 (0.014)
| 400 | 10 | 0.120 (0.003) | 0.121 (0.003) | 0.142(0.037) | 0.124 (0.007) | 0.128 (0.005)
| 400 | 50 | 0.129 (0.003) | 0.122 (0.002) | 0.205 (0.049) | 0.139 (0.023) | 0.127 (0.004)
Scenario 2 | 200 | 10 | 0.054 (0.003) | 0.050 (0.002) | 0.069 (0.044) | 0.055 (0.005) | 0.052 (0.005)
$\sigma^{2}=2.25$ | 200 | 50 | 0.061 (0.005) |  0.051 (0.002) | 0.187 (0.078) | 0.177 (0.081) | 0.057 (0.006)
| 400 | 10 | 0.051 (0.002) | 0.049 (0.002) | 0.051 (0.010) | 0.051 (0.004) | 0.049 (0.005)
| 400 | 50 | 0.061 (0.003) | 0.051 (0.002) | 0.095(0.050) | 0.060 (0.003) | 0.051 (0.003)
Scenario 2 | 200 | 10 | 0.123 (0.003) |  0.120 (0.003) | 0.177 (0.052) | 0.144 (0.035) | 0.131 (0.011)
$\sigma^{2}=9$ | 200 | 50 | 0.133 (0.005) |  0.121 (0.002) | 0.199 (0.053) | 0.249 (0.010) | 0.145 (0.016)
| 400 | 10 | 0.121 (0.003) | 0.120 (0.002) | 0.143 (0.038) | 0.123 (0.005) | 0.127 (0.009)
| 400 | 50 | 0.131 (0.004) | 0.121 (0.002) | 0.198 (0.053) | 0.153 (0.040) | 0.125 (0.005)

For all methods we compared, we set $S$ to be the fitted value of a polynomial
regression of Y with X as predictors (i.e., set
$S(x)=\hat{E}\left[Y|x\right]$), as there are strong prognostic effects. From
Table 1, the L-O-Learning has advantages in the linear bound settings and the
G-O-Learning has advantages in the nonlinear settings. In particular, compared
to other methods, our methods are more robust to high covariate dimensionality
(large $d$), low sample size (small $n$), and low signal-to-noise ratio (large
$\sigma^{2}$) settings. In terms of computational time, the RF and L-O-
Learning/G-O-Learning are comparable, and they are much faster than the SVM.

###  5.2 Study of HbA1c Control for Diabetic Patients

A fundamental aspect of managing diabetes is achieving tight control of blood
sugar levels to reduce the risk of diabetes complications. Current diabetes
guidelines recognize that treatment recommendations regarding tight control of
HbA1c may not be appropriate for older patients, especially those with
comorbid conditions, due to increased risk for adverse events (American
Diabetes Association,, 2018). For example, complex patients may be at risk for
negative outcomes when treated aggressively (Ismail-Beigi et al.,, 2011;
Riddle and Karl,, 2012). These increased risks have even more significance if
the benefits of tight control are limited due to shortened life expectancy. If
the long-term benefits are not relevant and the short-term risks are
significant, there may be little reason to pursue tight control of HbA1c. As a
proof-of-concept analysis, we would like to recommend one-sided PDI for the
upper bound of HbA1c such that the number hospitalization events are
minimized.

We use a de-identified electronic health record (EHR) dataset from a large
health care system which contains laboratory records, information on primary
care visits, and medication history of diabetes patients between 01/2003 to
12/2011. We use $45$ relevant features and $8126$ subjects whose observed
HbA1c measures that are between $4\%$ and $10\%$ to learn the PDI of HbA1c.
Details about the data set can be found in the supplement. For the purpose of
notation consistency, we take the negative value of number of hospitalization
events as the outcome. Hence a larger value of $Y$ represents a preferred
outcome. By this coding, we also set $S=-0.5$ and $\alpha=0.5$, which leads to
the interpretation of PDI such that if a patient controls his/her HbA1c less
than the upper bound, they would likely not experience diabetes-related
hospitalization (i.e., with the probability of hospitalization less than
$0.5$) in the next 90-days after HbA1c is measured. Since the observed outcome
is a count, any choice of $S$ that is between $-1$ and $0$ will lead to the
same result and interpretation as setting $S=-0.5$.

We randomly select $1000$ patients as a training data set to learn the policy
and use the remaining $7126$ patients as a testing data set. We repeat this
random splitting $100$ times and report the averaged results. We apply all
methods compared in our simulations to this real data. The tuning parameters
are selected by using the same 5-fold cross-validation strategy.

To demonstrate the result, we plotted the distributions of the observed HbA1c
level and the estimated upper bounds of samples in one testing data set in
Figure 2. While the majority of patients have observed HbA1c between 5.5% and
7.5%, their estimated upper bounds of the L-O-Learning and G-O-Learning are
generally between 6.5% and 8%, which suggests that some patients may be over-
controlling their HbA1c levels. In comparison, the PDI recommendations from
the comparison methods are highly tilted towards the upper limit, especially
the SVM. It is unlikely that physicians will recommend such high HbA1c values
to patients. Furthermore, the grid search range of HbA1c levels for the
indirect methods is between 4% and 10%, and indirect methods recommended a
high proportion of patients to have 10% as the HbA1c levels upper bound
(logistic regression: 30% of all patients; SVM: 70%; RF: 50%). If we change
the grid search range to be between 4% and 9%, the indirect methods would
recommend a high proportion of patients to have 9% as the upper bound. This
indicates the policies estimated from indirect methods can have unwanted
dependence on the grid search range. Such patterns of the estimated upper
bounds of HbA1c levels shown in Figure 2 are similar across all $100$
replications.

Figure 2: Density plot of observed and estimated upper bounds of HbA1c PDIs.

Furthermore, the upper bounds learned from our methods (especially the G-O-
Learning) have much better agreement with the current medical guidance
compared to other methods. In particular, the 2018 American Diabetes
Association (ADA) Standards of Medical Care in Diabetes advises the following
upper HbA1c bounds for diabetes patients (American Diabetes Association,,
2018): 6.5% for people who can achieve this goal without experiencing a lot of
hypoglycemia episodes or other negative effects of having lower blood glucose
levels; 7.5% for all children with diabetes; 8% for people with a history of
severe hypoglycemia, people who have had diabetes for many years or have
difficulty achieving tighter control, and people with limited remaining life
expectancy. In our dataset, 88% of the patients are over 60 years old, 50%
over 73 years old, and 20% over 80 years old. Most of these patients have
diabetes for many years. Hence, the distributions of upper bounds
recommendations from our methods coincide well with the standards given the
age distribution.

##  6 CONCLUSION

We have proposed a framework for learning individualized dose intervals with
desired theoretical properties. Our method extends a previous method that can
only learn the policy for an optimal dosage (Chen et al.,, 2016). We proposed
a DC algorithm to efficiently solve the underlying non-convex optimization
problems with linear or nonlinear policies learning.

Note that both direct and indirect methods proposed by us can learn PDI,
meanwhile, we advocate using the direct method since it can recommend dose
intervals without the need of fitting complex regression models. Learning PDI
relies on the function restriction due to Conditions C1 and C2 in Section 4,
and such a restriction applies to both direct and indirect methods. If the
conditions are not satisfied, the PDI may not be well defined or it will be a
union of several discontinuous intervals.

We also believe that $S$ should be chosen to have good interpretation. In our
HbA1c example, choosing $S$ between $-1$ to $0$ leads to the same
interpretation of the resulting PDI. For the applications where there is no
obvious choice of $S$, we could set it to be a number between $\alpha\times
100\%$-th percentile of $Y(a_{L})|X$ and $\alpha\times 100\%$-th percentile of
$Y(a_{U})|X$ to ensure the existence of PDI.

Lastly, learning causally meaningful PDIs using observed data relies on the
no-unmeasured confounder assumption. In our HbA1c example, we worked closely
with our physician collaborator to include curated variables from EHR and
medical claims to ensure this assumption. At the same time, sensitivity
analysis may help examine the robustness of the PDI policy. We want to
emphasize that before implementing an estimated policy for PDI in practice, it
is critically important to verify the superiority of the estimated policy over
existing clinical practice via randomized trials.

### Acknowledgments

We would like to thank all reviewers for their helpful comments that greatly
improved this paper. Research reported in this work was partially funded
through a Patient-Centered Outcomes Research Institute (PCORI) Award
(ME-2018C2-13180). The views in this work are solely the responsibility of the
authors and do not necessarily represent the views of the Patient-Centered
Outcomes Research Institute (PCORI), its Board of Governors or Methodology
Committee.

## References

  * American Association of Clinical Endocrinologists, (2019) American Association of Clinical Endocrinologists (2019).  Management of common comorbidities of diabetes.
  * American Diabetes Association, (2018) American Diabetes Association (2018).  Standards of medical care in diabetes—2018 abridged for primary care providers.  Clinical diabetes: a publication of the American Diabetes Association, 36(1):14.
  * Austin, (2018) Austin, P. C. (2018).  Assessing the performance of the generalized propensity score for estimating the effect of quantitative or continuous exposures on binary outcomes.  Statistics in Medicine, 37(11):1874–1894.
  * Bica et al., (2020) Bica, I., Jordon, J., and van der Schaar, M. (2020).  Estimating the effects of continuous-valued interventions using generative adversarial networks.  arXiv preprint arXiv:2002.12326.
  * Borghans et al., (2012) Borghans, I., Kleefstra, S. M., Kool, R. B., and Westert, G. P. (2012).  Is the length of stay in hospital correlated with patient satisfaction?  International Journal for Quality in Health Care, 24(5):443–451.
  * Breiman, (2001) Breiman, L. (2001).  Random forests.  Machine learning, 45(1):5–32.
  * Breiman et al., (1984) Breiman, L., Friedman, J., Stone, C. J., and Olshen, R. A. (1984).  Classification and regression trees.  CRC press.
  * Chakraborty and Moodie, (2013) Chakraborty, B. and Moodie, E. (2013).  Statistical methods for dynamic treatment regimes.  Springer.
  * 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.
  * Collins and Varmus, (2015) Collins, F. S. and Varmus, H. (2015).  A new initiative on precision medicine.  New England journal of medicine, 372(9):793–795.
  * Eberts et al., (2013) Eberts, M., Steinwart, I., et al. (2013).  Optimal regression rates for svms using gaussian kernels.  Electronic Journal of Statistics, 7:1–42.
  * Flack and Adekola, (2020) Flack, J. M. and Adekola, B. (2020).  Blood pressure and the new acc/aha hypertension guidelines.  Trends in cardiovascular medicine, 30(3):160–164.
  * Fong et al., (2018) Fong, C., Hazlett, C., Imai, K., et al. (2018).  Covariate balancing propensity score for a continuous treatment: Application to the efficacy of political advertisements.  The Annals of Applied Statistics, 12(1):156–177.
  * Freedman et al., (2010) Freedman, B. I., Shenoy, R. N., Planer, J. A., Clay, K. D., Shihabi, Z. K., Burkart, J. M., Cardona, C. Y., Andries, L., Peacock, T. P., Sabio, H., et al. (2010).  Comparison of glycated albumin and hemoglobin a1c concentrations in diabetic subjects on peritoneal and hemodialysis.  Peritoneal Dialysis International, 30(1):72–79.
  * Huling et al., (2021) Huling, J. D., Greifer, N., and Chen, G. (2021).  Independence weights for causal inference with continuous exposures.  arXiv preprint arXiv:2107.07086.
  * Imbens and Rubin, (2015) Imbens, G. W. and Rubin, D. B. (2015).  Causal inference in statistics, social, and biomedical sciences.  Cambridge University Press.
  * International Warfarin Pharmacogenetics Consortium, (2009) International Warfarin Pharmacogenetics Consortium (2009).  Estimation of the warfarin dose with clinical and pharmacogenetic data.  New England Journal of Medicine, 360(8):753–764.
  * Ismail-Beigi et al., (2011) Ismail-Beigi, F., Moghissi, E., Tiktin, M., Hirsch, I. B., Inzucchi, S. E., and Genuth, S. (2011).  Individualizing glycemic targets in type 2 diabetes mellitus: implications of recent clinical trials.  Annals of internal medicine, 154(8):554–559.
  * Kallus and Santacatterina, (2019) Kallus, N. and Santacatterina, M. (2019).  Kernel optimal orthogonality weighting: A balancing approach to estimating effects of continuous treatments.  arXiv preprint arXiv:1910.11972.
  * Kallus and Zhou, (2018) Kallus, N. and Zhou, A. (2018).  Policy evaluation and optimization with continuous treatments.  In International Conference on Artificial Intelligence and Statistics, pages 1243–1251.
  * Kosorok and Laber, (2019) Kosorok, M. R. and Laber, E. B. (2019).  Precision medicine.  Annual Review of Statistics and Its Application, 6:263–286.
  * Lanckriet and Sriperumbudur, (2009) Lanckriet, G. R. and Sriperumbudur, B. K. (2009).  On the convergence of the concave-convex procedure.  In Advances in Neural Information Processing Systems, pages 1759–1767.
  * Le Thi Hoai and Tao, (1997) Le Thi Hoai, A. and Tao, P. D. (1997).  Solving a class of linearly constrained indefinite quadratic problems by DC algorithms.  Journal of global optimization, 11(3):253–285.
  * Liang and Yu, (2020) Liang, M. and Yu, M. (2020).  A semiparametric approach to model effect modification.  Journal of the American Statistical Association, 0(0):1–13.
  * Little et al., (2011) Little, R. R., Rohlfing, C. L., and Sacks, D. B. (2011).  Status of hemoglobin a1c measurement and goals for improvement: from chaos to order for improving diabetes care.  Clinical chemistry, 57(2):205–214.
  * Lou et al., (2018) Lou, Z., Shao, J., and Yu, M. (2018).  Optimal treatment assignment to maximize expected outcome with multiple treatments.  Biometrics, pages 506–516.
  * Luedtke and van der Laan, (2016) Luedtke, A. R. and van der Laan, M. J. (2016).  Comment on: ‘personalized dose finding using outcome weighted learning’.  Journal of the American Statistical Association, 111(516):1526–1530.
  * Mancinelli et al., (2000) Mancinelli, L., Cronin, M., and Sadée, W. (2000).  Pharmacogenomics: the promise of personalized medicine.  Aaps Pharmsci, 2(1):29–41.
  * Moodie et al., (2012) Moodie, E. E., Chakraborty, B., and Kramer, M. S. (2012).  Q-learning for estimating optimal dynamic treatment rules from observational data.  Canadian Journal of Statistics, 40(4):629–645.
  * Moodie et al., (2009) Moodie, E. E., Platt, R. W., and Kramer, M. S. (2009).  Estimating response-maximized decision rules with applications to breastfeeding.  Journal of the American Statistical Association, 104(485):155–165.
  * Nitin, (2010) Nitin, S. (2010).  Hba1c and factors other than diabetes mellitus affecting it.  Singapore Med J, 51(8):616–622.
  * Radin, (2014) Radin, M. S. (2014).  Pitfalls in hemoglobin a1c measurement: when results may be misleading.  Journal of general internal medicine, 29(2):388–394.
  * Rich et al., (2014) Rich, B., Moodie, E. E., and Stephens, D. A. (2014).  Simulating sequential multiple assignment randomized trials to generate optimal personalized warfarin dosing strategies.  Clinical trials, 11(4):435–444.
  * Riddle and Karl, (2012) Riddle, M. C. and Karl, D. M. (2012).  Individualizing targets and tactics for high-risk patients with type 2 diabetes: practical lessons from accord and other cardiovascular trials.  Diabetes Care, 35(10):2100–2107.
  * Schulz and Moodie, (2020) Schulz, J. and Moodie, E. E. (2020).  Doubly robust estimation of optimal dosing strategies.  Journal of the American Statistical Association, pages 1–13.
  * Scott et al., (2017) Scott, J. G., Berglund, A., Schell, M. J., Mihaylov, I., Fulp, W. J., Yue, B., Welsh, E., Caudell, J. J., Ahmed, K., Strom, T. S., et al. (2017).  A genome-based model for adjusting radiotherapy dose (gard): a retrospective, cohort-based study.  The lancet oncology, 18(2):202–211.
  * Shawe-Taylor et al., (2004) Shawe-Taylor, J., Cristianini, N., et al. (2004).  Kernel methods for pattern analysis.  Cambridge university press.
  * Sondhi et al., (2020) Sondhi, A., Arbour, D., and Dimmery, D. (2020).  Balanced off-policy evaluation in general action spaces.  In International Conference on Artificial Intelligence and Statistics, pages 2413–2423. PMLR.
  * Steinwart and Christmann, (2008) Steinwart, I. and Christmann, A. (2008).  Support vector machines.  Springer Science & Business Media.
  * Tibshirani, (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society. Series B (Methodological), pages 267–288.
  * Tübbicke, (2020) Tübbicke, S. (2020).  Entropy balancing for continuous treatments.  arXiv preprint arXiv:2001.06281.
  * Vapnik, (2013) Vapnik, V. (2013).  The nature of statistical learning theory.  Springer science & business media.
  * Wu and Liu, (2007) Wu, Y. and Liu, Y. (2007).  Robust truncated hinge loss support vector machines.  Journal of the American Statistical Association, 102(479):974–983.
  * Xu et al., (2015) Xu, Y., Yu, M., Zhao, Y.-Q., Li, Q., Wang, S., and Shao, J. (2015).  Regularized outcome weighted subgroup identification for differential treatment effects.  Biometrics, 71(3):645–653.
  * Zhao et al., (2009) Zhao, Y., Kosorok, M. R., and Zeng, D. (2009).  Reinforcement learning design for cancer clinical trials.  Statistics in medicine, 28(26):3294–3315.
  * Zhou et al., (2021) Zhou, W., Zhu, R., and Zeng, D. (2021).  A parsimonious personalized dose-finding model via dimension reduction.  Biometrika, 108(3):643–659.
  * Zhu et al., (2020) Zhu, L., Lu, W., Kosorok, M. R., and Song, R. (2020).  Kernel assisted learning for personalized dose finding.  In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’20, page 56–65.

Supplementary Material:
Policy Learning for Optimal Individualized Dose Intervals

##  Appendix A DC ALGORITHM FOR ONE-SIDED PDI

Let $f_{L}(X_{i};\mathbf{v})=\sum_{j=1}^{n}v_{j}k(X_{i},X_{j})+v_{0}$ and
$\mathbf{v}=(v_{1},.....,v_{n},v_{0})$. As claimed in Section 3.2 in the main
paper, the DC algorithm repeatedly updates

$\begin{split}\mathbf{v}^{t+1}=\arg\min_{\mathbf{v}}\left(\mathcal{L}_{1}(\mathbf{v})-[\nabla\mathcal{L}_{2}(\mathbf{v}^{t})]^{\top}(\mathbf{v}-\mathbf{v}^{t})\right).\end{split}$
(9)

We derive the form of the quadratic programming by first showing that

$\begin{split}\nabla\mathcal{L}_{2}(\mathbf{v})&={\nabla}_{\mathbf{v}}\Bigg{(}\sum_{i=1}^{n}\frac{1}{P(A_{i}|X_{i})}\Big{[}\big{(}1-\alpha\big{)}I(Y_{i}>S)\Psi_{\epsilon,2}\big{(}f_{L}(X_{i}),A_{i}\big{)}+\alpha
I(Y_{i}\leq
S)\Psi_{2,\epsilon}\big{(}A_{i},f_{L}(X_{i})\big{)}\Big{]}\Bigg{)}\\\
&=\sum_{i=1}^{n}\frac{1}{P(A_{i}|X_{i})}\Big{[}\big{(}1-\alpha\big{)}I(Y_{i}>S)Q_{i}(\beta)-\alpha
I(Y_{i}\leq
S)\tilde{Q}_{i}(\beta)\Big{]}\nabla_{\mathbf{v}}f_{L}(X_{i})\end{split}$

where
$\nabla_{\mathbf{v}}f_{L}(X_{i})=(k(X_{i},X_{1}),\ldots,k(X_{i},X_{n}),1)$,
$Q_{i}(\mathbf{v})=I\big{(}f_{L}(X_{i};\mathbf{v})-A_{i}>\epsilon\big{)}/\epsilon$
and
$\tilde{Q}_{i}(\mathbf{v})=I\big{(}A_{i}-f_{L}(X_{i};\mathbf{v})>\epsilon\big{)}/\epsilon$

Following that, we have

$\displaystyle[\nabla\mathcal{L}_{2}(\mathbf{v}^{t})]^{\top}\mathbf{v}$
$\displaystyle=\sum_{i=1}^{n}\frac{1}{P(A_{i}|X_{i})}\Big{[}\big{(}1-\alpha\big{)}I(Y_{i}>S)Q_{i}(\mathbf{v}^{t})-\alpha
I(Y_{i}\leq
S)\tilde{Q}_{i}(\mathbf{v}^{t})\Big{]}\big{\\{}\sum_{j=1}^{n}v_{j}k(X_{i},X_{j})+v_{0}\big{\\}}$

By plugging in the above results into the DC procedure,

$\displaystyle\mathbf{v}^{t+1}$ $\displaystyle=\arg\min_{\mathbf{v}}\hskip
9.75853pt\mathcal{L}_{1}(\mathbf{v})-[\nabla\mathcal{L}_{2}(\mathbf{v}^{t})]^{\top}(\mathbf{v}-\mathbf{v}^{t})=\arg\min_{\beta}\hskip
9.75853pt\mathcal{L}_{1}(\mathbf{v})-[\nabla\mathcal{L}_{2}(\mathbf{v}^{t})]^{\top}\mathbf{v}$
$\displaystyle=\arg\min_{\mathbf{v}}\hskip
9.75853pt\sum_{i=1}^{n}\frac{1}{P(A_{i}|X_{i})}\Bigg{(}\big{(}1-\alpha\big{)}I(Y_{i}>S)\Big{[}\Psi_{1,\epsilon}\big{(}f_{L}(X_{i};\mathbf{v}),A_{i}\big{)}-Q_{i}(\mathbf{v}^{t})(\sum_{j=1}^{n}v_{j}k(X_{i},X_{j})+v_{0})\Big{]}$
$\displaystyle\hskip 9.75853pt+\alpha I(Y_{i}\leq
S)\Big{[}\Psi_{1,\epsilon}\big{(}A_{i},f_{L}(X_{i};\mathbf{v})\big{)}+\tilde{Q}_{i}(\mathbf{v}^{t})(\sum_{j=1}^{n}v_{j}k(X_{i},X_{j})+v_{0})\Big{]}\Bigg{)}+\frac{\lambda_{n}}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}v_{i}v_{j}k(X_{i},X_{j})$
(10)

To solve the above optimization problem, we now show that it is essentially a
quadratic programming problem which can be easily solved using standard
optimization packages.

Let $H_{i}=\frac{(1-\alpha)I(Y_{i}>S)}{P(A_{i}|X_{i})\lambda_{n}}$ and
$\tilde{H}_{i}=\frac{\alpha(Y_{i}\leq S)}{P(A_{i}|X_{i})\lambda_{n}}$.
Minimizing Equation (10) is equivalent to

$\displaystyle\arg\min_{\mathbf{v}}$
$\displaystyle\sum_{i=1}^{n}H_{i}\xi_{i}+\sum_{i=1}^{n}\tilde{H}_{i}\tilde{\xi}_{i}+\sum_{i=1}^{n}(\tilde{H}_{i}\tilde{Q}_{i}^{t}-H_{i}Q_{i}^{t})(\sum_{j=1}^{n}v_{j}k(X_{i},X_{j})+v_{0})+\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}v_{i}v_{j}k(X_{i},X_{j})$
$\displaystyle s.t.\;\;\;$
$\displaystyle\xi_{i}-[(\sum_{j=1}^{n}v_{j}k(X_{i},X_{j})+v_{0})-A_{i}]/\epsilon\geq
0,\hskip 48.78522pt\xi_{i}\geq 0$
$\displaystyle\tilde{\xi}_{i}-[A_{i}-(\sum_{j=1}^{n}v_{j}k(X_{i},X_{j})+v_{0})]/\epsilon\geq
0,\hskip 48.78522pt\tilde{\xi}_{i}\geq 0$

Applying the Lagrangian multiplier, the problem above is equivalent to

$\displaystyle\max_{\alpha_{i},\tilde{\alpha}_{i},\mu_{i},\tilde{\mu}_{i}}\min_{\mathbf{v},\xi}\sum_{i=1}^{n}H_{i}\xi_{i}+\sum_{i=1}^{n}\tilde{H}_{i}\tilde{\xi}_{i}+\sum_{i=1}^{n}(\tilde{H}_{i}\tilde{Q}_{i}^{t}-H_{i}Q_{i}^{t})(\sum_{j=1}^{n}v_{j}k(X_{i},X_{j})+v_{0})+\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}v_{i}v_{j}k(X_{i},X_{j})$
$\displaystyle-\sum_{i=1}^{n}\alpha_{i}(\xi_{i}-[(\sum_{j=1}^{n}v_{j}k(X_{i},X_{j})+v_{0})-A_{i}]/\epsilon)-\sum_{i=1}^{n}\tilde{\alpha}_{i}(\tilde{\xi}_{i}-[A_{i}-(\sum_{j=1}^{n}v_{j}k(X_{i},X_{j})+v_{0})]/\epsilon)-\sum_{i=1}^{n}\mu_{i}\xi_{i}-\sum_{i=1}^{n}\tilde{\mu}_{i}\tilde{\xi}_{i}$
$\displaystyle s.t.$ $\displaystyle\hskip
24.39261pt\sum_{i=1}^{n}v_{i}k(X_{i},X_{j})+\sum_{i=1}^{n}(\tilde{H}_{i}\tilde{Q}_{i}^{t}-H_{i}Q_{i}^{t}+\alpha_{i}/\epsilon-\tilde{\alpha}_{i}/\epsilon)k(X_{i},X_{j})=0\;\;\;\forall
j=1,\dots,n$ $\displaystyle\hskip
24.39261pt\sum_{i=1}^{n}(\tilde{H}_{i}\tilde{Q}_{i}^{t}-H_{i}Q_{i}^{t}+\alpha_{i}/\epsilon-\tilde{\alpha}_{i}/\epsilon)=0$
$\displaystyle\hskip 24.39261ptH_{i}-\alpha_{i}-\mu_{i}=0\;\;\;\forall
i=1,\dots,n$ $\displaystyle\hskip
24.39261pt\tilde{H}_{i}-\tilde{\alpha}_{i}-\tilde{\mu}_{i}=0\;\;\;\forall
i=1,\dots,n$ (11)

Incorporating the constraint that $\alpha_{i}\geq 0,\;\tilde{\alpha}_{i}\geq
0,\;\mu_{i}\geq 0,\;\tilde{\mu}_{i}\geq 0$ and denoting
$\beta_{i}=\tilde{H}_{i}\tilde{Q}_{i}^{t}-H_{i}Q_{i}^{t}+\alpha_{i}/\epsilon-\tilde{\alpha}_{i}/\epsilon$,
we have the final quadratic programming problem as

$\displaystyle\min_{\beta}\hskip 24.39261pt$
$\displaystyle\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\beta_{i}k(X_{i},X_{j})\beta_{j}+\sum_{i=1}^{n}\beta_{i}A_{i}$
$\displaystyle s.t.\hskip 24.39261pt$
$\displaystyle\tilde{H}_{i}\tilde{Q}_{i}^{t}-H_{i}Q_{i}^{t}-\tilde{H}_{i}/\epsilon\leq\beta_{i}\leq\tilde{H}_{i}\tilde{Q}_{i}^{t}-H_{i}Q_{i}^{t}+H_{i}/\epsilon\;\;\;\forall
i=1,\dots,n$ $\displaystyle\sum_{i=1}^{n}\beta_{i}=0$ (12)

The $t+1$ step solution of $(v_{1},\ldots,v_{n})$ is thus attained from

$\displaystyle\hskip
24.39261pt\sum_{i=1}^{n}v_{i}k(X_{i},X)=-\sum_{i=1}^{n}\beta_{i}k(X_{i},X)$

As for the intercept $v_{0}$, although there exists an analytic solution by
applying the KKT conditions, we solve it using a line search for simplicity.
We iteratively solve the quadratic programming problem until convergence.

##  Appendix B PROOF OF THEOREMS FOR ONE-SIDED PDI

###  B.1 Proof of Theorem 4.1

###### Proof.

First, let $f^{*}\in F\equiv\\{f:\forall
x\in\mathcal{X},P(Y(f(x))>S)=\alpha\\}$. For any $f^{\prime}$ such that
$f^{\prime}\neq f^{*}$, denote $\mathcal{X}_{0}=\\{x:x\in\mathcal{X},\exists
f\in F\;s.t.\;f^{\prime}(x)=f(x)\\}$,
$\mathcal{X}_{1}=\\{x:x\in\mathcal{X},\;f^{\prime}(x)>\max_{f\in F}f(x)\\}$,
and $\mathcal{X}_{2}=\\{x:x\in\mathcal{X},\;f^{\prime}(x)<\min_{f\in
F}f(x)\\}$

Under the ignorability assumption,
$P(Y>S|A=f^{\prime}(x),x)=P(Y(f^{\prime}(x))>S|x)>\alpha$ for
$x\in\mathcal{X}_{1}$, while
$P(Y>S|A=f^{\prime}(x),x)=P(Y(f^{\prime}(x))>S|x)<\alpha$ for
$x\in\mathcal{X}_{2}$.

Then

$\displaystyle R(f^{\prime})-R(f^{*})$ $\displaystyle=$ $\displaystyle
E\Big{[}\frac{1}{P(A\,|\,X)}\bigg{(}\big{(}1-\alpha\big{)}I\big{(}Y>S\big{)}\Big{(}I\big{(}A<f^{\prime}(X)\big{)}-I\big{(}A<f^{*}(X)\big{)}\Big{)}$
$\displaystyle+\alpha I\big{(}Y\leq
S\big{)}\Big{(}I\big{(}A>f^{\prime}(X)\big{)}-I\big{(}A>f^{*}(X)\big{)}\Big{)}\bigg{)}\Big{]}$
$\displaystyle=$ $\displaystyle
E\Bigg{[}\int_{f^{*}(X)}^{f^{\prime}(X)}\Big{(}\big{(}1-\alpha\big{)}P\big{(}Y>S|A=a,X\big{)}$
$\displaystyle\hskip 100.27054pt-\alpha P\big{(}Y\leq
S|A=a,X\big{)}\Big{)}da\bigg{|}X\in\mathcal{X}_{1}\Bigg{]}P(X\in\mathcal{X}_{1})$
(A1)
$\displaystyle+E\Bigg{[}\int_{f^{\prime}(X)}^{f^{*}(X)}\Big{(}-\big{(}1-\alpha\big{)}P\big{(}Y>S|A=a,X\big{)}$
$\displaystyle\hskip 100.27054pt+\alpha P\big{(}Y\leq
S|A=a,X\big{)}\Big{)}da\bigg{|}X\in\mathcal{X}_{2}\Bigg{]}P\big{(}X\in\mathcal{X}_{2}\big{)}$
(A2)
$\displaystyle+E\Bigg{[}\int_{f^{*}(X)}^{f^{\prime}(X)}\Big{(}\big{(}1-\alpha\big{)}\alpha-\alpha\big{(}1-\alpha\big{)}\Big{)}da\bigg{|}X\in\mathcal{X}_{0}\Bigg{]}P\big{(}X\in\mathcal{X}_{0}\big{)}$
(A3)

Here

$\begin{split}(A3)=&0\\\
(A1)=&E\left[\int_{f^{*}(X)}^{f^{\prime}(X)}\Big{(}P\big{(}Y>S|A=a,X\big{)}-\alpha\Big{)}da\bigg{|}X\in\mathcal{X}_{1}\right]P\big{(}X\in\mathcal{X}_{1}\big{)}\\\
=&E\left[\int_{\max_{f\in\mathcal{F}}f(X)}^{f^{\prime}(X)}\Big{(}P\big{(}Y>S|A=a,X\big{)}-\alpha\Big{)}da\bigg{|}X\in\mathcal{X}_{1}\right]P\big{(}X\in\mathcal{X}_{1}\big{)}\\\
\geq&0\\\ >&0\qquad\text{ if }\quad P(X\in\mathcal{X}_{1})\neq 0\\\
\end{split}$

$\begin{split}(A2)=&E\left[\int_{f^{\prime}(X)}^{f^{*}(X)}\Big{(}-P\big{(}Y>S|A=a,X\big{)}+\alpha\Big{)}da\bigg{|}X\in\mathcal{X}_{2}\right]P(X\in\mathcal{X}_{2})\\\
=&E\left[\int_{f^{\prime}(X)}^{\min_{f\in\mathcal{F}}f(X))}\Big{(}-P\big{(}Y>S|A=a,X\big{)}+\alpha\Big{)}da\bigg{|}X\in\mathcal{X}_{2}\right]P(X\in\mathcal{X}_{2})\\\
\geq&0\\\ >&0\qquad\text{ if }\quad P(X\in\mathcal{X}_{2})\neq 0\\\
\end{split}$

Therefore, $R(f^{\prime})-R(f^{*})\geq 0$. And $R(f^{\prime})-R(f^{*})=0$ if
and only if $P(X\in\mathcal{X}_{1})=P(X\in\mathcal{X}_{2})=0$, i.e,
$f^{\prime}\in F$. Hence any possible optimal left bound function
$f_{L,opt}(x)=\arg\min_{f}R(f)$ belongs to $F$. In addition, according to the
assumptions (C1) and (C2) of partial monotonicity for PDI, we have that
$P\left(Y(f_{L,opt}(x))>S\right)=\alpha$ and
$P\left(Y(a^{\prime}(x))>S\right)\geq\alpha$ for $\forall x\in\mathcal{X}$,
where $a^{\prime}(x)$ is an arbitrary measurable function s.t.
$f_{L,opt}(x)\leq a^{\prime}(x)\leq a_{U}$. Theorem 4.1 is proved.

∎

###  B.2 Proof of Theorem 4.2

###### Proof.

For any measurable function $f:\mathcal{X}\rightarrow\mathbb{R}$, we have

$\displaystyle|R(f)-R_{\epsilon_{n}}(f)|$ $\displaystyle=$ $\displaystyle
E\bigg{|}\frac{1}{\epsilon_{n}P(A\,|\,X)}\bigg{(}\big{(}1-\alpha\big{)}I(Y>S)\max(\epsilon_{n}-(f(X)-A)_{+},0)$
$\displaystyle+\alpha I(Y\leq
S)\max(\epsilon_{n}-(A-f(X))_{+},0)\bigg{)}\bigg{|}$ $\displaystyle=$
$\displaystyle
E_{X}\bigg{|}\frac{1}{\epsilon_{n}}\int_{a\in[a_{L},a_{U}]}\bigg{(}\big{(}1-\alpha\big{)}P(Y>S|A=a,X)\max(\epsilon_{n}-(f(X)-a)_{+},0)$
$\displaystyle\qquad\qquad+\alpha P(Y\leq
S|A=a,X)\max(\epsilon_{n}-(a-f(X))_{+},0)\bigg{)}da\;\bigg{|}$
$\displaystyle=$ $\displaystyle
E_{X}\bigg{|}\frac{\big{(}1-\alpha\big{)}}{\epsilon_{n}}\int_{a\in[f(X)-\epsilon_{n},f(X)]}P(Y>S|A=a,X)(\epsilon_{n}-(f(X)-a))da$
$\displaystyle\qquad\qquad+\frac{\alpha}{\epsilon_{n}}\int_{a\in[f(X),f(X)+\epsilon_{n}]}P(Y\leq
S|A=a,X)(\epsilon_{n}-(a-f(X)))da\;\bigg{|}$ $\displaystyle=$ $\displaystyle
E_{X}\bigg{|}\frac{\big{(}1-\alpha\big{)}}{\epsilon_{n}}\int_{z\in[0,\epsilon_{n}]}P(Y>S|A=f(X)-z,X)(\epsilon_{n}-z)dz$
$\displaystyle\qquad\qquad+\frac{\alpha}{\epsilon_{n}}\int_{z\in[0,\epsilon_{n}]}P(Y\leq
S|A=f(X)+z,X)(\epsilon_{n}-z)dz\;\bigg{|}$ $\displaystyle\leq\epsilon_{n}$

Therefore Theorem 4.2 follows. ∎

###  B.3 Proof of Theorem 4.3

###### Proof.

According to Theorem 4.2, we have

$\displaystyle R(\hat{f}_{L})-R(f_{opt})$ $\displaystyle\leq
R_{\epsilon_{n}}(\hat{f}_{L})-R_{\epsilon_{n}}(f_{L,opt})+2C\epsilon_{n}$
$\displaystyle\leq
R_{\epsilon_{n}}(\hat{f}_{L})-R_{\epsilon_{n}}(f_{L}^{*})+2C\epsilon_{n}$
$\displaystyle\leq\bigg{(}\lambda_{n}||\hat{f}_{L}||_{k}^{2}+R_{\epsilon_{n}}(\hat{f}_{L})-R_{\epsilon_{n}}^{*}\bigg{)}+\bigg{(}2C\epsilon_{n}\bigg{)}$
$\displaystyle=(I)+(II)$ (13)

where $f_{L}^{*}$ is the minimizer of $R_{\epsilon_{n}}(\cdot)$ and
$R_{\epsilon_{n}}^{*}=R_{\epsilon_{n}}(f_{L}^{*})$.

To bound (13), we need to first bound $(I)$. As in Chen et al., (2016), we
refer to Theorem 7.23 in Steinwart and Christmann, (2008). In order to use the
oracle inequality in the theorem, there are four conditions to be satisfied.

  * •

(B1) The loss function $L(\cdot)$ has a supremum bound $L(\cdot)\leq B$ for a
constant $B>0$.

  * •

(B2) The loss function $L(\cdot)$ is locally Lipschitz continuous and can be
clipped at a constant $M>0$ such that $\tilde{f}=I(|f|\leq M)f+I(|f|\geq M)M$.

  * •

(B3) The variance bound $E_{P}\big{(}L\circ\tilde{f}-L\circ
f^{*}_{L,P}\big{)}^{2}\leq V\Big{(}E_{P}\big{(}L\circ\tilde{f}-L\circ
f^{*}_{L,P}\big{)}\Big{)}^{v}$ is satisfied for a constant $v\in[0,1]$, $V\geq
B^{2-v}$ and all $f\in\mathcal{H}$.

  * •

(B4) For fixed $n\geq 1$, there exist constants $p\in(0,1)$ and $a\geq B$ such
that the entropy number $E_{D_{X}\sim
P_{X}^{n}}e_{i}(id:\mathcal{H}\rightarrow L_{2}(D_{X}))\leq
ai^{-\frac{1}{2p}}$, $i\leq 1$.

To verify Condition (B1), recall that the relaxed loss function of the lower
bound of the one-sided PDI is

$\displaystyle
L_{\epsilon}(X,A,f_{L}(X))=\frac{1}{P(A\,|\,X)}\Big{\\{}\big{(}1-\alpha\big{)}I(Y>S)\Psi_{\epsilon}({f}_{L}(X),A)+\alpha
I(Y\leq S)\Psi_{\epsilon}(A,{f}_{L}(X))\Big{\\}}$

Assuming the inverse probability ${1}/{P(A\,|\,X)}$ is bounded by a constant
$B$, then $L_{\epsilon}(X,A,f_{L}(X))$ is naturally bounded by $B$ as well.

To verify Condition (B2), similar to (B1), $L_{\epsilon}(X,A,f_{L}(X))$ is
Lipschitz continuous with a Lipschitz constant $B$. Also, $f_{L}$ can be
clipped to have a smaller risk. Suppose that $M$ is an upper bound of the
absolute value of a reasonable range of dose, say,
$M=\max\\{|a_{L}|,|a_{U}|\\}$, and $\tilde{f}_{L}=I(|f_{L}|\leq
M)f_{L}+I(|f_{L}|\geq M)M$. It naturally follows that $R(\tilde{f}_{L})\leq
R(f_{L})$, since any unreasonably large dose recommendation introduces a
larger risk.

Condition (B3) is satisfied with $v=0$ and $V=4B^{2}$, because

$\displaystyle E\big{[}\big{(}L_{\epsilon}\circ\tilde{f}_{L}-L_{\epsilon}\circ
f_{L}^{*}\big{)}^{2}\big{]}\leq
2E\big{[}\big{(}L_{\epsilon}\circ\tilde{f}_{L}\big{)}^{2}+\big{(}L_{\epsilon}\circ
f_{L}^{*}\big{)}^{2}\big{]}\leq 4B^{2}$

The Gaussian Kernel used in Sections 5.1 and 5.2 is one type of benign
kernels. According to Theorem 7.34 of Steinwart and Christmann, (2008), (B4)
is satisfied with the constant
$a=c_{\epsilon,p}\gamma_{n}^{-\frac{(1-p)(1+\zeta)d}{2p}}$, where
$d/(d+\tau)<p<1$ and $\zeta>0$ are two constants.

Since Conditions (B1)-(B4) are satisfied in our case, applying Theorem 7.23
from Steinwart and Christmann, (2008) yields

$\displaystyle(I)$ $\displaystyle\leq
9*\bigg{\\{}\lambda_{n}||f_{L}^{0}||_{k}^{2}+R_{\epsilon_{n}}(f_{L}^{0})-R_{\epsilon_{n}}^{*}\bigg{\\}}+K_{0}\left[\frac{1}{\gamma_{n}^{(1-p)(1+\zeta)d}\lambda_{n}^{p}n}\right]^{\frac{1}{2-p}}+36\sqrt{2}B\Big{(}\frac{\tau}{n}\Big{)}^{\frac{1}{2}}+\frac{15B\tau}{n}$
(14)

with probability $1-3e^{-\tau}$.

In order to bound
$\big{(}\lambda_{n}||f_{L}^{0}||_{k}^{2}+R_{\epsilon_{n}}(f_{L}^{0})-R_{\epsilon_{n}}^{*}\big{)}$,
the approximation error, we refer to Section 2 in Eberts et al., (2013). As
(14) holds for any $f_{L}^{0}\in\mathcal{H}_{\gamma}$, we construct a specific
$f_{L}^{0}$ to facilitate the proof. First, similar to Equation (8) in Eberts
et al., (2013), we define a function $\mathcal{K}(x)=\sum_{j=1}^{r}{r\choose
j}(-1)^{1-j}\frac{1}{j^{d}}(\frac{2}{\gamma^{2}})^{d/2}\mathcal{K}_{j\gamma/\sqrt{2}}(x)$
where $\mathcal{K}_{\gamma}=\exp(-\gamma^{2}||x||_{2}^{2})$ for all $x\in
R^{d}$, and subsequently define $f_{L}^{0}$ via convolution as follows,

$\displaystyle
f_{L}^{0}=\mathcal{K}*f_{L}^{*}=\int_{R^{d}}\mathcal{K}(X-t)f_{L}^{*}(t)dt,\;\;x\in
R^{d}$

If we assume $f^{*}_{L}\in L_{2}(R^{d})\cap L_{\infty}(R^{d})$, then from
Theorem 2.3 in Eberts et al., (2013), it can be shown that
$f_{L}^{0}\in\mathcal{H}_{\gamma}$, where $\mathcal{H}_{\gamma}$ is the RKHS
with the Gaussian kernel. In addition, Theorem 2.2 provides the upper bound
for the access risk of $f_{L}^{0}$, which can be incorporated to yield the
following results.

$\displaystyle\lambda_{n}||f_{L}^{0}||_{k}^{2}+R_{\epsilon_{n}}(f_{L}^{0})-R_{\epsilon_{n}}^{*}$
$\displaystyle=\lambda_{n}||\mathcal{K}*f_{L}^{*}||_{k}^{2}+R_{\epsilon_{n}}(\mathcal{K}*f_{L}^{*})-R_{\epsilon_{n}}^{*}$
$\displaystyle\leq\lambda_{n}(\gamma_{n}\sqrt{\pi})^{-d}(2^{r}-1)^{2}||f_{L}^{*}||^{2}_{L_{2}(R^{d})}+R_{\epsilon_{n}}(\mathcal{K}*f_{L}^{*})-R_{\epsilon_{n}}^{*}$
$\displaystyle\leq\lambda_{n}(\gamma_{n}\sqrt{\pi})^{-d}(2^{r}-1)^{2}||f_{L}^{*}||^{2}_{L_{2}(R^{d})}+B|\mathcal{K}*f_{L}^{*}-f_{L}^{*}|_{L_{1}(P_{X})}$
$\displaystyle\leq\lambda_{n}(\gamma_{n}\sqrt{\pi})^{-d}(2^{r}-1)^{2}||f_{L}^{*}||^{2}_{L_{2}(R^{d})}+BC_{r,1}||g||_{L_{p}(P_{X})}\omega_{r,L_{1}(R^{d})}(f^{*}_{L},\gamma_{n}/2)$
(15)

Given the fact that $f_{opt}\in B_{1,\infty}^{\delta}(R^{d})$, if we further
assume $f_{L}^{*}\in B_{1,\infty}^{\delta}(R^{d})$, i.e.,
$B^{\delta}_{1,\infty}(R^{d})=\\{f\in
L_{\infty}(R^{d}):\sup_{t>0}(t^{-\delta}\omega_{r,L_{1}(R^{d})})(f,t)<\infty\\}$,
then $w_{r,L_{1}}(R^{d})(f_{L}^{*},\gamma_{n}/2)<c_{0}\gamma^{\delta}$, where
$c_{0}$ is a constant. Plugging this into (15), we obtain

$\displaystyle\lambda_{n}||f_{L}^{0}||_{k}^{2}+R_{\epsilon_{n}}(f_{L}^{0})-R_{\epsilon_{n}}^{*}\leq
c_{1}\lambda_{n}\gamma_{n}^{-d}+c_{2}\gamma_{n}^{\delta}$

After combining all the parts together,

$\displaystyle R(\hat{f}_{L})-R(f_{opt})\leq
c_{1}\frac{\lambda_{n}}{\gamma_{n}^{d}}+c_{2}\gamma_{n}^{\delta}+c_{3}\frac{1}{\gamma_{n}^{\frac{(1-p)(1+\zeta)d}{2-p}}\lambda_{n}^{\frac{p}{2-p}}n^{\frac{1}{2-p}}}+c_{4}\frac{\tau^{\frac{1}{2}}}{n^{\frac{1}{2}}}+c_{5}\frac{\tau}{n}+c_{6}\epsilon_{n}$

By properly choosing the constant, i.e.

$\displaystyle\gamma_{n}\propto\bigg{(}\frac{1}{n}\bigg{)}^{\frac{1}{2\delta+d}}\hskip
33.42691pt\lambda_{n}\propto\bigg{(}\frac{1}{n}\bigg{)}^{\frac{\delta+d}{2\delta+d}}\hskip
33.42691pt\epsilon_{n}=\mathcal{O}\bigg{(}n^{-\frac{\delta}{2\delta+d}}\bigg{)}$

we have the convergence rate:

$\displaystyle
R(\hat{f}_{L})-R(f_{opt})=\mathcal{O}\bigg{(}n^{-\frac{\delta}{2\delta+d}}\bigg{)}$

∎

##  Appendix C THEORETICAL RESULTS FOR TWO-SIDED PDI

To ensure the two-sided PDI exists, we need the following assumptions similar
to the one-sided PDI case:

(C3) For any $x$, the boundaries $a_{L}$ and $a_{U}$ for the dose range
satisfy $P\\{Y(a_{L})>S\big{|}x\\}\leq\alpha$ and
$P\big{\\{}Y(a_{U})>S|x\big{\\}}\leq\alpha$.
(C4) There exists a function $a^{*}(\cdot)\in[a_{L},a_{U}]$ such that
$P\\{Y(a^{*}(x))>S|x\\}\geq\alpha$
(C5) For any $x$, there are two numbers $a_{1,x}<a_{2,x}$ with
$a_{1,x}\in[a_{L},a_{U}]$ and $a_{2,x}\in[a_{L},a_{U}]$ such that
$P\\{Y(a_{1,x})>S|x\\}=\alpha$ and $P\\{Y(a_{2,x})>S|x\\}=\alpha$.

######  Corollary C.0.1.

Let $\left[f_{L,opt},f_{U,opt}\right]=\arg\min_{f_{L}\leq
f_{U}}R(f_{L},f_{U})$. Then, under the assumptions (C3)-(C5),
$\left[f_{L,opt},f_{U,opt}\right]$ has the following properties:
(1)
$P\big{\\{}Y\big{(}f_{L,opt}(x)\big{)}>S\,|\,x\big{\\}}=\alpha\;\text{and}\;P\big{\\{}Y\big{(}f_{U,opt}(x)\big{)}>S\,|\,x\big{\\}}=\alpha,\;\forall
x$.
(2) For any measurable function $a(\cdot)$ s.t. $f_{L,opt}(x)\leq a(x)\leq
f_{U,opt}(x)$, $\forall x\in\mathcal{X}$, the potential outcome $Y(a(x))$
satisfies $P(Y(a(x))>S|x)\geq\alpha$.
(3) For any $a_{1}(x)<f_{L,opt}(x)$, $P\big{(}Y(a_{1}(x))>S|x\big{)}<\alpha$,
and for any $a_{2}(x)>f_{U,opt}(x)$, $P\big{(}Y(a_{2}(x))>S|x\big{)}<\alpha$.

###### Proof.

First,

$F_{L}\equiv\\{f:\forall x\in\mathcal{X},P(Y(f(x))>S|x)=\alpha\quad
and\quad\exists a^{*}<f(x),\;P(Y(a^{*}(x))>S|x)\leq\alpha\\}$
$F_{U}\equiv\\{f:\forall x\in\mathcal{X},P(Y(f(x))>S|x)=\alpha\quad
and\quad\exists a^{*}>f(x),\;P(Y(a^{*}(x))>S|x)\leq\alpha\\}$

Let $f^{*}_{L}\in F_{L}$ and $f^{*}_{U}\in F_{U}$, then for arbitrary
$f^{\prime}_{L}$ and $f^{\prime}_{U}$, the sample space of $\mathcal{X}$ can
be decomposed as
$\mathcal{X}=\mathcal{X}_{0}\cup\mathcal{X}_{1}\cup\mathcal{X}_{2}\cup\mathcal{X}_{3}\cup\mathcal{X}_{4}\cup\mathcal{X}_{5}\cup\mathcal{X}_{6}$
where

$\displaystyle\mathcal{X}_{0}=\\{x:x\in\mathcal{X},\quad\exists f_{L}\in
F_{L}\;s.t.\;f^{\prime}_{L}(x)=f_{L}(x)\quad and\quad\exists f_{U}\in
F_{U}\;s.t.\;f^{\prime}_{U}(x)=f_{U}(x)\\}$
$\displaystyle\mathcal{X}_{1}=\\{x:x\in\mathcal{X},\quad\max_{f_{L}\in
F_{L}}f_{L}(x)<f^{\prime}_{L}(x)\leq f^{\prime}_{U}(x)<\min_{f_{U}\in
F_{U}}f_{U}(x)\\}$ $\displaystyle\mathcal{X}_{2}=\\{x:x\in\mathcal{X},\quad
f^{\prime}_{L}(x)<\min_{f_{L}\in F_{L}}f_{L}(x)\leq\max_{f_{U}\in
F_{U}}f_{U}(x)<f^{\prime}_{U}(x)\\}$
$\displaystyle\mathcal{X}_{3}=\\{x:x\in\mathcal{X},\quad
f^{\prime}_{L}(x)\leq\min_{f_{L}\in F_{L}}f_{L}(x)\leq
f^{\prime}_{U}(x)\leq\max_{f_{U}\in F_{U}}f_{U}(x)\\}$
$\displaystyle\mathcal{X}_{4}=\\{x:x\in\mathcal{X},\quad\min_{f_{L}\in
F_{L}}f_{L}(x)\leq f^{\prime}_{L}(x)\leq\max_{f_{U}\in F_{U}}f_{U}(x)\leq
f^{\prime}_{U}(x)\\}$ $\displaystyle\mathcal{X}_{5}=\\{x:x\in\mathcal{X},\quad
f^{\prime}_{L}(x)\leq f^{\prime}_{U}(x)<\min_{f_{L}\in
F_{L}}f_{L}(x)\leq\max_{f_{U}\in F_{U}}f_{U}(x)\\}$
$\displaystyle\mathcal{X}_{6}=\\{x:x\in\mathcal{X},\quad\min_{f_{L}\in
F_{L}}f_{L}(x)\leq\max_{f_{U}\in F_{U}}f_{U}(x)<f^{\prime}_{L}(x)\leq
f^{\prime}_{U}(x)\\}$

Then

$\displaystyle R(f^{\prime})-R(f^{*})=$ $\displaystyle
E\bigg{[}\frac{1}{P(A\,|\,X)}\bigg{(}\big{(}1-\alpha\big{)}I\big{(}Y>S\big{)}\Big{(}I(A\notin[f^{\prime}_{L}(X),f^{\prime}_{U}(X)])-I(A\notin[f^{*}_{L}(X),f^{*}_{U}(X)])\Big{)}$
$\displaystyle+\alpha I\big{(}Y\leq
S\big{)}\Big{(}I(A\in[f^{\prime}_{L}(X),f^{\prime}_{U}(X)])-I(A\in[f^{*}_{L}(X),f^{*}_{U}(X)])\Big{)}\bigg{)}\bigg{]}$
$\displaystyle=$ $\displaystyle
E\left[\dots\big{|}X\in\mathcal{X}_{1}\right]P(X\in\mathcal{X}_{0})$ (M0)
$\displaystyle+E\left[\dots\big{|}X\in\mathcal{X}_{1}\right]P(X\in\mathcal{X}_{1})$
(M1)
$\displaystyle+E\left[\dots\big{|}X\in\mathcal{X}_{2}\right]P(X\in\mathcal{X}_{2})$
(M2)
$\displaystyle+E\left[\dots\big{|}X\in\mathcal{X}_{3}\right]P(X\in\mathcal{X}_{3})$
(M3)
$\displaystyle+E\left[\dots\big{|}X\in\mathcal{X}_{4}\right]P(X\in\mathcal{X}_{4})$
(M4)
$\displaystyle+E\left[\dots\big{|}X\in\mathcal{X}_{5}\right]P(X\in\mathcal{X}_{5})$
(M5)
$\displaystyle+E\left[\dots\big{|}X\in\mathcal{X}_{6}\right]P(X\in\mathcal{X}_{6})$
(M6)

where $(M3)\geq 0$, $(M3)>0$ if $P(X\in\mathcal{X}_{3})>0$, and $(M4)\geq 0$,
$(M4)>0$ if $P(X\in\mathcal{X}_{3})>0$, according to the results of Theorem
4.1. These results can be directly applied since it is easily seen that within
$\mathcal{X}_{3}$ and $\mathcal{X}_{4}$, assumptions for Theorem 4.1 are
satisfied. Obviously $(M0)=0$ and

$\displaystyle(M1)=$ $\displaystyle
E\Bigg{[}\int_{f^{*}_{L}(X)}^{f^{\prime}_{L}(X)}\Big{(}\big{(}1-\alpha\big{)}P(Y>S|A=a,X)-\alpha
P(Y\leq S|A=a,X)\Big{)}da$
$\displaystyle+\int^{f^{*}_{U}(X)}_{f^{\prime}_{U}(X)}\Big{(}\big{(}1-\alpha\big{)}P(Y>S|A=a,X)-\alpha
P(Y\leq
S|A=a,X)\Big{)}da\bigg{|}X\in\mathcal{X}_{1}\Bigg{]}P(X\in\mathcal{X}_{1})$
$\displaystyle=$ $\displaystyle E\bigg{[}\int_{\max_{f_{L}\in
F_{L}}f_{L}(X)}^{f^{\prime}_{L}(X)}\Big{(}P(Y>S|A=a,X)-\alpha\Big{)}da$
$\displaystyle+\int^{\min_{f_{U}\in
F_{U}}f_{U}(X)}_{f^{\prime}_{U}(X)}\Big{(}P(Y>S|A=a,X)-\alpha\Big{)}da\bigg{|}X\in\mathcal{X}_{1}\bigg{]}P(X\in\mathcal{X}_{1})$

Therefore $(M1)\geq 0$ and $>0$ if $P(X\in\mathcal{X}_{1})\neq 0$. Similarly,

$\displaystyle(M2)=$ $\displaystyle
E\bigg{[}\int^{f^{*}_{L}(X)}_{f^{\prime}_{L}(X)}\Big{(}-\big{(}1-\alpha\big{)}P(Y>S|A=a,X)+\alpha
P(Y\leq S|A=a,X)\Big{)}da$
$\displaystyle+\int_{f^{*}_{U}(X)}^{f^{\prime}_{U}(X)}\Big{(}-\big{(}1-\alpha\big{)}P(Y>S|A=a,X)+\alpha
P(Y\leq
S|A=a,X)\Big{)}da\bigg{|}X\in\mathcal{X}_{2}\bigg{]}P(X\in\mathcal{X}_{2})$
$\displaystyle=$ $\displaystyle E\bigg{[}\int^{\min_{f_{L}\in
F_{L}}f_{L}(X)}_{f^{\prime}_{L}(X)}\Big{(}-P(Y>S|A=a,X)+\alpha\Big{)}da$
$\displaystyle+\int_{\max_{f_{U}\in
F_{U}}f_{U}(X)}^{f^{\prime}_{U}(X)}\Big{(}-P(Y>S|A=a,X)+\alpha\Big{)}da\bigg{|}X\in\mathcal{X}_{2}\bigg{]}P(X\in\mathcal{X}_{2})$

Therefore $(M2)\geq 0$ and $>0$ if $P(X\in\mathcal{X}_{2})\neq 0$.

$\displaystyle(M5)=$ $\displaystyle
E\Bigg{[}\int_{f^{\prime}_{L}(X)}^{f^{\prime}_{U}(X)}\Big{(}-\big{(}1-\alpha\big{)}P(Y>S|A=a,X)+\alpha
P(Y\leq S|A=a,X)\Big{)}da$
$\displaystyle+\int^{f^{*}_{L}(X)}_{f^{*}_{U}(X)}\Big{(}\big{(}1-\alpha\big{)}P(Y>S|A=a,X)-\alpha
P(Y\leq
S|A=a,X)\Big{)}da\bigg{|}X\in\mathcal{X}_{5}\Bigg{]}P(X\in\mathcal{X}_{5})$
$\displaystyle=$ $\displaystyle
E\bigg{[}\int_{f^{\prime}_{L}(X)}^{f^{\prime}_{U}(X)}\Big{(}-P(Y>S|A=a,X)+\alpha\Big{)}da$
$\displaystyle+\int^{f^{*}_{L}(X)}_{f^{*}_{U}(X)}\Big{(}P(Y>S|A=a,X)-\alpha\Big{)}da\bigg{|}X\in\mathcal{X}_{5}\bigg{]}P(X\in\mathcal{X}_{5})$

Therefore $(M5)\geq 0$ and $>0$ if $P(X\in\mathcal{X}_{2})\neq 0$.

$\displaystyle(M6)=$ $\displaystyle
E\Bigg{[}\int_{f^{\prime}_{L}(X)}^{f^{\prime}_{U}(X)}\Big{(}-\big{(}1-\alpha\big{)}P(Y>S|A=a,X)+\alpha
P(Y\leq S|A=a,X)\Big{)}da$
$\displaystyle+\int^{f^{*}_{L}(X)}_{f^{*}_{U}(X)}\Big{(}\big{(}1-\alpha\big{)}P(Y>S|A=a,X)$
$\displaystyle\hskip 100.27054pt-\alpha P(Y\leq
S|A=a,X)\Big{)}da\bigg{|}X\in\mathcal{X}_{6}\Bigg{]}P(X\in\mathcal{X}_{6})$
$\displaystyle=$ $\displaystyle
E\bigg{[}\int_{f^{\prime}_{L}(X)}^{f^{\prime}_{U}(X)}\Big{(}-P(Y>S|A=a,X)+\alpha\Big{)}da$
$\displaystyle+\int^{f^{*}_{L}(X)}_{f^{*}_{U}(X)}\Big{(}P(Y>S|A=a,X)-\alpha\Big{)}da\bigg{|}X\in\mathcal{X}_{6}\bigg{]}P(X\in\mathcal{X}_{6})$

Therefore $(M6)\geq 0$ and $>0$ if $P(X\in\mathcal{X}_{6})\neq 0$.
Consequently, $R(f^{\prime}_{U},f^{\prime}_{U})-R(f^{*}_{L},f^{*}_{U})\geq 0$,
and $R(f^{\prime}_{U},f^{\prime}_{U})-R(f^{*}_{L},f^{*}_{U})=0$ if and only if
$P(X\notin\mathcal{X}_{0})=0$. Hence for any optimal two-sided bound functions
$[f_{L,opt},f_{U,opt}]=\arg\min_{f_{L}\leq f_{U}}R(f_{L},f_{U})$, we have
$f_{L,opt}\in F_{L}$ and $f_{L,opt}\in F_{U}$. According to the assumptions
(C3), (C4) and (C5), we have that
$P\left(Y\big{(}f_{L,opt}(x)\big{)}>S|x\right)=P\left(Y\big{(}f_{U,opt}(x)\big{)}>S|x\right)=\alpha$
and $P\left(Y(a^{\prime}_{x})>S|x\right)\geq\alpha$ for $\forall
x\in\mathcal{X}$, where $a^{\prime}(x)$ is an arbitrary measurable function
s.t. $f_{L,opt}(x)\leq a^{\prime}(x)\leq f_{U,opt}(x)$. ∎

######  Corollary C.0.2.

For any interval of two measurable functions $\left[f_{L}(x),f_{U}(x)\right]$,
where $f_{L}(x)\leq f_{U}(x),\forall x\in\mathcal{X}$, we have
$|R(f_{L},f_{U})-R_{\epsilon_{n}}(f_{L},f_{U})|\leq C\epsilon_{n}$, where $C$
is a constant.

###### Proof.

For any measurable function $f_{L}:\mathcal{X}\rightarrow\mathbb{R}$,
$f_{U}:\mathcal{X}\rightarrow\mathbb{R}$ and $f_{L}(X)\leq f_{U}(X)\forall
x\in\mathcal{X}$, we have

$\displaystyle|$ $\displaystyle R(f)-R_{\epsilon_{n}}(f)|$
$\displaystyle=E_{X}\bigg{|}\frac{\big{(}1-\alpha\big{)}}{\epsilon_{n}}\int_{a\in[f_{L}(X)-\epsilon_{n},f_{L}(X)]}P(Y>S|A=a,X)(\epsilon_{n}-(f_{L}(X)-a))da$
$\displaystyle\qquad\qquad+\frac{\alpha}{\epsilon_{n}}\int_{a\in[f_{L}(X),f_{L}(X)+\epsilon_{n}]}P(Y\leq
S|A=a,X)(\epsilon_{n}-(a-f_{L}(X)))da$
$\displaystyle\qquad\qquad+\frac{\alpha}{\epsilon_{n}}\int_{a\in[f_{U}(X)-\epsilon_{n},f_{U}(X)]}P(Y\leq
S|A=a,X)(\epsilon_{n}-(f_{U}(X)-a))da$
$\displaystyle\qquad\qquad+\frac{\big{(}1-\alpha\big{)}}{\epsilon_{n}}\int_{a\in[f_{U}(X),f(X)_{R}+\epsilon_{n}]}P(Y>S|A=a,X)(\epsilon_{n}-(a-f_{U}(X)))da\;\bigg{|}$
$\displaystyle=E_{X}\bigg{|}\frac{\big{(}1-\alpha\big{)}}{\epsilon_{n}}\int_{z\in[0,\epsilon_{n}]}P(Y>S|A=f_{L}(X)-z,X)(\epsilon_{n}-z)dz$
$\displaystyle\qquad\qquad+\frac{\alpha}{\epsilon_{n}}\int_{z\in[0,\epsilon_{n}]}P(Y\leq
S|A=f_{L}(X)+z,X)(\epsilon_{n}-z)dz$
$\displaystyle\qquad\qquad+\frac{\alpha}{\epsilon_{n}}\int_{z\in[0,\epsilon_{n}]}P(Y\leq
S|A=f_{U}(X)-z,X)(\epsilon_{n}-z)dz$
$\displaystyle\qquad\qquad+\frac{\big{(}1-\alpha\big{)}}{\epsilon_{n}}\int_{z\in[0,\epsilon_{n}]}P(Y>S|A=f_{U}(X)+z,X)(\epsilon_{n}-z)dz\;\bigg{|}$
$\displaystyle=2\epsilon_{n}.$

∎

We provide a convergence rate by borrowing the results from Theorem 4.3, and
replace the convergence of the two-sided PDI estimator with two estimators of
one-sided PDIs. In order to do so, extra assumptions have to be made. It would
be an interesting future work if better rates under weaker assumptions can be
achieved.

######  Corollary C.0.3.

Assume that $f_{L,opt}\in B^{\delta}_{1,\infty}(\mathbb{R}^{d})$ and
$f_{U,opt}\in B^{\delta}_{1,\infty}(\mathbb{R}^{d})$ where
$B^{\delta}_{1,\infty}(\mathbb{R}^{d})=\\{f\in
L_{\infty}((\mathbb{R}^{d})):\sup_{t>0}(t^{-\delta}w_{(r,L_{1})}(f,t))<\infty\\}$
and $w$ is the modulus of continuity. In addition, assume that there exists a
measurable function $f_{M,opt}(x)$, which is known or can be estimated
consistently, such that $f_{L,opt}(x)<f_{M}(x)<f_{U,opt}(x)$ for all
$x\in\mathcal{X}$. Then for any $\eta>0$, $d/(d+\tau)<p<1$, $\tau>0$,
$0<c<1$,and properly chosen $\gamma_{n}$, $\lambda_{n}$ and $\epsilon_{n}$, we
have

$R(\hat{f}_{L,n},\hat{f}_{U,n})-R(f_{L,opt},f_{U,opt})=\mathcal{O}_{p}\Bigg{(}n^{-\frac{\delta}{2\delta+d}}\Bigg{)}$
(16)

with probability $c\big{(}1-3e^{-\tau}\big{)}^{2}$. Here $d$ is the dimension
of $\mathcal{X}$.

###### Proof.

The proof of Corollary C.0.3 consists of two steps. In the first step, we
assume that there exists a function or a consistent estimator of a function
that separates the dataset into two subsets which we call as left and right
subsets, in each of which Conditions (C1) and (C2) are satisfied. In the
second step, the integrated excessive risk of the two-sided PDI estimator is
bounded by the sum of excessive risks of two independent one-sided PDI
estimators.

Suppose there is a function $f_{M,opt}(x)$, such that
$f_{L,opt}(x)<f_{M,opt}(x)<f_{U,opt}(x)$ for all $x\in\mathcal{X}$. Such a
function is typically unknown in practice, hence we consider the case when
$f_{M,opt}(x)$ is estimated from data. One choice of $f_{M,opt}$ is
$f_{\max,opt}=\arg\max_{f}E\big{[}Y|A=f(X),X\big{]}$, the optimal individual
dose rule (IDR) in Chen et al., (2016). Theoretical properties of the
estimator $\hat{f}_{\max}$ have been established in Chen et al., (2016) and
also discussed in Luedtke and van der Laan, (2016). Despite the fact that it
is suggested in Luedtke and van der Laan, (2016) that the O-learning based
estimator might convergence in nearly $\mathcal{O}(n^{-1/2})$, the rate
established there was based on the risk (or value function) instead of the
rule $f_{\max,opt}$ itself. To facilitate proving the properties of the two-
sided PDI, we make the following further assumption for the estimated optimal
dose function $\hat{f}_{\max}$.

$\displaystyle f_{L,opt}(x)<\hat{f}_{\max}(x)<f_{U,opt}(x),\quad\forall
x\in\mathcal{X}$ (17)

In the second step, two subsets are defined where the lower bound function and
the upper bound function can be estimated separately using the corresponding
subset, assuming there is a $f_{M,opt}$, known or estimated, such that
$f_{L,opt}(x)<f_{M,opt}(x)<f_{U,opt}(x)$ for all $x\in\mathcal{X}$. Denote the
original dataset as $\mathcal{D}=\\{\big{(}X_{i},A_{i},Y_{i}\big{)},\;1\leq
i\leq n\\}$ and then we can subsequently define the left subset as
$\mathcal{D_{L}}=\\{\big{(}X_{i},A_{i},Y_{i}\big{)},\;\forall
i\;s.t.\;A_{i}\leq f_{M,opt}(X_{i})\\}$ and the right subset as
$\mathcal{D_{U}}=\\{\big{(}X_{i},A_{i},Y_{i}\big{)},\;\forall
i\;s.t.\;A_{i}>f_{M,opt}(X_{i})\\}$. These two subsets are independent as
there is no shared observations.

Now define $\hat{f}_{L}^{\prime}=\min_{f_{L}}\hat{R}_{\epsilon}(f_{L})$ on
$\mathcal{D_{L}}$ and
$\hat{f}_{U}^{\prime}=\min_{f_{U}}\hat{R}_{\epsilon}(f_{U})$ on
$\mathcal{D_{U}}$, while
$\big{\\{}\hat{f}_{L},\hat{f}_{U}\big{\\}}=\min_{f_{L},f_{U}}\hat{R}_{\epsilon}(f_{L},f_{U})$
on the dataset of $\mathcal{D}$. Similarly, define
$f_{L,*}^{\prime}=\min_{f_{L}}R_{\epsilon}(f_{L})$ on the (sub)population of
$\mathcal{D_{L}}$ and $f_{U,*}^{\prime}=\min_{f_{U}}R_{\epsilon}(f_{U})$ on
the (sub)population of $\mathcal{D_{U}}$, while
$\big{\\{}f_{L,*},f_{U,*}\big{\\}}=\min_{f_{L},f_{U}}R_{\epsilon}(f_{L},f_{U})$
on the original population. To simplify the problem, we assume that
$f_{L,*}(x)-f_{U,*}(x)\leq-2\epsilon$, since this inequality will hold when
$\epsilon\rightarrow 0$ as $n\rightarrow 0$. By the definitions of
$\mathcal{D_{L}}$, $\mathcal{D_{U}}$ and $\mathcal{D}$, we have
$\big{\\{}f_{L,*}^{\prime},f_{U,*}^{\prime}\big{\\}}=\big{\\{}f_{L,*},f_{U,*}\big{\\}}$,
i.e., the population estimates of the lower and upper bounds estimated jointly
are the same as the population estimates of the lower and upper bounds
estimated individually while assuming the other part is known. In addition, we
make the following assumption of the uniform convergence of the empirical risk
function.

$\displaystyle\sup_{f_{L},f_{U}\in\mathcal{H}}\big{|}\hat{R}_{\epsilon}(f_{L},f_{U})-R_{\epsilon}(f_{L},f_{U})\big{|}\leq\kappa\sqrt{\frac{\log(n)}{n}}\;\text{with
probability }c$ (18)

where $\kappa$ is a constant related to $c$ and the complexity of the RKHS.

We will show how the loss function of the two-sided PDI can be decomposed as
the sum of two one-sided PDI losses plus a quantity not related to the dose.
Recall that the loss function for the lower bound of the one-sided PDI is

$\displaystyle
L_{\epsilon}(X,A,f_{L}(X))=\frac{1}{P(A\,|\,X)}\Big{\\{}\big{(}1-\alpha\big{)}I(Y>S))\Psi_{\epsilon}({f}_{L}(X),A)+\alpha
I(Y\leq S))\Psi_{\epsilon}(A,{f}_{L}(X))\Big{\\}}$

The loss function for the upper bound of the one-sided PDI is

$\displaystyle
L_{\epsilon}(X,A,f_{U}(X))=\frac{1}{P(A\,|\,X)}\Big{\\{}\big{(}1-\alpha\big{)}I(Y>S))\Psi_{\epsilon}(A,{f}_{U}(X))+\alpha
I(Y\leq S))\Psi_{\epsilon}({f}_{U}(X),A)\Big{\\}}.$

And the loss function for the two-sided PDI is

$\displaystyle
L_{\epsilon}(X,A,f_{L}(X),f_{U}(X))=\frac{1}{P(A\,|\,X)}\Big{\\{}\big{(}1-\alpha\big{)}I(Y>S))\Psi^{out}_{\epsilon}({f}_{L}(X),A,f_{U}(X))$
$\displaystyle\hskip 146.34824pt+\alpha I(Y\leq
S))\Psi^{in}_{\epsilon}(f_{L}(X),A,{f}_{U}(X))\Big{\\}}.$

It follows that the difference between the two-sided loss and the sum of the
one-sided losses is a constant which does not contain the bound functions,
under the assumption that $f_{L,*}(x)-f_{U,*}(x)\leq-2\epsilon$.

$\displaystyle
L_{\epsilon}(X,A,f_{L}(X),f_{U}(X))-L_{\epsilon}(X,A,f_{L}(X))-L_{\epsilon}(X,A,f_{U}(X))$
$\displaystyle=\frac{1}{P(A\,|\,X)}\Big{\\{}\big{(}1-\alpha\big{)}I(Y>S)\Psi_{\epsilon}(A,f_{U}(X))+\alpha
I(Y\leq
S)\big{(}\Psi_{\epsilon}(f_{U}(X),A)-1\big{)}\Big{\\}}-L_{\epsilon}(X,A,f_{U}(X))$
$\displaystyle=\frac{-1}{P(A\,|\,X)}\Big{\\{}\alpha I(Y\leq S)\Big{\\}}$

Combined with Corollary C.0.2, we have

$\displaystyle R(\hat{f}_{L},$
$\displaystyle\hat{f}_{U})-R(f_{L,opt},f_{U,opt})\leq
R_{\epsilon_{n}}(\hat{f}_{L},\hat{f}_{U})-R_{\epsilon_{n}}(f_{L,opt},f_{U,opt})+2C\epsilon_{n}$
$\displaystyle\leq
R_{\epsilon_{n}}(\hat{f}_{L},\hat{f}_{U})-R_{\epsilon_{n}}(f_{L,*},f_{U,*})+2C\epsilon_{n}$
$\displaystyle\leq
R_{\epsilon_{n}}(\hat{f}_{L}^{\prime},\hat{f}_{U}^{\prime})-R_{\epsilon_{n}}(f_{L,*},f_{U,*})+2C\epsilon_{n}+2\kappa\sqrt{\frac{\log(n)}{n}}$
$\displaystyle=R_{\epsilon_{n}}(\hat{f}_{L}^{\prime})-R_{\epsilon_{n}}(f_{L,*}^{\prime})+R_{\epsilon_{n}}(\hat{f}_{U}^{\prime})-R_{\epsilon_{n}}(f_{U,*}^{\prime})+2C\epsilon_{n}+2\kappa\sqrt{\frac{\log(n)}{n}}$
(19)

with probability $c$ as $n\rightarrow\infty$. The probability $c$ comes from
the inequality (18). The first inequality of (19) is guaranteed by Theorem
C.0.1. The second inequality is from the definition of $f_{L,*}$ and
$f_{U,*}$. The third inequality is from the assumption (18). The coefficient
is $2\times\kappa$ as the uniform bound has to be applied twice because the
definitions of $\hat{f}_{L}$ and $\hat{f}_{U}$ only guarantee
$\hat{R}_{\epsilon_{n}}(\hat{f}_{L},\hat{f}_{U})<\hat{R}_{\epsilon_{n}}(\hat{f}_{L}^{\prime},\hat{f}_{U}^{\prime})$
but not
$R_{\epsilon_{n}}(\hat{f}_{L},\hat{f}_{U})<R_{\epsilon_{n}}(\hat{f}_{L}^{\prime},\hat{f}_{U}^{\prime})$.
The last equality is from the decomposition of the two-sided loss and the
equivalence of $\big{\\{}f_{L,*}^{\prime},f_{U,*}^{\prime}\big{\\}}$ and
$\big{\\{}f_{L,*},f_{U,*}\big{\\}}$.

The theorem follows by applying (14) twice, because $\hat{f}_{L}^{\prime}$ and
$\hat{f}_{U}^{\prime}$ are independently estimated on two non-overlapping
datasets $\mathcal{D_{L}}$ and $\mathcal{D_{U}}$. That is,
$R_{\epsilon_{n}}(\hat{f}_{L}^{\prime})-R_{\epsilon_{n}}(f_{L,*}^{\prime})$ is
independent with
$R_{\epsilon_{n}}(\hat{f}_{U}^{\prime})-R_{\epsilon_{n}}(f_{U,*}^{\prime})$.
The uniform bound (18) and Theorem 4.3 hold simultaneously with a probability
of at least $c(1-3e^{-\tau})^{2}$. Here $\epsilon_{n}$ only shows up in the
term $2C\epsilon_{n}$. Hence by choosing
$\epsilon_{n}=\mathcal{O}\big{(}n^{-\frac{\delta}{2\delta+d}}\big{)}$, it will
not impact the convergence rate. The term $2\kappa\sqrt{\frac{\log(n)}{n}}$
corresponds to the convergence rate of $\mathcal{O}\big{(}\sqrt{\frac{\log
n}{n}}\big{)}$. Hence, (19) is dominated either by the first two risk terms or
by $2\kappa\sqrt{\frac{\log(n)}{n}}$. ∎

##  Appendix D ADDITIONAL SIMULATIONS WITH NO CONFOUNDING COVARIATES

In this section, we consider scenarios similar to those in the main paper
except that there is no confounding. In particular, the treatment $A$ of the
training data follows a uniform distribution $\text{Unif}(-2,2)$ (i.e.,
independent of $X$) while other data generating mechanisms are the same as in
the main paper. We repeat each setting for $100$ times and report results by
evaluating the empirical risks on corresponding independent testing data sets
with sample sizes of $10000$. Similar to the settings with confounding, our
methods perform better than the competitive methods.

Table 2: Empirical risk $R(\hat{f}_{L})$ (with SD) under settings without confounding covariates | n | d | L-O-Learning | G-O-Learning | Logistic | SVM | RF
---|---|---|---|---|---|---|---
Scenario 1 | 200 | 10 | 0.137 (0.008) | 0.132 (0.011) | 0.154 (0.024) | 0.139 (0.011) | 0.136 (0.008)
$\sigma^{2}=9$ | 200 | 50 | 0.157 (0.014) | 0.133(0.011) | 0.185(0.023) | 0.181 (0.043) | 0.135 (0.012)
| 400 | 10 | 0.128 (0.006) | 0.126 (0.006) | 0.132(0.011) | 0.129 (0.006) | 0.130 (0.009)
| 400 | 50 | 0.154 (0.006) | 0.128 (0.007) | 0.174 (0.022) | 0.142 (0.008) | 0.129 (0.008)
Scenario 2 | 200 | 10 | 0.136 (0.008) |  0.130 (0.011) | 0.153 (0.022) | 0.138 (0.011) | 0.133 (0.012)
$\sigma^{2}=9$ | 200 | 50 | 0.159 (0.013) | 0.134 (0.010) | 0.194 (0.024) | 0.226 (0.037) | 0.139 (0.017)
| 400 | 10 | 0.129 (0.005) | 0.123 (0.006) | 0.132 (0.012) | 0.127 (0.005) | 0.124 (0.005)
| 400 | 50 | 0.156 (0.007) | 0.126 (0.008) | 0.180 (0.026) | 0.144 (0.008) | 0.127 (0.004)

##  Appendix E MORE BACKGROUND ABOUT THE HBA1C CONTROL STUDY

Hemoglobin is an iron-containing oxygen which transports protein in the red
blood cells. HbA1c, the most abundant minor hemoglobin in the human body, is
formed when glucose accumulates in red blood cells and binds to the
hemoglobin. This process occurs slowly and continuously over the lifespan of
red blood cells, which is 120 days on average. This makes HbA1c an ideal
biomarker of long-term glycemic control. Patients who are susceptible to high
blood glucose are usually recommended to have their HbA1c levels regularly
measured and recorded.

HbA1c test, since became commercially available in 1978, has been one of the
standard tools for monitoring diabetes progression. American Diabetes
Association (ADA) recommended using A1c measurement in 1988. The Diabetes
Control and Complications Trial (DCCT) demonstrated its importance as a
predictor of diabetes-related outcomes, and the ADA started recommending
specific A1c targets in 1994 (Little et al.,, 2011). In 2010, ADA added
HbA1c$\geq 6.5\%$ as a criterion for diabetes diagnosis. The prevalence of
HbA1c test can be partially explained by its clinic convenience. As HbA1c is
unaffected by acute perturbations in glucose levels, there is no need for
fasting or timed samples. For people without diabetes, the normal range for
the HbA1c level is between 4% and 5.6%. HbA1c levels between 5.7% and 6.4%
imply a higher chance of getting diabetes. Levels of 6.5% or higher indicate
the presence of diabetes.

Despite the popularity of HbA1c tests, HbA1c may not reflect the real
progression of diabetes given various conditions that patients are predisposed
to. Interfering factors may lead to false results of HbA1c tests (Radin,,
2014). In some cases, the direction of bias is predictable while in other
cases not. For instance, HbA1c can be falsely elevated by iron deficiency
anemias, vitamin B-12 anemias, folate deficiency anemias, asplenia and other
conditions associated with decreased red cell turnover (Nitin,, 2010). HbA1c
measures may be lower for patients with conditions that shorten the life of
the red blood cell or increase its turnover rate (Freedman et al.,, 2010),
including acute and chronic blood loss, hemolytic anemia, splenomegaly and
end-stage renal disease. While in some other cases, the direction of bias is
not always easily predictable. For example, the complex interplay of glycemic
control and other treatments, such as erythropoietin therapy and the treatment
of uremia may influence HbA1c level in a more case-by-case way (Radin,, 2014).

Due to the susceptibility of HbA1c measures to various medical and
physiological factors, HbA1c measures reflect glycemia differently for
different patients. Therefore, the recommendation of a single HbA1c control
level for all the patients may not always be desirable. A lower glycemic level
is not always better, and over-controlling glycemia can lead to faint or
compromised life quality. Especially for patients who are older, with severe
complications, or in the phase of early post-surgery rehabilitation, HbA1c is
not supposed to be controlled for them as low as for younger and healthier
patients.

As a response to the progress in medical research, the idea of glucose control
has been evolving over time. In the past, an HbA1c of 7% was considered the
golden rule of health for everyone. In recent years, however, the importance
of a patient-centered approach to managing HbA1c levels has been recognized,
which may better correspond to the patient’s needs for diabetes management and
their personal conditions and preferences. To fill this gap, we introduce the
HbA1c dataset to study the relationship between HbA1c control and the events
of diabetic patients. Specifically, the sample is a subset of a type II
diabetes dataset which is drawn from a large Midwestern multi-specialty
physician and patient group. Electronic Health Record (EHR) datasets are
merged with Medicare claim data to obtain laboratory records, information on
primary care visits, and medication history. All patients have been tracked
for a specific time period from the first quarter of 2003 to the fourth
quarter of 2011. The full observational records are not available for the
majority of patients.

Patients identified with diabetes must have at least 1 inpatient or SNF
Medicare claim or at least two carrier claims for diabetes, which should be
less than 2 years apart. The diagnosis is defined by one or more of the
following categories: the ICD-9 codes, Diabetes mellitus (250.xx),
Polyneuropathy in diabetes (357.2), Diabetes retinopathy (362.0x), and
Diabetic cataract (366.41). There are several types of patients excluded from
the dataset, including those who are without continuous coverage with Medicare
Parts A & B for a baseline year or for at least one subsequent quarter,
without Medicare railroad benefits or were not enrolled in a Medicare HMO.
Each patient needs to have at least 5 quarter time periods in one of the 3
insurers from 2003 to 2011 including 4 baseline quarters and a measurement
period of at least 1 subsequent quarter. Also, the patients should have A1c
records over the 5 quarter period which are available in the health system. In
the sample of analysis, 51.3% of patients are female and 91.8% of the patients
are white. The ages of participants at the cohort start time range from 18 to
102 years old, with an average of 65.7 and a standard deviation of 13.7.

We regard the HbA1c measurement as the treatment dose in this study.
Observations are centered around the availability of HbA1c measurements. Each
observation is associated with a unique HbA1c measurement from an individual
patient following the baseline period. The outcome variable is the number of
negative events, which include hospitalizations and Emergency Department (ED)
visits that cover 90 days after the HbA1c measurement. In addition, we also
include a list of short-term covariates that we average within 90 days prior
to the HbA1c measurement as well as a list of long-term covariates that are
averaged from up to 900 days prior to the HbA1c measurement. After removing
subjects with missing baseline covariates, $8126$ patients remain. An
illustration of the data structure can be found in Figure 3. Each observation
is generated from an outcome period (90 days after the HbA1c measure), a
short-term covariate period (90 days before the HbA1c measure) and a long-term
covariate period (including the baseline information).

Figure 3: Illustration of Data Structure

The covariates in the data set fall in the following categories.

  1. 1.

Demographics: age, gender, race

  2. 2.

Long-term measurements: height, the indicator of disability at baseline,
length of enrollment, average of HbA1c in the past 91-900 days, the indicator
of Medicaid, number of measurements recorded in past the 91-900 days,
Sulfonylureas intake, Insulin intake, and other medications intake

  3. 3.

Short-term measurements: number of HbA1c measurements in the past 1-90 days;
average of HbA1c in the past 1-90 days, number of hospitalization, number of
ED-visits; low-density lipoprotein cholesterol, systolic blood pressure,
diastolic blood pressure, BMI, weight, cardiovascular disease, ischemic heart
disease, congestive heart failure, hepatocellular carcinoma, hypoglycemia and
injury, infections; and 18 other comorbidities.