# Multi-stage optimal dynamic treatment regimes for survival outcomes with
dependent censoring

Hunyong Cho, Shannon T. Holloway, and Michael R. Kosorok Department of
Biostatistics, University of North Carolina at Chapel Hill,
<EMAIL_ADDRESS>of Statistics, North Carolina State University
Department of Biostatistics, Department of Statistics and Operations Research,
University of North Carolina at Chapel Hill

###### Abstract

We propose a reinforcement learning method for estimating an optimal dynamic
treatment regime for survival outcomes with dependent censoring. The estimator
allows the treatment decision times to be dependent on the failure time and
conditionally independent of censoring, supports a flexible number of
treatment arms and treatment stages, and can maximize either the mean survival
time or the survival probability at a certain time point. The estimator is
constructed using generalized random survival forests, and its consistency is
shown using empirical process theory. Simulations and leukemia data analysis
results suggest that the new estimator brings higher expected outcomes than
existing methods in various settings. An R package dtrSurv is available on
CRAN.

Keywords: Precision medicine; Dynamic treatment regime; Survival analysis;
Reinforcement learning; Random forest; Conditionally independent censoring;
Empirical process.

##  1 Introduction

Multi-stage treatments are becoming more prevalent in medicine. Chronic
conditions such as cancer, auto-immune disease, HIV, and heart disease, often
require multiple treatments over a long period. For instance, cancer patients
may receive multiple rounds of chemo- and/or biological therapies (Habermann
et al., 2006; Huang et al., 2020). Auto-immune disease patients might receive
additional treatments in response to relapses (Edwards & Cambridge, 2006;
Hogan & Radhakrishnan, 2013). Traditionally such chronic diseases have been
treated following a one-size-fits-most approach by which the standard of care
is selected as the treatment or treatment plan that will likely benefit most
patients with a similar condition, i.e., the treatment that is optimal only in
the sense of the marginal effects.

However in recent years, a more personalized approach to the treatment of
chronic diseases has become of interest. Under this approach patients benefit
from treatment rules that incorporate the patient heterogeneity of treatment
effects and thereby treatment decisions are based on the individual
characteristics of each patient. Dynamic treatment regimes formalize this
data-driven approach to treatment by optimizing the overall outcome of
interest and providing treatment suggestions that are based on the patient’s
information available at each point of the decision making process (Murphy,
2003; Kosorok & Moodie, 2015; Kosorok & Laber, 2019; Tsiatis et al., 2019). In
part because they utilize all available information, such optimized dynamic
treatment regimes are often found to be more beneficial than a traditional set
of fixed treatment plans(Kidwell, 2015). Further, because optimal dynamic
treatment regimes focus on the overall or long-term outcome, they are more
advantageous than applying multiple single-stage treatment rules that only
optimize the stage-level outcomes.

Finding an optimal dynamic treatment regime is a reinforcement learning
problem (Sutton & Barto, 2018), wherein an action by an agent, or a physician,
changes both the immediate reward and the environment of the next stage and
the learner, or the statistician, searches for the rule of actions that brings
the most overall reward. Like many statistical learning methods, reinforcement
learning problems are often solved using the Q-learning algorithm (Watkins,
1989; Zhu et al., 2015; Zhao et al., 2011; Zhang et al., 2017; Qian & Murphy,
2011), which defines the value as the expected cumulative sum of discounted
rewards and optimizes the set of rules using backward recursion.

In precision medicine applications, the rewards are often defined in the
context of survival outcomes, such as mean survival time or survival
probabilities, which pose unique challenges. For survival outcomes, the
failure time is often censored, and any approach should take into account the
missing information. Several methods have been developed for this purpose, for
example, redistribution of the failure probability as in the Kaplan-Meier
Efron (1967) and Robertson & Uppuluri (1984) or inverse probability weighting
(Robins et al., 1994; Robins & Rotnitzky, 1992; Wahed & Tsiatis, 2006;
Orellana et al., 2010). Further, the treatment timing may not be fixed and may
even depend on a patient’s status or previous treatments having association
with the failure time. In such scenarios, the three events at each stage are
failure, proceeding to the next treatment, and censoring, which can be
restated in a competing risks framework, in which the three dependent events
compete and only the earliest event is observed. However, the full joint
distribution of the events is not identifiable without knowledge of their
dependence structure (Tsiatis, 1975). Thus, the dependency between the events
needs to be carefully considered. For example, a cancer patient whose physical
condition has improved has a longer expected survival time and at the same
time may want to take another round of chemotherapy earlier than expected. If
this is the case, the survival time and the treatment time have a negative
association.

There are two additional challenges that arise in the analysis of survival
outcomes that must be addressed. First, the total number of treatments
received may vary among patients. For example, some patients may experience
the failure or be censored before taking all planned treatments. And finally,
there could be different optimization criteria for survival times. For some
patients and clinicians, the average survival time is of primary interest,
while the six-month survival probability may be of greater concern for others.
A method that can address different types of endpoints would be useful in
practice.

Several methods regarding the estimation of dynamic treatment regimes for
survival outcomes have been proposed in the literature. Goldberg & Kosorok
(2012) first developed a dynamic treatment regime estimator for survival
outcomes using a Q-learning framework. To address the problems introduced by a
different number of treatments as well as by censoring, the authors proposed
modifying the survival data so that each observation has the same number of
treatment stages without missing values and showed that the value of a policy
in the original problem can be represented as an expectation of the modified
data. After the data modification, the problem was put into a standard
Q-learning framework using inverse probability weighting to account for
censoring. However, censoring was assumed to be completely independent of all
covariates and event times. The estimated Q-function and the resulting
decision rules would be thus subject to bias for data with dependent
censoring.

Huang et al. (2014) addressed a similar problem using backward recursion.
Their motivating problem was a recurrent disease clinical trial where patients
receive an initial treatment followed by a salvage treatment if patients
experience either treatment resistance or relapse. Simoneau et al. (2019)
extended the dynamic treatment regime estimator for continuous outcomes via
weighted least squares (Wallace & Moodie, 2015) into the censored time-to-
event data setting. Both of these methods use the accelerated failure time
model for the failure time distribution, which carries a risk of model
misspecification. Thus, a more flexible model is needed to allow for more
relaxed distributional assumptions. Wahed & Thall (2013) proposed a dynamic
treatment regime estimator using full specification of the likelihood and Xu
et al. (2016) developed a Bayesian alternative, where the disease progression
dynamic was modeled using a dependent Dirichlet process prior with Gaussian
process measure. While the latter approach provides a more flexible modeling
framework than the former, the class of regimes for both methods is refined to
a fixed number of pre-defined treatment sequences.

Each of the existing methods described above has its limitations. First, not
all methods allow for a flexible number of treatment stages and/or treatment
levels. For example, the methods in Wahed & Thall (2013) and Huang et al.
(2014) were designed for a two-stage decision problem, and the Simoneau et al.
(2019) approach is restricted to binary treatments. Second, the outcomes of
interest in most methods are limited to the mean survival time. Jiang et al.
(2017) proposed an optimal dynamic treatment regime estimator that maximizes
the $t$-time survival probability. However, none of the methods permit
choosing the criterion except Wahed & Thall (2013). Third, parametric models
or strong structural assumptions are employed for all of these methods, with
the exception of Goldberg & Kosorok (2012). Strong modeling assumptions, such
as the accelerated failure time model or the proportional hazards models used
by the existing methods, are subject to model-misspecification, may limit the
policy class, and, consequently, may invoke loss in the value of the optimal
policy. Finally, the censoring assumptions of Goldberg & Kosorok (2012) are
restrictive. Censoring is associated with failure time in many medical
applications, but it is often reasonable to assume independence of censoring
after conditioning on patient historical information.

Table 1: Assumptions of the existing and the proposed methods

method | $Q$ | $|\mathcal{A}^{(q)}|$ | criterion | $T$ | $C$ | $\\{\bm{\pi}\\}$
---|---|---|---|---|---|---
the new method | finite | finite | $E[T\wedge\tau],S(t)$ | NP | CI | flexible
Goldberg & Kosorok (2012) | finite | finite | $E[T\wedge\tau]$ | NP | independent | flexible
Huang et al. (2014) | 2+ | 2+ | $E[T]$ | AFT | CI | linear
Jiang et al. (2017) | finite | 2 | $S(t)$ | PH | CI | linear
Simoneau et al. (2019) | finite | 2 | $E[T]$ | AFT | CI | linear
Wahed & Thall (2013) | finite | finite |  $E[T]$, $S(t)$ | AFT | CI | fixed
Xu et al. (2016) | finite | finite | $E[T]$ | BNP | CI | fixed

$Q$, the number of stages; $|\mathcal{A}^{(q)}|$, the number of treatment arms
at stage $q$; criterion, the target value to be optimized; $T$, failure time
model; $C$, censoring assumption; $\\{\bm{\pi}\\}$, policy class; 2+, can be
generalized to more than two; NP, non-parametric; AFT, accelerated failure
time; PH, proportional hazards; BNP, Bayesian non-parametric; CI, conditional
independence; fixed, a fixed number of distinct treatment rules.

In this manuscript, we develop a general dynamic treatment regime estimator
for censored time-to-failure outcomes that addresses all of the limitations
discussed for the existing estimators. Our estimator maximizes either the mean
survival time or the survival probability at a certain time using backward
recursion. As the objective may not be expressed as the sum of intermediate
rewards, a standard Q-learning algorithm is not applicable. Instead, the
conditional survival probability information, rather than the summarized
Q-values, is appended to the previous stage information. A generalized random
survival forest is developed for this task, where survival curves for each
individual, instead of the observed survival or censoring time, is fed into
the random survival forest. A general implementation of our method is
available on CRAN, R package dtrSurv.

The key contributions of our work are the flexibility of the proposed method
and the exposition of its theoretical properties. The proposed method affords
significantly more flexibility than existing methods. The target value can be
either the mean survival or the survival probability at a certain time. The
method allows for a flexible number of treatment stages as well as a flexible
number of treatment arms at each stage. Censoring can be conditionally
independent. The conditional survival probability is modeled nonparametrically
using a random survival forest-based algorithm. And finally, the proposed
method does not require estimating the joint distribution of failure time and
treatment time nor does it require any assumptions on their dependency, though
these are allowed. We further show the consistency of the random survival
forests and the estimated regime values.

Before we give an outline of the paper, we briefly review random survival
forests, an important component of the proposed method. A random survival
forest is a nonparametric survival regression estimator that provides weighted
Kaplan-Meier curves for each point in the feature space with weights given by
adaptive kernel estimates. In general, a random forest is an ensemble of trees
each of which are generated by recursively partitioning the data into two
daughter nodes. In the Breiman et al. (1984) Classification and Regression
Trees, or CART, approach, partitions are made by searching across all
variables and data midpoints, and the best cut point that maximizes the
heterogeneity between the resulting daughter nodes based on the mean
differences is taken. In Breiman (2001) random forests, Breiman et al. (1984)
trees are built and averaged based on multiple bootstrap samples considering
only a random subset of the variables when splitting each node. This
randomization—bootstrapping and random partitioning—brings diversity among
trees, and thus the variability of the estimators can be less than that
obtained in the CART approach. Geurts et al. (2006) added an alternative
randomization mode by restricting the candidate threshold values to a randomly
chosen cut point for each variable.

Survival trees (random survival forests) generalize trees (random forests) by
allowing censored outcomes as inputs and providing covariate-conditional
survival probability estimates as outputs. See Gordon & Olshen (1985) and
Segal (1988) for survival trees and Hothorn et al. (2005), Ishwaran et al.
(2008), and Zhu & Kosorok (2012) for random survival forests. Due to
censoring, heterogeneity between two nodes cannot be ascertained by a simple
mean; alternate tests such as the log-rank test (Peto & Peto, 1972; Gehan,
1965; Tarone & Ware, 1977) must be used. At each terminal node of a survival
tree, the survival probability is estimated by the Kaplan-Meier procedure. A
random survival forest is obtained by averaging the survival trees.

The remainder of the manuscript is organized as follows. In Section 2, we give
notation and describe the dynamic treatment rule estimator. Consistency of the
estimator is derived in Section 3, and its performance is illustrated through
simulations in Section 4. Section 5 is devoted to the application of the
proposed method to a two-stage leukemia clinical trial. We conclude the paper
by discussing future research directions in Section 6.

##  2 The method

###  2.1 Data setup and notation

We assume that each of $n$ patients can have a maximum of $Q$ visits, or
treatment stages, with a maximum study length $\tau>0$. At the $q$th stage,
$q=1,2,...,Q$, the $i$th patient, among the $n^{(q)}\leq n$ available
patients, receives treatment $A^{(q)}_{i}\in\mathcal{A}^{(q)}$, if he or she
has survived and has not dropped out by the beginning of the stage, where
$\mathcal{A}^{(q)}$ is the finite treatment space for the $q$th stage,
$q=1,2,...,Q$. Throughout this manuscript, we often drop the subject index $i$
when it does not cause confusion. Our interest lies in estimating the survival
distribution $S_{\bm{\pi}}$ of the overall failure time $T$ of patients if
they followed a dynamic treatment regime
$\bm{\pi}=(\pi^{(1)},\pi^{(2)},...\pi^{(Q)})^{\top}$, and then finding the
‘best’ rule $\bm{\pi}_{*}$, which is the rule that optimizes a certain
criterion $\phi$. For $\phi$, we consider
$\phi_{\mu}(S)=-\int_{t>0}(t\wedge\tau)dS(t)$ and $\phi_{\sigma,t}(S)=S(t)$
for some $0<t\leq\tau$. Without loss of generality, we assume that a prolonged
time to event is preferred; thus the objective is to maximize $\phi(S)$ over
the regimes.

We allow the time to next treatment $U^{(q)}_{i}$ of the $i$th patient to
depend on the historical information $\bm{H}^{(q)}_{i}\in\mathcal{H}^{(q)}$
available at the beginning of the $q$th stage and the treatment assignment
$A^{(q)}_{i}$. At the $q$th stage, we observe either the $i$th patient’s
failure, advancement to next treatment, or censoring. The times from the
beginning of the stage to each of these events are denoted as $T^{(q)}_{i}$,
$U^{(q)}_{i}$, and $C^{(q)}_{i}$, respectively. The length of the stage
$X^{(q)}_{i}$ for the $i$th patient is defined by the minimum of these times,
and the hypothetically uncensored stage length is defined as
$V^{(q)}_{i}=T^{(q)}_{i}\wedge U^{(q)}_{i}$. And, if $q=Q$,
$V^{(Q)}_{i}=T^{(Q)}_{i}$. We denote the $q$th stage censoring and treatment
indicators as $\delta_{i}^{(q)}=1(V_{i}^{(q)}\leq C_{i}^{(q)})$ and
$\gamma_{i}^{(q)}=1(T^{(q)}_{i}\leq U^{(q)}_{i})$, respectively. Note that
when $C^{(q)}_{i}<V^{(q)}_{i}$, $\gamma_{i}^{(q)}$ is not observable. The
baseline time at stage $q$ is defined as
$B^{(q)}_{i}=\sum_{k=1}^{q-1}X^{(k)}_{i},q>1$, and $B_{i}^{(1)}=0$. Let
$T_{i},$ $C_{i}$, $X_{i}=T_{i}\wedge C_{i},$ and $\delta_{i}=1(T_{i}\leq
C_{i})$ denote the overall failure time, overall censoring time, overall
observed time, and overall censoring indicator, respectively. We denote
$\bm{Z}^{(q)}_{i}\in\mathcal{Z}^{(q)}$ as the covariate information of the
$i$th patient that is newly available at the beginning of the $q$th stage.
Thus, $\bm{H}^{(q)}_{i}$ may include the historical information such as
$A^{(k)}_{i}$, $B^{(k)}_{i}$, $\bm{Z}^{(k)}_{i}$, $k=1,2,..,q-1$ and
$\bm{Z}^{(q)}_{i}$; where $d^{(q)}$ denotes the dimension of
$\bm{H}^{(q)}_{i}~{}$ for all $i=1,2,...,n^{(q)};q=1,2,..,Q$. The number
$n^{(q)}$ of patients eligible for treatment assignment at the beginning of
stage $q$ is $\sum_{i=1}^{n}1(U^{(q)}_{i}\leq T^{(q)}_{i}\wedge C^{(q)}_{i})$.
The notation is summarized in Table 3 in Appendix 1.

###  2.2 Overview of the method

As is the case in Q-learning, our estimator optimizes the intermediate stage
decision assuming that all the later decisions are optimal, and thus backward
recursion is used. Our method is based on the observation that the remaining
life $L^{(q)}$ at stage $q$ is defined recursively as a convolution of
$L^{(q+1)}$ and the current stage length $V^{(q)}$ and that

$\displaystyle S^{(q)}(t\mid\bm{H}^{(q)},A^{(q)},\delta^{(q)}=1)=$
$\displaystyle\int
S^{(q+1)}(t-V^{(q)}\mid\bm{H}^{(q+1)},A^{(q+1)},\delta^{(q)}=1)$
$\displaystyle\hskip
28.45274ptdP(V^{(q)},\bm{H}^{(q+1)},A^{(q+1)}\mid\bm{H}^{(q)},A^{(q)},\delta^{(q)}=1),$
(1)

where $S^{(q)}(t\mid\cdot)=1$ for all $t<0$ and
$S^{(q)}(t\mid\gamma^{(q-1)}=1)=1(t<0)$, for $q=Q-1,Q-2,...,1.$ In other
words, assuming no censoring events, the survival probability $S^{(q)}$ of the
remaining life $L^{(q)}$ at stage $q$ can be obtained by augmenting the
survival probability of the remaining life $L^{(q+1)}$ of the next stage by
the current stage length $V^{(q)}$ and then marginalizing over all events that
occur later than the beginning of stage $q$. Under the covariate-independent
censoring assumption, the contribution of the censored cases is incorporated
into the equation by redistributing the probability to the right (Efron,
1967):

$\displaystyle S^{(q)}(t\mid\bm{H}^{(q)},A^{(q)})$ $\displaystyle=$
$\displaystyle\int\bigg{[}\delta^{(q)}S^{(q+1)}(t-X^{(q)}\mid\bm{H}^{(q+1)},A^{(q+1)})$
$\displaystyle~{}~{}~{}~{}~{}+(1-\delta^{(q)})\Big{\\{}1(X^{(q)}>t)+1(X^{(q)}\leq
t)\frac{S^{(q)}(t\mid\bm{H}^{(q+1)},A^{(q+1)})}{S^{(q)}(X^{(q)}\mid\bm{H}^{(q+1)},A^{(q+1)})}\Big{\\}}\bigg{]}$
$\displaystyle~{}~{}~{}\times
dP(X^{(q)},\bm{H}^{(q+1)},A^{(q+1)},\delta^{(q)}\mid\bm{H}^{(q)},A^{(q)}).$
(2)

Assuming the standard causal inference assumptions, which we will formally
introduce in Subsection 3.2 (Assumption 8–10), we can translate
$S^{(q)}(t\mid\bm{H}^{(q)},A^{(q)}=a)$ into the counterfactual quantity
$S^{(q)}_{a}(t\mid\bm{H}^{(q)})$. Thus, the regime is estimated using the
following procedure: 1. estimating the conditional survival distribution
$S^{(q)}(t\mid\bm{H}^{(q)},A^{(q)})$ of the remaining life $L^{(q)}$ for
$q=Q$; 2. optimizing the survival distribution over $A^{(q)}$ to have
$\hat{\pi}^{(q)}$ and $\hat{S}^{(q)}_{*}(t\mid\bm{H}^{(q)})$ for $q=Q$; 3.
augmenting the estimated optimal curves by the previous stage lengths, i.e.,
$\hat{S}^{(q)}_{*}(t-X^{(q-1)}\mid\bm{H}^{(q)},A^{(q)})$ for $q=Q$; and, 4.
repeating steps 1-3 for $q=Q-1,Q-2,...,1$. This procedure is summarized at a
high level in Algorithm 1.

Result: A dynamic treatment regime estimate
$\hat{\bm{\pi}}=(\hat{\pi}^{(1)},\hat{\pi}^{(2)},...,\hat{\pi}^{(Q)})^{\top}$

for  _stage $q=Q,Q-1,...,1$_ do

Obtain $\hat{S}^{(q)}(\cdot\mid\bm{H}^{(q)},A^{(q)})$ via generalized random
survival forest (Algorithm 2);

Obtain
$\hat{\pi}^{(q)}(\bm{h})=\arg\max_{a\in\mathcal{A}^{(q)}}\phi^{(q)}\Big{\\{}\hat{S}^{(q)}(\cdot\mid\bm{H}^{(q)}=\bm{h},A^{(q)}=a)\Big{\\}}$;

Obtain
$\hat{S}^{(q)}_{*}(\cdot\mid\bm{H}^{(q)})=\hat{S}^{(q)}(\cdot\mid\bm{H}^{(q)},A^{(q)}=\hat{\pi}^{(q)}(\bm{H}^{(q)}))$;

Algorithm 1 the proposed dynamic treatment regime estimator

At the beginning of the last stage, there are $n^{(Q)}$ patients who have
neither experienced failure nor been censored. The data for these patients are
used to estimate the terminal stage covariate-conditional survival probability
$S^{(Q)}(\cdot\mid\bm{H}^{(Q)},A^{(Q)})$ for each treatment arm using random
survival forests described in detail in Subsection 2.4. Then the optimal
regime for the final stage is estimated by
$\hat{\pi}^{(Q)}=\arg\max_{A}\phi^{(Q)}\Big{\\{}\hat{S}^{(Q)}(\cdot,A)\Big{\\}},$
where $\hat{S}$ is the estimate of $S$ and $\phi^{(q)},q=2,..,Q$ is defined
adaptively for the intermediate and terminal survival probabilities.
Specifically,
$\phi_{\mu}^{(q)}\Big{\\{}S^{(q)}(\cdot\mid\bm{H}^{(q)})\Big{\\}}=B^{(q)}-\int_{t>X^{(q)}}(t\wedge\tau)dS^{(q)}(t\mid\bm{H}^{(q)})$
and
$\phi_{\sigma,t}^{(q)}\Big{\\{}S^{(q)}(\cdot\mid\bm{H}^{(q)})\Big{\\}}=S^{(q)}(t-B^{(q)}\mid\bm{H}^{(q)})$.

We denote the estimated optimized survival distribution of the $q$th stage as
$\hat{S}^{(q)}_{*}=\hat{S}^{(q)}\\{\cdot,\hat{\pi}^{(q)}(\cdot)\\}$. The
optimized survival information of the final stage is carried back to the
previous stage. Specifically, the $i$th patient who is available at the
$(Q-1)$st stage and has a stage length $X_{i}^{(Q-1)}$ is given a
probabilistic augmentation so that the resulting survival probability of
remaining life at the second to the last stage is
$\hat{S}^{(Q-1)}_{*}(t-X_{i}^{(Q-1)}\mid\bm{H}^{(Q)}_{i})$. For those who have
already experienced failure or censoring by the end of the stage, no
augmentation is required, as either their remaining life has been determined
already during the stage or the censoring indicator and the stage length
contain all the necessary information.

Next, the conditional distribution of remaining life from the beginning of the
$(Q-1)$st stage is estimated by the random survival forests based on the
$n^{(Q-1)}$ patients. The distinction between the last stage and all the other
stages is that given data, the patients’ remaining life conditioning on the
current historical information is observed in a scalar form at the last stage,
while the remaining life for the other stages is not observed but is given in
the form of a stochastic process. So a modification to the random survival
forest estimator is required and is discussed in Subsection 2.4. Then
estimation of $S^{(q)},\pi^{(q)},$ and $S^{(q)}_{*}$ is recursively performed
until we reach the first stage and obtain
$\hat{\bm{\pi}}=(\hat{\pi}^{(1)},\hat{\pi}^{(2)},...,\hat{\pi}^{(Q)})^{\top}$.

###  2.3 A few aspects of the optimization: Backward recursion and composite
values

We discuss a few aspects of the optimization problem. First, we go over why
backward recursion can be used in optimizing the truncated survival mean
values ($\phi_{\mu}$) or survival probabilities ($\phi_{\sigma,t}$). For some
criteria, such as the median $\phi_{med}(S)=S^{-1}(0.5)$, this backward
recursion does not guarantee an overall optimum. However, for the mean
survival time $\phi_{\mu}$ and the survival probability $\phi_{\sigma,t}$ at
time $t$, backward recursion is a legitimate method according to Proposition
1, which is proved in Appendix 2.

######  Proposition 1 (Legitimacy of backward recursion).

Pick either $\phi_{\mu}$ or $\phi_{\sigma,t}$ for $\phi$. Let $\\{\bm{\pi}\\}$
be an arbitrary class of decision rules and
$\bm{\pi}_{*}=(\pi^{(1)}_{*},\pi^{(2)}_{*},...,$ $\pi^{(Q)}_{*})^{\top}$
satisfy
$\phi\Big{\\{}S(\cdot\mid\bm{H}^{(1)},\bm{\pi}_{*})\Big{\\}}\geq\phi\Big{\\{}S(\cdot\mid\bm{H}^{(1)},\bm{\pi})\Big{\\}}$
for all $\bm{\pi}\in\\{\bm{\pi}\\}.$ Then for any $q=1,2,...,Q-1,$

$\pi^{(q)}_{*}(\bm{H}^{(q)})=\arg\max_{a^{(q)}\in\mathcal{A}^{(q)}}\phi^{(q)}\Big{\\{}S^{(q)}(\cdot\mid\bm{H}^{(q)},\pi^{(q+1)}_{*},\pi^{(q+2)}_{*},...,\pi^{(Q)}_{*})\Big{\\}}.$

Second, we study an alternative to the $\phi_{\sigma,t}$ criterion. When
$\phi_{\sigma,t}(S)=S(t)$ is chosen as the optimization criterion, the
treatment rule can not discriminate treatment arms at any point later than
$t$. For example, when an investigator is interested in the six month survival
probability and a patient takes her third round of treatment seven months
after her onset of the disease, the optimal dynamic treatment rule could
recommend any treatment option randomly, since it does not affect the
patient’s six month survival probability. Thus a composite criterion of first
maximizing $S(t)$ and second maximizing $E[T\wedge\tau]$ could be more
beneficial than the simple criterion $\phi_{\sigma,t}$. The composite
criterion values $\phi_{\sigma,\mu,t}(S)=\\{S(t),E[T\wedge\tau]\\}$ are
ordered lexicographically; i.e.,
$\phi_{\sigma,\mu,t}(S_{\bm{\pi}_{1}})>\phi_{\sigma,\mu,t}(S_{\bm{\pi}_{2}})$
if and only if either of the following is true:

$\displaystyle S_{\bm{\pi}_{1}}(t)>S_{\bm{\pi}_{2}}(t),$ or $\displaystyle
S_{\bm{\pi}_{1}}(t)=S_{\bm{\pi}_{2}}(t)\text{ and
}E[S_{\bm{\pi}_{1}}\wedge\tau]>E[S_{\bm{\pi}_{2}}\wedge\tau].$

###  2.4 A generalized random survival forest

We develop a generalized random survival forest that takes survival
probabilities as the outcome variable. The high level algorithm is provided in
Algorithm 2. Suppose we are interested in estimating the marginal survival
probability of a node $\mathcal{N}$ at stage $q$, where each individual’s
failure time information is given by the survival probability of the remaining
life and the censoring status:
$\\{(S_{i}^{(q)},\delta_{i}^{(q)})\\}_{i=1}^{n_{\mathcal{N}}}$ where
$n_{\mathcal{N}}$ is the size of node $\mathcal{N}$ and $S_{i}^{(q)}$ is
shorthand for $\hat{S}^{(q+1)}(t-X^{(q)}_{i}\mid\bm{H}^{(q+1)}_{i})$. Then we
modify the Kaplan-Meier estimator as

$\displaystyle\hat{S}^{(q)}(t\mid\mathcal{N})$ $\displaystyle=\prod_{s\geq
0}^{t}\left\\{1+\frac{\sum_{i\in\mathcal{N}}\delta_{i}^{(q)}dS_{i}^{(q)}(s)}{\sum_{i\in\mathcal{N}}S_{i}^{(q)}(s^{-})}\right\\},$
(3)

where $s^{-}=s-\epsilon$ with an arbitrarily small positive value $\epsilon$.
This modified Kaplan-Meier estimator will be used in both tree partitioning
and terminal node survival probability estimation.

Result: generalized random survival forest for the $q$th stage
$\hat{S}^{(q)}(\cdot\mid\bm{H}^{(q)},A^{(q)})$

Parameters: the number of trees $n_{\text{tree}}$, the minimum terminal node
size $n_{\text{min}}$, the minimum number of events in terminal nodes
$n_{\text{min.E}}$;

Input data:
$\\{(X^{(q)}_{i},\gamma^{(q)}_{i},\delta^{(q)}_{i},\bm{H}^{(q)}_{i},A^{(q)}_{i})\\}_{i=1}^{n^{(q)}}$,
$\hat{S}^{(q+1)}_{*}(\cdot\mid\bm{H}^{(q+1)})$;

Stochastic augmentation:
$S^{(q)}_{i}=\hat{S}^{(q+1)}_{*}(t-X^{(q)}_{i}\mid\bm{H}^{(q+1)}=\bm{h}^{(q+1)}_{i})$
with $S^{(q)}_{i}(t)=1(t\leq X^{(q)}_{i})$ if $q=Q$ or if
$\delta^{(q)}_{i}=0$;

for  _treatment arm $a\in\mathcal{A}^{(q)}$_ do

for  _tree $b=1,2,...,n_{\text{tree}}$_ do

Recursively partition the feature space with the augmented data
$\\{(S^{(q)}_{i},\delta^{(q)}_{i},\bm{H}^{(q)}_{i})\\}_{\\{i:a^{(q)}_{i}=a\\}}$
via LR (4) or MD (5) until node size $<2n_{\text{min}}$ or number of failures
in each node $<2n_{\text{min.E}}$;

Obtain the terminal node survival probability
$\hat{S}^{(q)}_{b,a}(t\mid\bm{H}^{(q)})$ via modified Kaplan-Meier (3);

Obtain the random survival forest estimator
$\hat{S}^{(q)}(\cdot\mid\bm{H}^{(q)},A^{(q)}=a)=\frac{1}{n_{\text{tree}}}\sum_{b=1}^{n_{\text{tree}}}\hat{S}^{(q)}_{b,a}(t\mid\bm{H}^{(q)})$;

Algorithm 2 the generalized random survival forest estimator

Now we develop the splitting rules. Many random survival forest methods use
the two-sample log-rank test statistic to partition the population into
subgroups with a homogeneous survival distribution (Ishwaran et al., 2008;
Hothorn et al., 2005). When failure time is not observed as a scalar but is
given as a probabilistic process, the log-rank test can be generalized to the
generalized log-rank statistic (Cho et al., 2019);

$\displaystyle
LR(\mathcal{N}_{1},\mathcal{N}_{2})=\frac{\int_{0}^{\tau}\frac{Y_{2}(t){dN_{1}(t)}+Y_{1}(t){dN_{2}(t)}}{Y(t)}}{\sqrt{\int_{0}^{\tau}\frac{Y_{1}(t)Y_{2}(t)dN(t)(Y(t)-dN(t))}{Y(t)^{3}}}},$
(4)

where $\mathcal{N}_{1}$ and $\mathcal{N}_{2}$ are two candidate daughter
nodes,
$Y_{l}(t)=\frac{1}{|\mathcal{N}_{l}|}\sum_{i\in\mathcal{N}_{l}}S^{(q)}_{i}(t^{-}),N_{l}(t)=1-\frac{1}{|\mathcal{N}_{l}|}\sum_{i\in\mathcal{N}_{l}}S_{i}(t),l=1,2$,
$Y(t)=\lambda_{1}Y_{1}(t)+\lambda_{2}Y_{2}(t),$
$N(t)=\lambda_{1}N_{1}(t)+\lambda_{2}N_{2}(t),$
$\lambda_{l}=\frac{|\mathcal{N}_{l}|}{|\mathcal{N}_{1}|+|\mathcal{N}_{2}|},$
and $|\mathcal{N}|=\sum_{i=1}^{n^{(q)}}1(i\in\mathcal{N})$. The generalized
log-rank splitting rule chooses
$(\mathcal{N}_{1},\mathcal{N}_{2})=\arg\max_{\mathcal{N}_{1},\mathcal{N}_{2}}LR(\mathcal{N}_{1},\mathcal{N}_{2})$
as the best partition.

When the optimization criterion is the truncated mean survival time
$\phi_{\mu}$, splitting can also be done by finding the split point that
maximizes the difference of the truncated mean survival times. Define

$\displaystyle
MD(\mathcal{N}_{1},\mathcal{N}_{2})=\bigg{|}\int_{t>0}(t\wedge\tau)d\hat{S}(t\mid\mathcal{N}_{1})-\int_{t>0}(t\wedge\tau)d\hat{S}(t\mid\mathcal{N}_{2})\bigg{|},$
(5)

where $\hat{S}$ is the modified Kaplan-Meier estimator defined in Equation
(3). The mean difference splitting rule finds
$(\mathcal{N}_{1},\mathcal{N}_{2})=\arg\max_{\mathcal{N}_{1},\mathcal{N}_{2}}MD(\mathcal{N}_{1},\mathcal{N}_{2})$
as the best split.

The recursive partitioning continues until each terminal node size is less
than a predefined threshold $k$. At each terminal node $\mathcal{N}$ of the
$b$th survival tree, the survival probability estimate
$\hat{S}^{(q)}_{b}(t\mid\mathcal{N})$ is given by the modified Kaplan-Meier
estimator in (3), $b=1,2,...,n_{\text{tree}}$, where $n_{\text{tree}}$ is the
number of trees in the forest. The random survival forest estimate is then
given by

$\displaystyle\hat{S}^{(q)}(t\mid\bm{h}^{(q)})$
$\displaystyle=\frac{1}{n_{\text{tree}}}\sum_{b=1}^{n_{\text{tree}}}\sum_{l=1}^{\\#\\{\mathcal{N}_{b,l}\\}_{l}}1(\bm{h}^{(q)}\in\mathcal{N}_{b,l})\hat{S}^{(q)}(t\mid\mathcal{N}_{b,l}),$

where $\\#\\{\mathcal{N}_{b,l}\\}_{l}$ is the size, or the number of terminal
nodes, of the $b$th tree.

##  3 Theoretical properties

We show consistency of the proposed dynamic treatment regime estimator
$\hat{\bm{\pi}}$ using empirical process theory. We first derive that the
generalized random survival forest is a uniformly—over the feature
space—consistent estimator. The uniformity is required since estimating the
optimal treatment policy of the $q$th stage relies on integrating the
estimated survival curve of ($q+1$)st stage across the newly accrued
information $(\bm{H}^{(q+1)},T^{(q)},U^{(q)})\mid\bm{H}^{(q)}$ between the two
stages. Then we show that the value of the estimated regime is consistent for
the value of the optimal regime, where the versions of the value are defined
in Subsection 3.2.

###  3.1 Consistency of the generalized random survival forests

We derive consistency of the generalized random survival forest outside of the
context of dynamic treatment regimes. We have the following assumptions on
data distribution that are commonly used to achieve consistency of random
forests or random survival forests (Meinshausen, 2006; Wager & Walther, 2015;
Cui et al., 2017; Wager & Athey, 2018). These assumptions are tailored for use
later in the dynamic treatment regime setting by adding stage indicators in
Subsection 3.2.

######  Assumption 1 (Lipschitz continuous survival probability and censoring
density).

There exist constants $L_{S}$ and $L_{G}$ such that
$|S(t\mid\bm{h}_{1})-S(t\mid\bm{h}_{2})|\leq
L_{S}\|\bm{h}_{1}-\bm{h}_{2}\|_{1}$ and
$|g(t\mid\bm{h}_{1})-g(t\mid\bm{h}_{2})|\leq
L_{G}\|\bm{h}_{1}-\bm{h}_{2}\|_{1}$ for all
$\bm{h}_{1},\bm{h}_{2}\in\mathcal{H}$ and $t\in[0,\tau]$, where $G$ is the
censoring survival distribution and $g$ is its derivative with respect to
time.

######  Assumption 2 (Weakly dependent covariate values).

The covariate information ($\bm{Z}$) is given as a $d_{Z}$-dimensional vector
residing in a unit hypercube, i.e., $\bm{Z}\in[0,1]^{d_{Z}}$, with a density
$f_{\bm{Z}}$ such that $\zeta^{-1}\leq f_{\bm{Z}}(\bm{z})\leq\zeta$ for all
but countably many points $\bm{z}\in[0,1]^{d_{Z}}$ and some constant
$\zeta\geq 1$.

######  Remark 1.

The distributional assumption was made for notational convenience and could be
relaxed. For example, the results in this paper still hold for covariates that
are a multivariate normal random vector truncated by some constants since
random forests are invariant to scale transformation. Also the results are
still valid with categorical covariates with a finite number of levels after
modification to Assumption 1: in this setting, the norm is calculated by
excluding the categorical covariate coordinates. The case of categorical
covariates will be discussed in the proof of Theorem 1.

We also have assumptions regarding the tree-based estimator terminal node size
and splitting rules. The first assumption forces the terminal node size to be
both sufficiently small and sufficiently large.

######  Assumption 3 (Terminal node size).

The maximum size $k$ of the terminal nodes grows at the following rate:

$k\asymp n^{\beta},\frac{1}{2}<\beta<1,$

where $a\asymp b$ implies both $a=\mathcal{O}(b)$ and $b=\mathcal{O}(a)$.

Next, we give definitions of regular trees and random split trees that are
often used for consistency proofs in the random forest literature
(Meinshausen, 2006; Wager & Walther, 2015; Cui et al., 2017; Wager & Athey,
2018).

######  Definition 1 (random-split and $\alpha$-regular trees and forests).

A tree is called a random-split tree if each feature is given a minimum
probability ($\varphi/d$) of being the split variable at each intermediate
node, where $0<\varphi<1$. A tree is $\alpha$-regular, if every daughter node
has at least $\alpha$ fraction of the training sample in the corresponding
parent node. A random forest is called a random-split ($\alpha$-regular)
forest, if each of its member trees is random-split ($\alpha$-regular).

######  Assumption 4 ($\alpha$-regular and random split trees).

Trees are $\alpha$-regular and random-split with a constant $0<\varphi<1$.

Now we give our uniform consistency results of the generalized random survival
forests. This theorem is largely based on the Z-estimation lemma and
Donskerness of the tree kernels. A detailed proof is relegated to Appendix 3.

######  Theorem 1 (Uniform consistency of a random survival forest).

Suppose Assumptions 1-4 hold. Let $\tilde{S}^{\prime}$ be an estimator of
$S^{\prime}$ that is uniformly consistent in both $\bm{h}\in\mathcal{H}$ and
$t\in[0,\tau]$, where $S^{\prime}$ satisfies $S(t\mid\bm{h})=\int
S^{\prime}(t-(T\wedge
U)\mid\bm{H}^{\prime})dP(T,U,\bm{H}^{\prime}\mid\bm{H}=\bm{h})$. Then, the
generalized random survival forest $\hat{S}(t\mid\bm{h})$ built based on
$\tilde{S}^{\prime}$ is uniformly consistent. Specifically,

$\sup_{t\in[0,\tau],\bm{h}\in\mathcal{H}}|\hat{S}(t\mid\bm{h})-S(t\mid\bm{h})|\to
0,$

in probability as $n\to\infty.$

Consistency of a random survival forest sequence in a dynamic treatment regime
setting requires adaptive versions of Assumptions 1-3 and an additional
condition that guarantees an arbitrarily large sample size for the last stage
with high probability. To achieve this condition, we assume that, at each
stage except the last one, each patient receives the next treatment before
experiencing failure or censoring with at least a certain minimum probability.

######  Assumption 5 (Lipschitz continuous survival probability and censoring
density at each stage).

There exist constants $L_{S}$ and $L_{G}$ such that
$|S^{(q)}(t\mid\bm{h}_{1})-S^{(q)}(t\mid\bm{h}_{2})|\leq
L_{S}\|\bm{h}_{1}-\bm{h}_{2}\|$ and
$|g^{(q)}(t\mid\bm{h}_{1})-g^{(q)}(t\mid\bm{h}_{2})|\leq
L_{G}\|\bm{h}_{1}-\bm{h}_{2}\|$ for all
$\bm{h}_{1},\bm{h}_{2}\in\mathcal{H}^{(q)}$, $t\in[0,\tau-B^{(q)}]$, and
$q=1,2,...,Q$, where $G^{(q)}$ is the censoring survival distribution at the
$q$th stage and $g^{(q)}$ is its derivative with respect to time.

######  Assumption 6 (Weakly dependent historical information at each stage).

The patient history information $\bm{H}^{(q)}$ at stage $q$ is given as a
$d^{(q)}$-dimensional vector lying in a subset $\mathcal{H}^{(q)}$ of
$[0,1]^{d^{(q)}}$, and its joint and marginal densities are bounded so that
$\zeta^{-1}\leq f_{\bm{H}^{(q)}}(\bm{h})\leq\zeta$ and $\zeta^{-1}\leq
f_{\bm{H}^{(q)},j}(\bm{h}_{j})\leq\zeta$ for all but a finite number of
$\bm{h}\in\mathcal{H}^{(q)}$ and $j=1,2,...,d^{(q)}$ for some constant
$\zeta\geq 1$.

Assumption 6 modifies Assumption 2 so that the support of the history space
can be less than a hyper-cube, since, for example, the baseline survival
times, $(B^{(1)}\leq B^{(2)}\leq B^{(3)}\leq...)$, have a dependent support
structure.

######  Assumption 7 (Probability of proceeding to the next stage).

There exists a constant $M\in(0,1)$ such that for all
$h^{(q)}\in\mathcal{H}^{(q)}$ and $q=1,2,...,Q-1$,
$\textnormal{pr}(U^{(q)}<T^{(q)}\wedge C^{(q)}\mid\bm{h}^{(q)})>M.$

######  Corollary 1 (Uniform consistency of a random survival forest
sequence).

Suppose Assumptions 3–7 hold. Let
$\hat{\bm{S}}=(\hat{S}^{(1)},...,\hat{S}^{(q)},...,\hat{S}^{(Q)})$ be the
sequence of the generalized random survival forest estimators of
$\bm{S}=(S^{(1)},...,S^{(q)},...,S^{(Q)})$ such that the $q$th stage random
survival forest is built based on $\hat{S}^{(q+1)}$ for $q=1,2,...,Q-1$. Then,

$\sup_{t\in[0,\tau],\bm{h}\in\mathcal{H},q\in\\{1,2,...,Q\\}}|\hat{S}^{(q)}(t\mid\bm{h})-S^{(q)}(t\mid\bm{h})|\to
0,$

in probability as $n\to\infty.$

The proof of Corollary 1 is given in Appendix 3.

###  3.2 Consistency of the dynamic treatment regime estimator

We show consistency of the dynamic treatment regime estimator using the causal
inference framework. With this framework, the theoretical properties are
applicable to broader settings than sequentially randomized experiments. A
counterfactual outcome is defined as an outcome that would have been obtained
if a person had a treatment option contrary to his/her actual treatment. To
denote the counterfactual outcomes pertaining to a counterfactual treatment or
treatment policy, we add a subscript to the corresponding random variables or
survival functions. For example, $S_{\pi}^{(q)}(t|\bm{H}^{(q)})$ is a
counterfactual survival probability of the remaining life if treatment was
given according to a treatment rule $\pi$ to a person with history
$\bm{H}^{(q)}$ at stage $q$. We assume the following three standard causal
inference assumptions (Hernán & Robins, 2020; Rubin, 2005; Cole & Frangakis,
2009) known as the stable unit treatment value assumption, sequential
ignorability, and positivity.

######  Assumption 8 (Stable unit treatment value assumption, or SUTVA).

Each person’s counterfactual dynamics, such as failure time, time to next
treatment, and censoring, are not affected by the treatments or history of
other patients. Moreover, each treatment has only one version, or if there are
multiple versions, their differences are irrelevant.

######  Assumption 9 (Sequential ignorability).

Treatment assignment is given independently of the counterfactual outcomes,
conditionally on the individual’s historical information, where the outcomes
include all future random quantities such as failure time, time to next
treatment, and censoring time.

######  Assumption 10 (Positivity).

For each stage $q=1,2,...,Q$, given historical information, the probability of
having each treatment $a\in\mathcal{A}^{(q)}$ is greater than a constant
$L>0$. Or,

$\textnormal{pr}(A^{(q)}=a|\bm{H}^{(q)}=\bm{h})>L$

for all $a\in\mathcal{A}^{(q)}$, $\bm{h}$ such that
$\textnormal{pr}(\bm{H}^{(q)}=\bm{h})>0.$

The following theorem states that the dynamic treatment regime estimator has a
value consistent for the value of the optimal regime. We define the value
$\Phi$ of a treatment regime $\bm{\pi}$ as the criterion value of the survival
probability if all individuals in a population follow the treatment regime. We
use $\Phi(\bm{\pi})=\phi(S_{\bm{\pi}})$ to denote the value of $\bm{\pi}$, and
$\Phi$ inherits the subscript of $\phi$.

######  Theorem 2 (Consistency of the dynamic treatment regime estimator).

Suppose Assumptions 5–10 hold, and the optimal dynamic treatment regime
estimator ($\hat{\bm{\pi}}$) is built based on the generalized random forest
grown under Assumptions 3 and 4. Then, for each
$\phi=\phi_{\mu},\phi_{\sigma,t}$,

$|\Phi(\bm{\pi})-\Phi(\hat{\bm{\pi}})|\to 0,$

in probability as $n\to\infty.$

The proof is provided in the Appendix 4.

##  4 Simulations

###  4.1 Simulation settings

We perform simulations to study the performance of the proposed method and
compare those results to that of existing methods. The overall scheme of the
simulations follows the pattern in Goldberg & Kosorok (2012), where they
mimicked a cancer clinical trial using the tumor size and wellness dynamics.
We adapt their data generating mechanism so that dependence among censoring,
tumor size, and wellness is incorporated in the mechanism, and deterministic
processes are relaxed into more realistic ones.

In the hypothetical trial, patients can go through up to $Q=3$ treatment
rounds during the trial with length of the trial, $\tau=10$. Each patient has
tumor size $\phi(t)$ and wellness $\omega(t)$ at time $t\in[0,10]$ that
affects the timing of failure or treatment. We use $[t^{(q)},t^{(q+1)}]$ to
denote the $q$th stage. For each stage, patients are given either a more
aggressive ($A=1$) or a less aggressive ($A=0$) treatment.

At the beginning of the trial, patients have uniformly distributed tumor sizes
and wellness: $\rho(0)\sim U(0.5,1),\omega(0)\sim U(0.5,1)$. At each treatment
stage, the more aggressive treatment immediately reduces the tumor size more
than the less aggressive one. However, the more aggressive treatment worsens
patient wellness more than the less aggressive one:

$\displaystyle\rho(t^{(q)}+)$
$\displaystyle=\frac{\rho(t^{(q)})}{\omega(t^{(q)})\\{10-6A\\}}\vee 0$
$\displaystyle\omega(t^{(q)}+)$ $\displaystyle=\omega(t^{(q)})-2^{A-2}.$

After each treatment ($t^{(q)}<t\leq t^{(q+1)}$), the tumor size and wellness
follow the Ornstein-Uhlenbeck processes with negative dependence between the
two:

$\displaystyle\frac{d\chi(t)}{dt}=-2\chi(t)+\eta_{\chi}(t),$

where $\chi=\rho,\omega$ and the second term obeys a drifted Wiener process
such that

$\displaystyle\frac{d}{dt}\begin{pmatrix}\eta_{\rho}(t)\\\
\eta_{\omega}(t)\end{pmatrix}\sim N\Bigg{(}\begin{pmatrix}0.1\\\
0.05\end{pmatrix},\begin{pmatrix}1&-0.5\\\ -0.5&1\end{pmatrix}\Bigg{)}.$

The positive means in the normal distribution imply that on average both tumor
size and wellness grow. When tumor size is greater than 1, the patient
receives the next treatment. In other words,
$t^{(q+1)}=\min_{t}\arg\max_{t>t^{(q)}}1(\omega(t)>1)$, for $q=1,2,..,Q-1$.
The patient failure time is governed by the hazard process:

$\displaystyle\lambda_{F}(c)$
$\displaystyle=\frac{1}{5}\Bigg{\\{}\frac{\rho(u)}{\omega(u)}e^{\bm{g}(H^{(q)})^{\top}\bm{\beta}_{F}(A^{(q)})}+1\Big{(}\omega(u)<0.1\Big{)}\Bigg{\\}},$
(6)

where $g(H^{(q)})=\\{1,\log(B^{(q)}+1),{\bm{Z}^{(q)}}^{\top}\\}^{\top}$,
$\bm{\beta}_{F}(a)$ is the hazard ratio for the failure time given treatment
$a$ that is defined in Table 2, and the indicator accelerates the failure time
when the wellness drops below the critical level, 0.1. The baseline covariates
at the $q$th stage follow a $p-$variate normal distribution conditioning on
the previous stage value:
$\bm{Z}^{(q)}=\frac{1}{2}\bm{Z}^{(q-1)}+\frac{1}{2}\bm{Z}_{0}^{(q)}$, where
$\bm{Z}_{0}^{(q)}\sim N_{p}(\bm{0}_{p},0.8I_{p}+0.2\cdot
1_{p}^{p}),q=1,2,...,Q$, $1_{p}^{p}$ is the $p$-dimensional square matrix of
ones, and $\bm{Z}^{(0)}=\bm{0}_{p}$.

For each stage, the censoring occurs according to the hazard function,

$\displaystyle\lambda_{C}(u)$
$\displaystyle=e^{\bm{g}(H^{(q)})^{\top}\bm{\beta}_{C}(A^{(q)})},$ (7)

where the censoring hazard ratio, $\bm{\beta}_{C}(a)$, is also defined in
Table 2.

Table 2: The failure time and censoring hazard coefficients and dimensions of
covariates in Scenarios 1 to 4.

| scenario 1 | scenario 2 | scenario 3 | scenario 4
---|---|---|---|---
| (base) |  (small $p$) |  (large $p$) | (moderate censoring)
$p$ | 5 | 2 | 10 | 5
$\bm{\beta}_{F}(1)$ | $(0,-3,2,2,1,-1,-1)^{\top}$ | $(0,-2,2,0)^{\top}$ | $(0,-2,\bm{2}_{4}^{\top},-\bm{1}_{4}^{\top},-\bm{2}_{2}^{\top})$ | $(0,-3,2,2,1,-1,-1)^{\top}$
$\bm{\beta}_{F}(0)$ | $(0,-1,1,1,1,1,1)^{\top}$ | $(0,-1,1,1)^{\top}$ | $(0,-1,\bm{1}_{4}^{\top},\bm{1}_{4}^{\top},\bm{0}_{2}^{\top})$ | $(0,-1,1,1,1,1,1)^{\top}$
$\bm{\beta}_{C}(1)$ | $(-3,.2,\bm{.2}_{3}^{\top},-\bm{.2}_{2}^{\top})^{\top}$ | $(-3,.2,.2,-.2)^{\top}$ | $(-3,.2,\bm{.2}_{4}^{\top},-\bm{.2}_{3}^{\top},\bm{0}_{3}^{\top})^{\top}$ | $(-2,.2,\bm{.2}_{3}^{\top},-\bm{.2}_{2}^{\top})^{\top}$
$\bm{\beta}_{C}(0)$ | $(-3,.1,\bm{.2}_{3}^{\top},\bm{.2}_{2}^{\top})^{\top}$ | $(-3,.1,.2,.2)^{\top}$ | $(-3,.1,\bm{.2}_{4}^{\top},\bm{.2}_{3}^{\top},\bm{0}_{3}^{\top})^{\top}$ | $(-2,.1,\bm{.2}_{3}^{\top},\bm{.2}_{2}^{\top})^{\top}$

$\bm{\beta}_{F}(\cdot)$ and $\bm{\beta}_{C}(\cdot)$, the log hazard ratios of
the counterfactual failure and censoring times in (6) and (7), respectively.

The treatment is allocated according to a Bernoulli trial with propensity

$\displaystyle\pi(\bm{H}^{(q)})=\Big{\\{}1+\exp(-g(\bm{H}^{(q)})^{\top}\bm{\beta}_{\pi})\Big{\\}}^{-1},$
(8)

where the coefficient is
$\bm{\beta}_{\pi}=(-1,1,-\frac{1}{2}\bm{1}_{p}^{\top})^{\top}$ for the
observational data setting and $\bm{\beta}_{\pi}=\bm{0}_{p+2}$ for the
randomized trial setting.

Simulations are run with two sample sizes, $n=300$ and $1000$. Thus, we have
in total $16$ distinct factorial settings (four hazard/censoring scenarios,
two propensity designs, and two sample sizes). For each setting, we have
$n_{rep}=200$ simulation replicates, and we estimate the optimal treatment
rule based on both the $\phi_{\mu}$ and $\phi_{\sigma,t=5,\mu}$ criteria. We
compare the proposed method with the methods in Goldberg & Kosorok (2012) and
Simoneau et al. (2019). We further compare the results with the observed
policy characterized by (8) and the estimated optimal zero-order model, where
all patients are given the same embedded treatment regime, which is the best
on average over all embedded regimes. In particular, in implementing the
Goldberg & Kosorok (2012) method, the linear $Q$-function space and the random
forest-based approximation space were separately used. The method of Huang et
al. (2014) is not included in these simulations because Simoneau et al. (2019)
can be viewed as a doubly robust extension of their method. After estimation,
the values—either the expected truncated survival time or the survival
probability at $t=5$—are evaluated by generating an uncensored sample of size
$n_{eval}=10000$ according to the estimated policies. The simulations were
complete using R package dtrSurv. For reproducible research, the code is
available at https://github.com/hunyong/survQlearn.

###  4.2 Simulation results

In Figure 1 LEFT, the simulation results show how the truncated mean survival
time ($\phi_{\mu}$) of each policy behaves under different settings. For most
of the settings, on average, all estimated regimes have greater values than
the zero-order model, and the zero-order model has a higher value than the
observed policy. This implies that the standard of care improves the overall
survival time of patients over the observed policies and that individualized
treatment rules can further enhance the outcomes.

Figure 1 LEFT also shows that among the individualized treatment rules, the
proposed method outperforms the other methods in many settings. Even when the
method exhibits lower or about the same performance as the other methods for
the small sample size, it often has higher performance than the others for the
larger sample size. See, for example, the high censoring cases (the last row)
of either the observational or the RCT setting. In the settings where the
dimension of the covariate is high (the third row), the method of Simoneau et
al. (2019) often has high values, although sometimes estimation is not
possible. This might be because of the doubly robustness property of the
method and the fact that the true optimal decision rule might be close to a
linear function, as suggested by the linear interaction of the stage-wide
hazard functions in (6). The method of Simoneau et al. (2019), however, often
does not provide estimates when either the sample size of the terminal stage
is small due to censoring or failure or the dimension of covariates is large.

Figure 1 RIGHT shows the simulation results in terms of the survival
probability at $t=5$, or $\phi_{\sigma,5}$. Note that the estimated regimes
for the proposed method are distinct between LEFT and RIGHT but are the same
for the other methods. The proposed estimator shows better performance in
terms of the policy values under most settings, giving us similar
interpretations as the previous results.

Figure 1: The estimated mean truncated survival time (LEFT) and the estimated
survival probability at time $t=5$ (RIGHT) under each estimated policy. Each
dot represents one simulation replicate. The box plots are the quartile
summary of the points and the black dots are the average values of each
method. Each row of boxes represents a data generation scenarios, and each
column of boxes represents a combination of a study design and a sample size.

##  5 Leukemia data example

We analyze acute myeloid leukemia patient survival time using the proposed and
existing methods. We use a data set collected from a clinical trial in which
210 acute myeloid leukemia patients were randomized to a frontline treatment
followed by a salvage treatment that was adaptively chosen based on the
patient status (Wahed & Thall, 2013; Xu et al., 2016). One of four initial
treatment options was given to each patient with an equal probability:
combination of fludarabine, cytosine arabinoside, and idarubicin (FAI),
combination of FAI and all-trans retinoic acid (ATRA) (FAIA), combination of
FAI + granulocyte colony stimulating factor (G-CSF) (FAIG), and combination of
FAI + ATRA + G-CSF (FAIAG). Response to the first stage treatment took one of
the following forms: complete remission followed by disease progression,
failure (either partial remission or remission with subsequent relapse), or
resistance to treatment. Surviving patients that did not attain complete
remission were given a salvage treatment based on physician judgement. The
salvage treatments were labeled with one of two categories depending on
whether or not they contained high dose ara-C (HDAC) or not. We assume that
the salvage treatment effects are homogeneous across different versions and
that the current data set contains most of the necessary information to
achieve Assumption 9.

The data set is composed of age at entry, type of leukemia characterized by
cytogentic abnormality, treatment assignment, response to the initial
treatment (failure, complete remission, or resistance), duration of each
treatment stage, and survival and censoring status at each treatment stage.
Because both age and type of leukemia are known to have high power in
predicting overall survival time, they are the main factors taken into account
by clinicians when prescribing treatment in practice thereby mitigating the
potential violation of Assumption 9. Patient age at entry ranges from 21 to 87
years and is 61.1 years on average. The patients’ cytogenetics abnormalities
were grouped into three categories according to their prognosis (Wahed &
Thall, 2013): 65 poor, 79 intermediate, and 66 good. Out of 210 patients, 198
experienced death during the study, and 12 were censored. The overall observed
survival time ($T\wedge C$) was on average 414 days and had quartiles of 45,
184, and 443 days. We set $\tau=450$ days. Figure 2 overviews the dynamics of
the study.

Figure 2: Treatment stages and patient status with the number of instances in
the parentheses.

We derive the optimal treatment rules using the proposed method, the Goldberg
& Kosorok (2012) method, the Simoneau et al. (2019) method, and the zero-order
model. The value of the observed policy—the combination of randomization and
physician’s judgement—is estimated and is compared with the other rules. As in
Section 4, the Goldberg & Kosorok (2012) method was implemented using two
different Q-function approximations: a linear model and a random forest.
Because the Simoneau et al. (2019) method was implemented using R package
DTRreg the current version (1.7) of which allows binary treatment arms only,
the method was applied after dichotomizing the initial treatment according to
inclusion of all-trans retinoic acid.

We compare the performance of the five models based on cross validation, where
$K=100$ training data sets of size 168 (80%) are sampled to estimate the
treatment rules without replacement and stratified according to the initial
treatment, and the remaining test sets are held out for evaluation. The value
of each treatment rule is estimated by the inverse probability (censoring)
weighting approach as follows.

Define

$\displaystyle W_{i}(\bm{\pi})$
$\displaystyle=\frac{\prod_{q=1}^{2}1(\pi^{(q)}(\bm{H}_{i}^{(q)})=A_{i}^{(q)})\delta_{i}}{\widehat{\textnormal{pr}}(A^{(1)}_{i})\widehat{\textnormal{pr}}(A^{(2)}_{i}\mid\bm{H}^{(2)}_{i})\prod_{q=1}^{2}\widehat{\textnormal{pr}}(T^{(q)}_{i}>C^{(q)}_{i}\mid\text{age}_{i})},$
(9)

where the propensity of the salvage treatment is modeled as
$\textnormal{pr}(A^{(2)}=a)=(1+\exp(-\bm{H}^{(2)}\bm{\beta}))^{-1},$
$\bm{H}^{(2)}$ is a matrix with columns of one, $B^{(2)}$, $A^{(1)}$, age,
cytogenetic leukemia type, and response to initial treatment. Because there
was an insufficient number of censored patients in the study, the censoring
propensity score is modeled using only the patient age at entry. Then the
value of each policy $\bm{\pi}$ and each criterion
$\phi=\phi_{\mu},\phi_{\sigma,t}$ is estimated by

$\displaystyle\hat{\phi}_{\mu}(S^{\bm{\pi}})$
$\displaystyle=\frac{\sum_{i:\text{test
set}}(T_{i}\wedge\tau)W_{i}(\bm{\pi})}{\sum_{i:\text{test
set}}W_{i}(\bm{\pi})}$ and $\displaystyle\hat{\phi}_{\sigma,t}(S^{\bm{\pi}})$
$\displaystyle=\frac{\sum_{i:\text{test
set}}1(T_{i}>t)W_{i}(\bm{\pi})}{\sum_{i:\text{test set}}W_{i}(\bm{\pi})},$

where for the second criterion, we use $t=180$ to evaluate the six month
survival probabilities. We estimate the value of the observed policy, i.e., a
hybrid of randomized rule for the initial treatment and the physician
judgement for the salvage treatment, under the same evaluation procedure.

Figure 3 summarizes the value estimates of the estimated polices. Each box and
scatter plot corresponds to a policy, and within each plot, each dot
represents a cross-validation replicate of the estimated value (truncated mean
survival time (LEFT) or six-month survival probability (RIGHT)) of the
corresponding policy. The results show that the proposed method achieves
higher values than the other estimated policies on average in terms of both
the truncated mean survival time and the six month survival probability. This
is true when compared with the zero-order model implying the advantage of the
proposed individualized treatment rules over the standard of care. However,
the estimated methods do not have higher values than the observed policy value
on average. Recalling that the second treatment assignments were done
according to physician judgement, this implies that their judgement was fairly
successful and that a larger sample size and more patient
information—available only to physicians but not included in this data
set—could enhance the performance of the methods. Despite the gap, however,
the proposed method yields the closest values to the observed policy values
among the competing methods.

Figure 3: Estimated truncated mean survival time (LEFT) and estimated six-
month survival probability (RIGHT) under the dynamic treatment regimes. Each
dot represents one cross-validation replicate. The box plots are the quartile
summary of the points and the black dots are the average values of each
method.

##  6 Discussion

In this paper, we introduced a new dynamic treatment regime estimator for
survival outcomes. The estimator is versatile in the sense that it allows for
a flexible number of treatment stages and arms, it gives a non-linear decision
rule that maximizes either the mean survival time or the survival probability,
and it permits dependent censoring. In this section, we discuss some
considerations when using our method and how the method can be further
extended.

The proposed estimator assumes that the distribution is well characterized by
the treatment stages rather than the natural passage of time. That is, the
cohort at each treatment stage might include patients that received the stage
treatment on day 10 as well as patients that received the stage treatment on
day 1000. Thus, this method is effective when the disease dynamic is
considered stationary between stage transitions or when the treatment effect
contains relatively stronger signal than that of the baseline disease
dynamics.

Another consideration is the choice of the study length, $\tau$, which is
chosen so that the survival probabilities are reliably estimated up to $\tau$.
For the stability of estimation, the chosen $\tau$ should not be too large;
however, reducing $\tau$ results in information loss. If all patients are
censored at a certain stage, the survival probability is not estimable. In
this case, a smaller $\tau$ should be chosen so that every stage has a
sufficient number of observed failure times.

One interesting extension of the proposed method is optimizing a specific
quantile of the survival distribution such as the median survival time.
However, this task introduces a unique challenge with respect to the backward
recursion. Specifically, the optimal decisions made for later stages may not
be optimal once earlier decisions moderate the distribution. Thus, an
exhaustive search may be needed to find an optimal policy. An extension of the
quantile-optimal dynamic treatment regime estimators developed by Wang et al.
(2018) and Linn et al. (2017) to right-censored data would be interesting
future work.

The proposed method assumes an unrestricted policy class. In practice,
however, clinicians and patients may prefer understandable and simple rules. A
linear rule, for example, is often of interest. By posing a Cox-type
proportional hazards assumption on top of our method, the resulting policy
class becomes a set of linear functions. The task then reduces to replacing
the generalized random survival forest with the generalized Cox model, where
the estimated survival probabilities, instead of the observed time with
censoring indicators, are entered as the outcome variable.

In establishing the uniform consistency results, we use the fact that under
the given assumptions, the maximum size of the terminal node becomes
arbitrarily small in probability uniformly, as argued in Wager & Walther
(2015). However, the uniform rate of convergence can be very low. As a random-
split tree grows deeper and has an increasingly large number of terminal
nodes, there exists, with a not-very-small probability, at least one terminal
node that is created without enough splits on a certain covariate. To enhance
the rate, splitting can potentially be done in a stratified fashion, by
assigning a certain minimum number of splits to each variable. The effect of
the tree building rules on the rate of convergence would be an interesting
future study especially in the presence of noise variables in high-dimensional
data settings.

## Acknowledgement

The authors thank Peter F. Thall, Abdus S. Wahed, and Yanxun Xu for providing
the leukemia data, and Donglin Zeng for bringing the composite criterion
optimization into the discussion. This research was funded in part by grant
P01 CA142538 from the National Cancer Institute .

## Supplementary material

Supplementary Material available at Biometrika online includes the extended
proofs of Theorems 1 and 2.

## Appendix 1

### Summary of notation

See Table 3 for the summary of notation.

Table 3: Summary of notation | (Basic notation)
---|---
$q$ |  Treatment stage. $1,2,...,Q$
$\tau$ | Study length
$n^{(q)}$ |  Number of training examples available at stage $q$
| (Predictors)
$A^{(q)}$ |  Action at stage $q$.
$\bm{Z}^{(q)}$ |  Covariate information newly available at the beginning of stage $q$.
$\bm{H}^{(q)}$ |  History available at the beginning of stage $q$ including $\bm{Z}^{(q)}$.
$d^{(q)}$ |  Dimension of history $\bm{H}^{(q)}$ at stage $q$.
$\pi^{(q)}(\bm{H}^{(q)})$ |  Decision rule at stage $q$ given $\bm{H}^{(q)}$.
| (Stage-wise outcomes)
$T^{(q)}$ |  Failure time from the beginning of stage $q$.
$U^{(q)}$ |  Time to next treatment from the beginning of stage $q$.
$C^{(q)}$ |  Censoring time from the beginning of stage $q$.
$V^{(q)}$ |  Uncensored length of stage $q$. $V^{(q)}=T^{(q)}\wedge U^{(q)}$.
$X^{(q)}$ |  Observed length of stage $q$. $X^{(q)}=V^{(q)}\wedge C^{(q)}$.
$\delta^{(q)}$ |  Censoring indicator of stage $q$. $\delta^{(q)}=1(V^{(q)}\leq C^{(q)}).$
$\gamma^{(q)}$ |  Treatment indicator of stage $q$. $\gamma^{(q)}=1(T^{(q)}\leq U^{(q)}).$
$B^{(q)}$ |  Baseline time at stage $q$. $B^{(q)}=\sum_{k=1}^{q-1}X^{(k)},q>1$, $B^{(1)}=0$.
$L^{(q)}$ |  Remaining life at stage $q$. $L^{(q)}=T-B^{(q)}$.
| (Overall outcomes)
$T$ |  Overall failure time. $T=\sum_{q=1}^{Q}V^{(q)}$
$C$ | Overall censoring time
$\delta$ |  Overall censoring indicator. $\delta=1(T\leq C)$
| (Theoretical settings)
$L_{S},L_{G}$ |  Lipschitz continuity constants of Assumption01.
$k$ | Maximum terminal node size.
$\beta$ |  Terminal node size rate in Assumption03.
$\alpha$ |  Regular split constant in Assumption01
| (Survival functions and values)
$S$ | Overall or generic failure survival function.
$S^{(q)}$ |  Survival function of remaining life at the beginning of stage $q$.
$G$ | Overall or generic censoring survival function.
$\phi_{\mu}(S)$ |  Mean truncated survival time. $\phi_{\mu}(S)=E_{S}[T\wedge\tau].$
$\phi_{\sigma,t}(S)$ |  Survival probability at time $t$. $\phi_{\sigma,t}(S)=S(t).$
$\phi^{(q)}_{\mu}(S^{(q)})$ |  Mean truncated survival time, given $B^{(q)}$. $\phi^{(q)}_{\mu}(S^{(q)})=B^{(q)}+E_{S^{(q)}}[L^{(q)}\wedge(\tau-B^{(q)})].$
$\phi^{(q)}_{\sigma,t}(S^{(q)})$ |  Survival probability at time $t$, given $B^{(q)}$. $\phi^{(q)}_{\sigma,t}(S^{(q)})=S^{(q)}(t-B^{(q)}).$

## Appendix 2

### Proof of Theorem 1

The following Z-estimator lemma (Kosorok, 2007) plays a key role in this
consistency proof. Let $\Psi:\Theta\mapsto\mathbb{L}$ and
$\Psi_{n}:\Theta\mapsto\mathbb{L}$ be maps between two normed spaces, where
$\Theta$ is a parameter space, $\mathbb{L}$ is some normed space,
$\|\cdot\|_{\mathbb{L}}$ denotes the uniform norm, $\Psi$ is a fixed map, and
$\Psi_{n}$ is a data-dependent map.

######  Lemma 1 (Consistency of Z-estimators).

Let $\Psi(S_{0})=0$ for some $S_{0}\in\Theta$, and assume
$\|\Psi(S_{n})\|_{\mathbb{L}}\to 0$ implies $\|S_{n}-S_{0}\|_{\mathbb{L}}\to
0$ for any sequence $\\{S_{n}\\}\in\Theta$. Then, if
$\|\Psi_{n}(\hat{S}_{n})\|_{\mathbb{L}}\to 0$ in probability for some sequence
of estimators $\hat{S}_{n}\in\Theta$ and
$\sup_{S\in\Theta}\|\Psi_{n}(S)-\Psi(S)\|_{\mathbb{L}}\to 0$ in probability,
$\|\hat{S}_{n}-S_{0}\|_{\mathbb{L}}\to 0$ in probability, as $n\to\infty$.

For our specific problem, let $\Theta$ be the space of all covariate-
conditional survival functions. Define a normed space
$\mathbb{L}=D_{[-1,1]}\\{[0,\tau]\times\mathbb{R}^{d}\\}$, where $D_{A}B$ is
the space of all right-continuous left-limits functions with range $A$ and
support $B$.

Define

$\displaystyle\psi_{S,t}=$ $\displaystyle\delta S^{\prime}(t-X\mid
H^{\prime})+(1-\delta)\\{\frac{S(t\mid H)}{S(X\mid H)}\wedge 1\\}-S(t\mid H),$
(10) $\displaystyle\tilde{\psi}_{S,t,\tilde{S}^{\prime}}=$
$\displaystyle\delta\tilde{S}^{\prime}(t-X\mid
H^{\prime})+(1-\delta)\\{\frac{S(t\mid H)}{S(X\mid H)}\wedge 1\\}-S(t\mid H)$
(11) $\displaystyle\Psi(S)=$
$\displaystyle\frac{P\psi_{S,t}\delta_{h}}{P\delta_{h}}\equiv P_{\cdot\mid
h}\psi_{S,t}$ and (12) $\displaystyle\Psi_{n}(S)=$
$\displaystyle\frac{\mathbb{P}_{n}\tilde{\psi}_{S,t}k_{h}}{\mathbb{P}_{n}k_{h}}\equiv\mathbb{P}_{n,{\cdot\mid
k_{h}}}\tilde{\psi}_{S,t,\tilde{S}^{\prime}},$ (13)

where $P$ is the population average of function values, i.e., $Pf=\int
f(h)dP(h)$, $\mathbb{P}_{n}$ is the sample average of function values, i.e.,
$\mathbb{P}_{n}f=\int
f(h)d\mathbb{P}_{n}(h)=\frac{1}{n}\sum_{i=1}^{n}f(H_{i})$, $S^{\prime}$ is a
fixed survival probability that is distinct from $S$, $\tilde{S}^{\prime}$ is
a data-dependent version of $S^{\prime}$ such that
$\lim_{n\to\infty}\tilde{S}^{\prime}\to S^{\prime}$ in probability. In the
context of dynamic treatment regimes, $S^{\prime}$ is the survival probability
of the remaining life of the next stage and $H^{\prime}$ is the historical
information available at the beginning of the next stage. Also,
$\delta_{h}=I(\cdot=h)$ is the unnormalized Dirac measure, and $k_{h}$ the
unnormalized forest kernel. To be more specific,
$k_{h}=\frac{1}{n_{\text{tree}}}\sum_{b=1}^{n_{\text{tree}}}1(h\in L_{b}(h))$,
where $n_{\text{tree}}$ is the number of trees in the forest, $L_{b}(h)$ is
the terminal node of the $b$th tree of the forest that contains the point $h$.
We used the term ‘unnormalized’ to mean that they are not multiplied by the
sample (or the population) size.

Then by Lemma 1 and the following four propositions (Propositions 2 to 5),
$\|\hat{S}-S_{0}\|_{\mathbb{L}}\to 0$ in probability as $n\to\infty$.

######  Proposition 2.

$\Psi(S_{0})=0$ for $S_{0}\in\Theta$,

######  Proposition 3.

Assume Assumptions 1 to 4. As $n\to\infty$, $\|\Psi(\hat{S})\|_{\mathbb{L}}\to
0$ implies $\|\hat{S}-S_{0}\|_{\mathbb{L}}\to 0$,

######  Proposition 4.

Assume Assumptions 1 to 4. Let $\hat{S}$ be the kernel-conditional modified
Kaplan-Meier estimator, which is the modified Kaplan-Meier in (3) applied to
the data with weights indicated by the forest kernel, $k_{h}$. As
$n\to\infty$, $\|\Psi_{n}(\hat{S})\|_{\mathbb{L}}\to 0$ in probability.

######  Proposition 5.

Assume Assumptions 1 to 4. As $n\to\infty$,
$\sup_{S\in\Theta}\|\Psi_{n}(S)-\Psi(S)\|_{\mathbb{L}}\to 0$ in probability.

The proofs of Propositions 2–4 are provided in Supplementary Material. The
proof of Proposition 5 relies on the entropy calculation and the Donsker
theorem. The strategy is to show that the tree kernels and consequently the
forest kernels are Donsker, and to show that the $\psi$ and $\tilde{\psi}$
functions are also Donsker, which we will show in Proposition 6 below. Then by
empirical process theory and with Assumptions 1, 2, 3, and 4, the uniform
consistency is derived. Hence, we provide the related lemma (6) regrading the
entropy of the tree kernels before we move on to the proof of Proposition 5.
The proof of Propositions 6 and the rest of the proof of Proposition 5 are
given in Supplementary Material.

######  Proposition 6 (Bounded entropy integral of the tree and forest
kernels).

The collections of the unnormalized tree and forest kernel functions are
Donsker, where the tree kernels $k^{tree}(\cdot)$ are axis-aligned random
hyper-rectangles, and the forest kernels $k^{forest}(\cdot)$ are the mean of
arbitrarily many ($n_{\text{tree}}$) tree kernels.

###### of Corollary 1.

Corollary 1 follows from recursion once we show that the $Q$th stage random
survival forest is consistent and $n^{(Q)}>cn$ for some constant $c$ almost
surely. For each $\bm{h}\in\mathcal{H}^{(Q)}$ and treatment
$a\in\mathcal{A}^{(Q)}$, Assumption 7 guarantees that the probability of a
patient making all planned visits is at least $M^{Q-1}$. That is, $n^{(Q)}\geq
M^{(Q)}n$ almost surely. Uniform consistency of the stage-wise random survival
forest again relies on arbitrarily small node size,
$\sup_{h\in\mathcal{H}}\sup_{h^{\prime}\in k_{h}}\|h-h^{\prime}\|_{1}$.
Although the support in Assumption 6 is different from that in Assumption 2,
the bounded marginal density in Assumption 2 is sufficient to follow the same
lines of proof given in Meinshausen (2006) and Wager & Walther (2015). Now, by
having $S^{\prime}(t\mid\bm{h})=\tilde{S}^{\prime}(t\mid\bm{h})=1(t\leq 0)$
for all $\bm{h}$, Theorem 1 yields uniform consistency of
$\hat{S}^{(Q)}(\cdot\mid\cdot)$. Consistency of the rest of the sequence then
follows as $n^{(q)}\geq n^{(Q)}$ for $q=Q-1,Q-2,...,1$. ∎

## Appendix 3

### Proof of Theorem 2

###### of Theorem 2.

The proof is largely based on the lines of proof used by Murphy (2005) and
Goldberg & Kosorok (2012). We tailor their proofs based on the Q-learning
framework into the context of this paper. We introduce notation for some
intermediate stage quantities that are the analogs to the intermediate stage
Q- and V- functions. Let
$\Xi_{S,\bm{\pi}}^{(q)}(\bm{H}^{(q)},a^{(q)})=\phi^{(q)}\\{S^{(q)}_{\bm{\pi}}(\cdot\mid\bm{H}^{(q)},A^{(q)}=a^{(q)})\\}$
and
$\Phi_{S,\bm{\pi}}^{(q)}(\bm{H}^{(q)},a^{(q)})=\phi^{(q)}\\{S^{(q)}_{\bm{\pi}}(\cdot\mid\bm{H}^{(q)},A^{(q)}=\pi^{(q)}(\bm{H}^{(q)})\\}$.
The dependency of $\Xi$ and $\Phi$ on the choice of $\phi$ is marked with the
corresponding subscripts, e.g.,
$\Phi_{\cdot,\mu}^{(q)}(\cdot)=\phi_{\mu}^{(q)}(\cdot)$, or is suppressed, if
it is clear by context.

We modify Lemmas 1 and 2 of Murphy (2005) to obtain Lemmas 2 and 3 below. Then
we conclude the proof of Theorem 2 by bounding
$E\Bigg{[}\bigg{\\{}\Xi_{S_{*},\bm{\pi}_{*}}^{(q)}(\bm{H}^{(q)},$
$\bm{A}^{(q-1)},A^{(q)})-\hat{\Xi}^{(q)}(\bm{H}^{(q)},\bm{A}^{(q-1)},A^{(q)})\bigg{\\}}^{2}\Bigg{]}$
of Lemma 3 using a quantity

$Err_{\Xi^{(q+1)}}(\Xi^{(q)})=E\bigg{\\{}\max_{a^{(q+1)}}\Xi^{(q+1)}(\bm{H}^{(q+1)},A^{(q+1)})-\Xi^{(q)}(\bm{H}^{(q)},A^{(q)})\bigg{\\}}^{2},$

where a similar strategy was used in Murphy (2005) and Goldberg & Kosorok
(2012) and we drop the arguments of $\Xi^{(q)}$ and $\Xi^{(q+1)}$ as long as
no confusion arises.

$\displaystyle
Err_{\hat{\Xi}^{(q+1)}}(\hat{\Xi}^{(q)})-Err_{\hat{\Xi}^{(q+1)}}(\Xi^{(q)}_{*})$
$\displaystyle=E[\hat{\Xi}^{(q)}]^{2}-E[\Xi^{(q)}_{*}]^{2}+2E\bigg{[}(\Xi^{(q)}_{*}-\hat{\Xi}^{(q)})\max_{a^{(q+1)}}\Xi^{(q+1)}_{*}\bigg{]}$
$\displaystyle~{}~{}-2E\bigg{[}(\max_{a^{(q+1)}}\Xi^{(q+1)}_{*}-\max_{a^{(q+1)}}\hat{\Xi}^{(q+1)})(\Xi^{(q)}_{*}-\hat{\Xi}^{(q)})\bigg{]}$
$\displaystyle=E[\hat{\Xi}^{(q)}]^{2}-E[\Xi^{(q)}_{*}]^{2}+2E\bigg{[}(\Xi^{(q)}_{*}-\hat{\Xi}^{(q)})\Xi^{(q)}_{*}\bigg{]}$
$\displaystyle~{}~{}-2E\bigg{[}(\max_{a^{(q+1)}}\Xi^{(q+1)}_{*}-\max_{a^{(q+1)}}\hat{\Xi}^{(q+1)})(\Xi^{(q)}_{*}-\hat{\Xi}^{(q)})\bigg{]}$
$\displaystyle=E[\hat{\Xi}^{(q)}-\Xi^{(q)}_{*}]^{2}-2E\bigg{[}(\max_{a^{(q+1)}}\Xi^{(q+1)}_{*}-\max_{a^{(q+1)}}\hat{\Xi}^{(q+1)})(\Xi^{(q)}_{*}-\hat{\Xi}^{(q)})\bigg{]}$
$\displaystyle\geq
E[\hat{\Xi}^{(q)}-\Xi^{(q)}_{*}]^{2}-2E^{1/2}\bigg{[}\max_{a^{(q+1)}}\Xi^{(q+1)}_{*}-\max_{a^{(q+1)}}\hat{\Xi}^{(q+1)}\bigg{]}^{2}E^{1/2}\bigg{[}\Xi^{(q)}_{*}-\hat{\Xi}^{(q)}\bigg{]}^{2}$
$\displaystyle\geq
E[\hat{\Xi}^{(q)}-\Xi^{(q)}_{*}]^{2}-2\sqrt{LE\bigg{[}\Xi^{(q+1)}_{*}-\hat{\Xi}^{(q+1)}\bigg{]}^{2}}E^{1/2}\bigg{[}\Xi^{(q)}_{*}-\hat{\Xi}^{(q)}\bigg{]}^{2}$
$\displaystyle\geq
E[\hat{\Xi}^{(q)}-\Xi^{(q)}_{*}]^{2}-\frac{1}{2}\\{4LE\bigg{[}\Xi^{(q+1)}_{*}-\hat{\Xi}^{(q+1)}\bigg{]}^{2}+E\bigg{[}\Xi^{(q)}_{*}-\hat{\Xi}^{(q)}\bigg{]}^{2}\\}$
$\displaystyle=E[\hat{\Xi}^{(q)}-\Xi^{(q)}_{*}]^{2}-2LE[\hat{\Xi}^{(q+1)}-\Xi^{(q+1)}_{*}]^{2},$

where in the second equality,

$E\bigg{[}\max_{a^{(q+1)}}\Xi^{(q+1)}_{*}=\phi^{(q+1)}\\{S^{(q+1)}_{*}(\cdot\mid\bm{H}^{(q+1)},A^{(q+1)})\\}\mid\bm{H}^{(q)},A^{(q)}\bigg{]}=\Xi^{(q)}_{*}(\bm{H}^{(q)},A^{(q)})$

was used; in the first and the last inequalities, the Cauchy-Schwarz
inequality and that fact that $ab\leq(a^{2}+b^{2})/2$ were used; and the
second inequality holds true, since

$\displaystyle\\{\max_{a^{(q+1)}}\Xi^{(q+1)}_{*}(a^{(q+1)})-\max_{a^{(q+1)}}\hat{\Xi}^{(q+1)}(a^{(q+1)})\\}^{2}$
$\displaystyle\leq\max_{a^{(q+1)}}\\{\Xi^{(q+1)}_{*}(a^{(q+1)})-\hat{\Xi}^{(q+1)}(a^{(q+1)})\\}^{2}$
$\displaystyle\leq
LE\\{\Xi^{(q+1)}_{*}(a^{(q+1)})-\hat{\Xi}^{(q+1)}(a^{(q+1)})\\}^{2}.$

Then we have

$E[\hat{\Xi}^{(q)}-\Xi^{(q)}_{*}]^{2}\leq
2\sum_{q^{\prime}=q}^{Q}(4L)^{q^{\prime}-q}\big{\\{}Err_{\hat{\Xi}^{(q+1)}}(\hat{\Xi}^{(q)})-Err_{\hat{\Xi}^{(q+1)}}(\Xi^{(q)}_{*})\big{\\}},$

$q=1,2,...,Q.$ Thus, by Lemma 3, the error bound of the value function is
bounded by

$\displaystyle\Phi(\bm{\pi}_{*})-\Phi(\hat{\bm{\pi}})$ $\displaystyle\leq$
$\displaystyle\sum_{q=1}^{Q}2L^{q/2}\sqrt{2\sum_{q^{\prime}=q}^{Q}(4L)^{q^{\prime}-q}\big{\\{}Err_{\hat{\Xi}^{(q+1)}}(\hat{\Xi}^{(q)})-Err_{\hat{\Xi}^{(q+1)}}(\Xi^{(q)}_{*})\big{\\}}}$
$\displaystyle\leq$
$\displaystyle\sum_{q=1}^{Q}\sqrt{2^{3-2q}\sum_{q^{\prime}=q}^{Q}L^{q^{\prime}}\big{\\{}Err_{\hat{\Xi}^{(q+1)}}(\hat{\Xi}^{(q)})-Err_{\hat{\Xi}^{(q+1)}}(\Xi^{(q)}_{*})\big{\\}}}$
(14)

Now we obtain the bound of difference in $Err$ values:

$\displaystyle
Err_{\hat{\Xi}^{(q+1)}}(\hat{\Xi}^{(q)})-Err_{\hat{\Xi}^{(q+1)}}(\Xi^{(q)}_{*})$
$\displaystyle=E\bigg{[}\big{\\{}(\max_{a^{(q+1)}}\hat{\Xi}^{(q+1)}-\hat{\Xi}^{(q)})+(\max_{a^{(q+1)}}\hat{\Xi}^{(q+1)}-\Xi^{(q)}_{*})\big{\\}}\big{\\{}\Xi^{(q)}_{*}-\hat{\Xi}^{(q)}\big{\\}}\bigg{]}$
$\displaystyle\leq 2c(\phi)E\big{[}\Xi^{(q)}_{*}-\hat{\Xi}^{(q)}\big{]}$
$\displaystyle\leq
2c(\phi)\sup_{\bm{h}^{(q)},a^{(q)}}\big{\\{}\Xi^{(q)}_{*}(\bm{h}^{(q)},a^{(q)})-\hat{\Xi}^{(q)}(\bm{h}^{(q)},a^{(q)})\big{\\}}$
$\displaystyle\leq
2c(\phi)\sup_{\bm{h}^{(q)},a^{(q)}}\Bigg{[}\phi^{(q)}\bigg{\\{}\bigg{|}S^{(q)}_{*}(\cdot\mid\bm{h}^{(q)},a^{(q)})-\hat{S}^{(q)}(\cdot\mid\bm{h}^{(q)},a^{(q)})\bigg{|}\bigg{\\}}\Bigg{]}$
$\displaystyle\leq
2c(\phi)\sup_{\bm{h}^{(q)},a^{(q)}}\Bigg{[}\phi^{(q)}\bigg{\\{}\sup_{t\in[0,\tau]}\bigg{|}S^{(q)}_{*}(\cdot\mid\bm{h}^{(q)},a^{(q)})-\hat{S}^{(q)}(\cdot\mid\bm{h}^{(q)},a^{(q)})\bigg{|}\bigg{\\}}\Bigg{]}$
$\displaystyle\leq
2c(\phi)\phi^{(q)}\bigg{\\{}\sup_{\bm{h}^{(q)},a^{(q)},t\in[0,\tau]}\bigg{|}S^{(q)}_{*}(\cdot\mid\bm{h}^{(q)},a^{(q)})-\hat{S}^{(q)}(\cdot\mid\bm{h}^{(q)},a^{(q)})\bigg{|}\bigg{\\}},$

where $c(\phi_{\mu})=\tau$ and $c(\phi_{\sigma,t})=1$, and the first
inequality comes from the triangular inequality and the fact that any $\Xi$
function is bounded by $[0,\tau]$ for $\phi=\phi_{\mu}$ and $[0,1]$ for
$\phi=\phi_{\sigma,t}$, and, in the fourth inequality, Fatou’s lemma was used.

Given that all assumptions of Corollary 1 hold, and that the space of history
$\mathcal{H}$ incorporates the historical actions, the theorem gives us

$\sup_{t\in[0,\tau],h^{(q)}\in\mathcal{H}^{(q)},a^{(q)}\in\mathcal{A}^{(q)},q\in\\{1,2,...,Q\\}}|\hat{S}^{(q)}(t\mid
h^{(q)},a^{(q)})-S_{*}^{(q)}(t\mid h^{(q)},a^{(q)})|\to 0,$

in probability as $n\to\infty$.

Since $\phi$ is either integration over the interval $[0,\tau]$ for
$\phi=\phi_{\mu}$ or evaluation at $t$ for $\phi=\phi_{\sigma,t}$, the error
is bounded by $c(\phi)$ times the supremum error of the survival curves. Thus,

$Err_{\hat{\Xi}^{(q+1)}}(\hat{\Xi}^{(q)})-Err_{\hat{\Xi}^{(q+1)}}(\Xi^{(q)}_{*})\to
0$

in probability as $n\to\infty$ and since in (of Theorem 2.)
$\Phi(\bm{\pi}_{*})-\Phi(\hat{\bm{\pi}})$ is bounded by finite summations and
bounded transformations, the desired result holds. ∎

######  Lemma 2.

For treatment regimes $\tilde{\bm{\pi}}$ and $\bm{\pi}$,

$\displaystyle\Phi(\tilde{\bm{\pi}})-\Phi({\bm{\pi}})=-E_{\bm{\pi}}\bigg{\\{}\sum_{q=1}^{Q}\Delta_{S_{*},\tilde{\bm{\pi}}}^{(q)}(\bm{H}^{(q)},A^{(q)})\bigg{\\}},$

where
$\Delta_{S,\bm{\pi}}^{(q)}(\bm{h}^{(q)},a^{(q)})=\Xi_{S,\bm{\pi}}(\bm{h}^{(q)},a^{(q)})-\Phi_{S,\bm{\pi}}(\bm{h}^{(q)})$
is the advantage of treatment $a^{(q)}$ at stage $q$, and $S_{*}$ is the true
survival probability. This statement is true for both criteria
$\phi=\phi_{\mu},\phi_{\sigma,t}$.

######  Lemma 3.

For all functions, $\hat{\Xi}^{(q)}$ satisfying
$\hat{\pi}^{(q)}(\bm{h}^{(q)})\in\arg\max_{a^{(q)}}\hat{\Xi}^{(q)}(\bm{h}^{(q)},a^{(q)}),q=1,2,...,Q$,
we have,

$\displaystyle\Phi(\bm{\pi}_{*})-\Phi(\hat{\bm{\pi}})$ $\displaystyle\leq$
$\displaystyle\sum_{q=1}^{Q}2L^{q/2}\sqrt{E\Bigg{[}\bigg{\\{}\Xi_{S_{*},\bm{\pi}_{*}}^{(q)}(\bm{H}^{(q)},\bm{A}^{(q-1)},A^{(q)})-\hat{\Xi}^{(q)}(\bm{H}^{(q)},\bm{A}^{(q-1)},A^{(q)})\bigg{\\}}^{2}\Bigg{]}},$

where $\bm{\pi}_{*}=(\pi^{(1)},...,\pi^{(Q)})$ is the optimal decision rule
for the survival probability $S_{*}$. This statement is true for both criteria
$\phi=\phi_{\mu},\phi_{\sigma,t}$ .

## References

  * Breiman (2001) Breiman, L. (2001).  Random forests.  Machine learning 45, 5–32.
  * Breiman et al. (1984) Breiman, L., Friedman, J., Olshen, R. & Stone, C. (1984).  Classification and regression trees. wadsworth int.  Group 37, 237–251.
  * Cho et al. (2019) Cho, H., Jewell, N. P. & Kosorok, M. R. (2019).  Interval censored recursive forests.  arXiv preprint arXiv:1912.09983 .
  * Cole & Frangakis (2009) Cole, S. R. & Frangakis, C. E. (2009).  The consistency statement in causal inference: a definition or an assumption?  Epidemiology 20, 3–5.
  * Cui et al. (2017) Cui, Y., Zhu, R., Zhou, M. & Kosorok, M. (2017).  Consistency of survival tree and forest models: splitting bias and correction.  arXiv preprint arXiv:1707.09631 .
  * Edwards & Cambridge (2006) Edwards, J. C. & Cambridge, G. (2006).  B-cell targeting in rheumatoid arthritis and other autoimmune diseases.  Nature Reviews Immunology 6, 394–403.
  * Efron (1967) Efron, B. (1967).  The two sample problem with censored data.  In Proceedings of the fifth Berkeley symposium on mathematical statistics and probability, vol. 4.
  * Gehan (1965) Gehan, E. A. (1965).  A generalized wilcoxon test for comparing arbitrarily singly-censored samples.  Biometrika 52, 203–224.
  * Geurts et al. (2006) Geurts, P., Ernst, D. & Wehenkel, L. (2006).  Extremely randomized trees.  Machine learning 63, 3–42.
  * Goldberg & Kosorok (2012) Goldberg, Y. & Kosorok, M. R. (2012).  Q-learning with censored data.  Annals of statistics 40, 529.
  * Gordon & Olshen (1985) Gordon, L. & Olshen, R. A. (1985).  Tree-structured survival analysis.  Cancer treatment reports 69, 1065–1069.
  * Habermann et al. (2006) Habermann, T. M., Weller, E. A., Morrison, V. A., Gascoyne, R. D., Cassileth, P. A., Cohn, J. B., Dakhil, S. R., Woda, B., Fisher, R. I., Peterson, B. A. et al. (2006).  Rituximab-chop versus chop alone or with maintenance rituximab in older patients with diffuse large b-cell lymphoma.  Journal of clinical oncology 24, 3121–3127.
  * Hernán & Robins (2020) Hernán, M. A. & Robins, J. M. (2020).  Causal inference: What if.
  * Hogan & Radhakrishnan (2013) Hogan, J. & Radhakrishnan, J. (2013).  The treatment of minimal change disease in adults.  Journal of the American Society of Nephrology 24, 702–711.
  * Hothorn et al. (2005) Hothorn, T., Bühlmann, P., Dudoit, S., Molinaro, A. & Van Der Laan, M. J. (2005).  Survival ensembles.  Biostatistics 7, 355–373.
  * Huang et al. (2020) Huang, H.-J., Yang, L.-Y., Tung, H.-J., Ku, F.-C., Wu, R.-C., Tang, Y.-H., Chang, W.-Y., Jung, S.-M., Wang, C.-C., Lin, C.-T. et al. (2020).  Management and clinical outcomes of patients with recurrent/progressive ovarian clear cell carcinoma.  Journal of the Formosan Medical Association 119, 793–804.
  * Huang et al. (2014) Huang, X., Ning, J. & Wahed, A. S. (2014).  Optimization of individualized dynamic treatment regimes for recurrent diseases.  Statistics in medicine 33, 2363–2378.
  * Ishwaran et al. (2008) Ishwaran, H., Kogalur, U. B., Blackstone, E. H., Lauer, M. S. et al. (2008).  Random survival forests.  The annals of applied statistics 2, 841–860.
  * Jiang et al. (2017) Jiang, R., Lu, W., Song, R. & Davidian, M. (2017).  On estimation of optimal treatment regimes for maximizing t-year survival probability.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 79, 1165–1185.
  * Kidwell (2015) Kidwell, K. M. (2015).  Chapter 2: Dtrs and smarts: Definitions, designs, and applications.  In Adaptive Treatment Strategies in Practice: Planning Trials and Analyzing Data for Personalized Medicine. SIAM, pp. 7–23.
  * Kosorok (2007) Kosorok, M. R. (2007).  Introduction to empirical processes and semiparametric inference.  Springer Science & Business Media.
  * Kosorok & Laber (2019) Kosorok, M. R. & Laber, E. B. (2019).  Precision medicine.  Annual review of statistics and its application 6, 263–286.
  * Kosorok & Moodie (2015) Kosorok, M. R. & Moodie, E. E. (2015).  Adaptive treatment strategies in practice: planning trials and analyzing data for personalized medicine.  SIAM.
  * Linn et al. (2017) Linn, K. A., Laber, E. B. & Stefanski, L. A. (2017).  Interactive q-learning for quantiles.  Journal of the American Statistical Association 112, 638–649.
  * Meinshausen (2006) Meinshausen, N. (2006).  Quantile regression forests.  Journal of Machine Learning Research 7, 983–999.
  * Murphy (2003) Murphy, S. A. (2003).  Optimal dynamic treatment regimes.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 65, 331–355.
  * Murphy (2005) Murphy, S. A. (2005).  A generalization error for q-learning.  Journal of Machine Learning Research 6, 1073–1097.
  * Orellana et al. (2010) Orellana, L., Rotnitzky, A. & Robins, J. M. (2010).  Dynamic regime marginal structural mean models for estimation of optimal dynamic treatment regimes, part i: main content.  The international journal of biostatistics 6.
  * Peto & Peto (1972) Peto, R. & Peto, J. (1972).  Asymptotically efficient rank invariant test procedures.  Journal of the Royal Statistical Society: Series A (General) 135, 185–198.
  * Qian & Murphy (2011) Qian, M. & Murphy, S. A. (2011).  Performance guarantees for individualized treatment rules.  Annals of statistics 39, 1180.
  * Robertson & Uppuluri (1984) Robertson, J. B. & Uppuluri, V. (1984).  A generalized kaplan-meier estimator.  The Annals of Statistics , 366–371.
  * Robins & Rotnitzky (1992) Robins, J. M. & Rotnitzky, A. (1992).  Recovery of information and adjustment for dependent censoring using surrogate markers.  In AIDS epidemiology. Springer, pp. 297–331.
  * Robins et al. (1994) Robins, J. M., Rotnitzky, A. & Zhao, L. P. (1994).  Estimation of regression coefficients when some regressors are not always observed.  Journal of the American statistical Association 89, 846–866.
  * Rubin (2005) Rubin, D. B. (2005).  Bayesian inference for causal effects.  Handbook of statistics 25, 1–16.
  * Segal (1988) Segal, M. R. (1988).  Regression trees for censored data.  Biometrics , 35–47.
  * Simoneau et al. (2019) Simoneau, G., Moodie, E. E., Nijjar, J. S., Platt, R. W. & Investigators, S. E. R. A. I. C. (2019).  Estimating optimal dynamic treatment regimes with survival outcomes.  Journal of the American Statistical Association , 1–24.
  * Sutton & Barto (2018) Sutton, R. S. & Barto, A. G. (2018).  Reinforcement learning: An introduction.  The MIT Press.
  * Tarone & Ware (1977) Tarone, R. E. & Ware, J. (1977).  On distribution-free tests for equality of survival distributions.  Biometrika 64, 156–160.
  * Tsiatis (1975) Tsiatis, A. (1975).  A nonidentifiability aspect of the problem of competing risks.  Proceedings of the National Academy of Sciences 72, 20–22.
  * Tsiatis et al. (2019) Tsiatis, A. A., Davidian, M., T, H. S. & Laber, E. B. (2019).  Dynamic Treatment Regimes: Statistical Methods for Precision Medicine.  CRC press.
  * Wager & Athey (2018) Wager, S. & Athey, S. (2018).  Estimation and inference of heterogeneous treatment effects using random forests.  Journal of the American Statistical Association 113, 1228–1242.
  * Wager & Walther (2015) Wager, S. & Walther, G. (2015).  Adaptive concentration of regression trees, with application to random forests.  arXiv preprint arXiv:1503.06388 .
  * Wahed & Thall (2013) Wahed, A. S. & Thall, P. F. (2013).  Evaluating joint effects of induction–salvage treatment regimes on overall survival in acute leukaemia.  Journal of the Royal Statistical Society: Series C (Applied Statistics) 62, 67–83.
  * Wahed & Tsiatis (2006) Wahed, A. S. & Tsiatis, A. A. (2006).  Semiparametric efficient estimation of survival distributions in two-stage randomisation designs in clinical trials with censored data.  Biometrika 93, 163–177.
  * Wallace & Moodie (2015) Wallace, M. P. & Moodie, E. E. (2015).  Doubly-robust dynamic treatment regimen estimation via weighted least squares.  Biometrics 71, 636–644.
  * Wang et al. (2018) Wang, L., Zhou, Y., Song, R. & Sherwood, B. (2018).  Quantile-optimal treatment regimes.  Journal of the American Statistical Association 113, 1243–1254.
  * Watkins (1989) Watkins, C. J. C. H. (1989).  Learning from delayed rewards.  Ph.D. thesis.
  * Xu et al. (2016) Xu, Y., Müller, P., Wahed, A. S. & Thall, P. F. (2016).  Bayesian nonparametric estimation for dynamic treatment regimes with sequential transition times.  Journal of the American Statistical Association 111, 921–950.
  * Zhang et al. (2017) Zhang, Y., Laber, E. B., Davidian, M. & Tsiatis, A. A. (2017).  Estimation of optimal treatment regimes using lists.  Journal of the American Statistical Association,(just-accepted) .
  * Zhao et al. (2011) Zhao, Y., Zeng, D., Socinski, M. A. & Kosorok, M. R. (2011).  Reinforcement learning strategies for clinical trials in nonsmall cell lung cancer.  Biometrics 67, 1422–1433.
  * Zhu & Kosorok (2012) Zhu, R. & Kosorok, M. R. (2012).  Recursively imputed survival trees.  Journal of the American Statistical Association 107, 331–340.
  * Zhu et al. (2015) Zhu, R., Zeng, D. & Kosorok, M. R. (2015).  Reinforcement learning trees.  Journal of the American Statistical Association 110, 1770–1784.