# The Optimal Dynamic Treatment Rule SuperLearner: Considerations,
Performance, and Application

Lina Montoya Email address for correspondence<EMAIL_ADDRESS>University of
North Carolina, Chapel Hill, Department of Biostatistics  Mark van der Laan
University of California, Berkeley, Division of Biostatistics  Alexander
Luedtke  University of Washington, Department of Statistics  Jennifer Skeem
University of California, Berkeley, Departments of Social Welfare and Public
Policy  Jeremy Coyle  University of California, Berkeley, Division of
Biostatistics  Maya Petersen  University of California, Berkeley, Division of
Biostatistics

(January 2021)

###### Abstract

The optimal dynamic treatment rule (ODTR) framework offers an approach for
understanding which kinds of patients respond best to specific treatments – in
other words, treatment effect heterogeneity. Recently, there has been a
proliferation of methods for estimating the ODTR. One such method is an
extension of the SuperLearner algorithm – an ensemble method to optimally
combine candidate algorithms extensively used in prediction problems – to
ODTRs. Following the “causal roadmap,” we causally and statistically define
the ODTR and provide an introduction to estimating it using the ODTR
SuperLearner. Additionally, we highlight practical choices when implementing
the algorithm, including choice of candidate algorithms, metalearners to
combine the candidates, and risk functions to select the best combination of
algorithms. Using simulations, we illustrate how estimating the ODTR using
this SuperLearner approach can uncover treatment effect heterogeneity more
effectively than traditional approaches based on fitting a parametric
regression of the outcome on the treatment, covariates and treatment-covariate
interactions. We investigate the implications of choices in implementing an
ODTR SuperLearner at various sample sizes. Our results show the advantages of:
(1) including a combination of both flexible machine learning algorithms and
simple parametric estimators in the library of candidate algorithms; (2) using
an ensemble metalearner to combine candidates rather than selecting only the
best-performing candidate; (3) using the mean outcome under the rule as a risk
function. Finally, we apply the ODTR SuperLearner to the “Interventions”
study, an ongoing randomized controlled trial, to identify which justice-
involved adults with mental illness benefit most from cognitive behavioral
therapy (CBT) to reduce criminal re-offending.

##  1 Introduction

The primary objective of a clinical trial is often to evaluate the overall,
average effect of a treatment on an outcome in a given population [1, 2, 3].
To accomplish this objective in the point treatment setting, baseline
covariate, treatment, and outcome data are often collected and the average
treatment effect (ATE) is estimated, quantifying the average impact of the
treatment in a population. Researchers may then interpret the impact of the
treatment as beneficial, neutral, or harmful. In this interpretation, the
treatment’s impact is one-size-fits-all; in other words, the effect of the
treatment is interpreted as the same for everyone in the study population.
But, it may be the case that an intervention tends to yield better outcomes
for certain kinds of people but not for others. For example, because justice-
involved people with mental illness are a heterogeneous group with diverse
symptoms, risk factors, and other treatment-relevant characteristics [4, 5],
assigning Cognitive Behavioral Therapy (CBT) may decrease the probability of
recidivism for individuals with high risk of recidivism but not low risk of
recidivism [6]. The ATE analysis may lead one to conclude that there is no
treatment effect in a given population, when there is, in fact, a differential
treatment effect for levels of variables.

Precision health aims to shift the question from “which treatment works best”
to “which treatment works best _for whom_?” (sometimes, it further asks: at
what time? And/or at what dose? [3]). The point of moving towards this
question is to move towards better subject outcomes. While a range of novel
study designs can help to address these questions by generating data in which
individualized treatment effects are unconfounded [3, 7, 8], data from classic
randomized controlled trials also provide a rich data source for discovering
treatment effect heterogeneity. Under the assumption of no unmeasured
confounding, the same methods can be applied to observational data.

One way of learning which treatment works best for whom is to estimate effects
within subgroups. Following our above example within the field of criminal
justice, one could split the sample into subjects who are likely versus
unlikely to re-offend, and look at the average effect of CBT on recidivism
within these two risk categories. Such a classic subgroup analysis helps to
move a step closer to understanding the treatment that works best for whom.
However, the need to restrict the number of tests performed and to pre-specify
analyses limits traditional subgroup analyses to comparing intervention
effects in a small set of subgroups in which heterogeneous treatment effects
are expected [2, 9]. In practice, the subject characteristics that are most
important for determining the best-suited intervention may not be clear based
on background knowledge. Further, effectively predicting the type of
intervention that a subject will best respond to may require accounting for a
wide range of subject characteristics and complex interactions between them.
For instance, identifying the subjects most likely to respond to CBT versus,
for example, treatment as usual (TAU) may require considering not only risk
level, but also age, educational attainment, sex, substance abuse,
psychological distress, and internal motivation to adhere to treatment – as
well as various interactions between these. In summary, the challenge is to
take a wide range of subject characteristics and flexibly learn how to best
combine them into a strategy or rule that assigns to each subject the specific
intervention that works best for him or her.

Estimating the optimal dynamic treatment rule (ODTR) for a given population
offers a formal approach for learning about heterogeneous treatment effects
and developing such a strategy. A dynamic treatment rule can be thought of as
a rule or algorithm where the input is subject characteristics and the output
is an individualized treatment choice for each subject [10, 11, 12, 13]. An
optimal dynamic treatment rule (also known as an optimal treatment regime,
optimal strategy, individualized treatment rule, optimal policy, etc.) is the
dynamic treatment rule that yields the best overall subject outcomes [14, 15].
In our criminal justice example, a dynamic treatment rule takes as input
subject characteristics such as age, criminal history, and education level and
outputs a treatment decision – either CBT or TAU. The ODTR is the dynamic
treatment rule under which the highest proportion of patients are not re-
arrested. It is the most effective and, if one incorporates cost or
constraints on resources [16], efficient way of allocating the interventions
at our disposal based on measured subject characteristics.

There have been major advances in estimating the ODTR within the fields of
statistics and computer science, with important extensions to the case where
treatment decisions are made at multiple points in time. Regression-based
approaches, such as Q-learning, learn the ODTR by modeling the outcome
regression (i.e., the expected outcome given treatment and covariates)
directly [14, 17, 18, 19, 20]. Robins and Murphy developed methods of
estimating the ODTR by modeling blip-to-reference functions (i.e., the strata-
specific effect of the observed treatment versus control) and regret functions
(i.e., the strata-specific loss incurred when given the optimal treatment
versus the observed treatment), respectively [14, 15, 21]. Direct-estimation
approaches to learning the ODTR, such as outcome weighted learning (OWL), aim
to search among a large class of candidate rules for the one that yields the
best expected outcome [22, 23, 24]. These are examples of broad classes of
ODTR estimators; within and outside of them there has been a proliferation of
methods to estimate the ODTR (see [3, 25, 26] for reviews of the state of the
art in estimating ODTRs and precision medicine).

Given the vast number of methods available for estimating the ODTR, the
question becomes: which approach to use? In some settings, some algorithms may
work better than others. SuperLearning [27] (or, more specific to prediction,
stacked regression [28]) was originally proposed as a method for data-
adaptively choosing or combining prediction algorithms. The basic idea is to
define a library of candidate algorithms and choose the candidate or the
combination of candidates that gives the best performance based on V-fold
cross-validation. This requires defining: (1) the algorithms to include in the
library, (2) a parametric family of weighted combinations of these algorithms,
the “metalearning” step [29], and (3) the choice of performance metric (i.e.,
risk) as the criterion for selecting the optimal combination of algorithms.
Given these three requirements, then one can estimate the risk for each
combination of algorithms using V-fold cross-validation, and choose the
combination with the lowest cross-validated risk. The SuperLearner framework
has been implemented extensively for prediction problems [30, 31, 32, 33], and
has been extended to the ODTR setting [34, 35]. In particular, Luedtke and van
der Laan showed that in the randomized controlled trial (RCT) and sequential
multiple assignment randomized trial (SMART) [25, 36, 7] settings, under the
assumption that the loss function is bounded, the ODTR SuperLearner estimator
will be asymptotically equivalent to the ODTR estimator chosen by the oracle
selector (that is, the ODTR estimator, among the candidate ODTR estimators,
that yields the lowest risk under the true data distribution [27]). This
implies that the ODTR SuperLearner will asymptotically do as well as or better
than any single candidate estimator in the library, provided that none of the
candidate algorithms are correctly specified parametric models. If there is a
well-specified parametric model in the library, the ODTR SuperLearner
estimator of the ODTR will achieve near parametric rates of convergence to the
true rule.

These theoretical results lay important groundwork for understanding the
asymptotic benefits to using the algorithm; however, less has been published
on how the ODTR SuperLearner performs in finite samples, the practical
implications of key choices when implementing the algorithm, and illustrations
of implementing this algorithm on real RCT data. In this paper, we provide an
introduction to the implementation of the ODTR SuperLearner in the point
treatment setting, and use simulations to investigate the tradeoffs inherent
in these user-supplied choices and how they may differ with varying sample
sizes. In particular, for sample sizes 1,000 and 300, we examine: (1) how to
select the candidate algorithms for estimating the ODTR; specifically, the
costs and benefits to expanding the library to include a wider set of diverse
ODTR algorithms, including simple parametric models versus more data adaptive
algorithms, and blip-based versus direct estimation algorithms; (2)
implications of the choice of parametric family for creating weighted
combinations for candidate ODTR learners (i.e., choice of metalearner); and,
(3) implications of the choice of risk function used to judge performance and
thereby select the optimal weighted combination of candidate learners.
Finally, we apply the ODTR SuperLearner to real data generated from the
Correctional Intervention for People with Mental Illness, or “Interventions,”
trial, an ongoing RCT in which justice-involved adults with mental illness
were either randomized to CBT or TAU. In applying the ODTR SuperLearner to
this sample, we aim to identify which people benefit most from CBT versus TAU,
in order to reduce recidivism.

The organization of this article is as follows. First, we step through the
causal roadmap (as described in [37]) for defining the true ODTR for a given
population. We focus on the case in which baseline covariates are measured, a
single binary treatment is randomized, and an outcome is measured. We then
give a brief introduction to some estimators of the ODTR, and in particular,
describe the SuperLearner approach for estimating the optimal rule that builds
on Luedtke and van der Laan’s work [34]. We investigate the implications of
the three sets of implementation choices outlined above in finite samples
using simulations (with corresponding R code illustrating implementation of
all estimators considered), and the performance under such options. Lastly, we
show results for the ODTR SuperLearner algorithm applied to the
“Interventions” Study. We close with concluding remarks and future directions.

##  2 Causal Roadmap and ODTR Framework

###  2.1 Data and Causal Model

Consider point-treatment data where $W\in\mathcal{W}$ are baseline covariates,
$A\in\\{0,1\\}$ is the treatment, and $Y\in\mathbb{R}$ is the outcome measured
at the end of the study. Our data can be described by the following structural
causal model (SCM), $\mathcal{M}^{F}$ [38]:

$\displaystyle W$ $\displaystyle=f_{W}(U_{W})$ $\displaystyle A$
$\displaystyle=f_{A}(W,U_{A})$ $\displaystyle Y$
$\displaystyle=f_{Y}(W,A,U_{Y})\text{ ,}$

where the full data $X=(W,A,Y)$ are endogenous nodes,
$U=(U_{W},U_{A},U_{Y})\sim P_{U}$ are unmeasured exogenous variables, and
$f=(f_{W},f_{A},f_{Y})$ are structural equations. If it is known that data
were generated from an RCT using simple randomization with equal probability
to each arm, then the above structural causal model would state that $Y$ may
be affected by both $W$ and $A$, but that $W$ does not affect $A$ (as in the
“Interventions” trial); this can be represented in the above model by letting
$U_{A}\sim Bernoulli(p=0.5)$ and $A=U_{A}$. In this point treatment setting, a
dynamic treatment rule is a function $d$ that takes as input some function $V$
of the measured baseline covariates $W$ and outputs a treatment decision:
$V\rightarrow d(V)\in\\{0,1\\}$. For the remainder of the paper, we consider
the case where $V=W$; in other words, we consider treatment rules that
potentially respond to all measured baseline covariates. However,
consideration of dynamic rules based on a more restrictive set of baseline
covariates is also of frequent practical interest, allowing, for example, for
consideration of dynamic rules based on measurements that can be more readily
attained; all methods described extend directly to this case. We denote the
set of all dynamic treatment rules as $\mathcal{D}$.

The “Interventions” data consist of baseline covariates $W$, which include
intervention site, sex, ethnicity, age, Colorado Symptom Index (CSI) score (a
measure of psychiatric symptoms), level of substance use, Level of Service
Inventory (LSI) score (a risk score to predict future recidivism that
summarizes risk factors like criminal history, educational and employment
problems, and attitudes supportive of crime), number of prior adult
convictions, most serious offense, Treatment Motivation Questionnaire (TMQ)
score (a measure of internal motivation for undergoing treatment), and
substance use level; the randomized treatment $A$, either a manualized
Cognitive Behavioral Intervention for people criminal justice system
(abbreviated CBT; $A=1$) or treatment as usual (TAU), primarily psychiatric or
correctional services ($A=0$); and a binary outcome $Y$ of recidivism, an
indicator that the person was not re-arrested within one year after study
enrollment. In Table 2 we show the distribution of the data.

###  2.2 Target Causal Parameter

Let $d(W)$ be a deterministic function that takes as input a vector of
baseline covariates, and gives as output a treatment assignment (in this case,
either $0$ or $1$). For a given rule $d$, we intervene on the above SCM to
derive counterfactual outcomes:

$\displaystyle W$ $\displaystyle=f_{W}(U_{W})$ $\displaystyle A$
$\displaystyle=d(W)$ $\displaystyle Y_{d(W)}$
$\displaystyle=f_{Y}(W,d(W),U_{Y})\text{ .}$

Here, $Y_{d(W)}$ is the counterfactual outcome for a subject if his/her
treatment $A$ were assigned using the dynamic treatment rule $d(W)$; to
simplify notation we refer to this counterfactual outcome as $Y_{d}$. The
counterfactual outcomes for a person if he/she were assigned treatment or
given control are denoted $Y_{1}$ and $Y_{0}$, respectively. Together, the
distribution of the exogenous variables $P_{U}$ and structural equations $f$
imply a distribution of the counterfactual outcomes, and the SCM provides a
model for the set of possible counterfactual distributions:
$P_{U,X}\in\mathcal{M}^{F}$.

Our target parameter of interest in this paper is the ODTR, defined as the
rule that, among all candidate rules $\mathcal{D}$, yields the best expected
outcomes. Using the convention that larger values of $Y$ correspond to better
outcomes, an ODTR is defined as a maximizer of $E_{P_{U,X}}[Y_{d}]$ over all
candidate rules

$\displaystyle d^{*}$
$\displaystyle\in\operatorname*{arg\,max}_{d\in\mathcal{D}}E_{P_{U,X}}[Y_{d}]\text{
.}$ (1)

Any such ODTR can be defined in terms of the conditional additive treatment
effect (CATE), namely $E_{P_{U,X}}[Y_{1}-Y_{0}|W]$, which is the effect of
treatment for a given value of covariates $W$. Any ODTR assigns treatment 1
and 0 to all strata of covariates for which the CATE is positive and negative,
respectively. If the CATE is 0 for a particular $W$ (i.e., there is no
treatment effect for that strata of $W$), the ODTR as defined above may have
more than one maximizing rule and therefore may be non-unique; this is why the
RHS of equation 1 above is a set [39]. An ODTR can take an arbitrary value for
strata at which the CATE is 0. If we assume that assigning treatment 0 is
preferable to assigning treatment 1 in the absence of a treatment effect, then
we would prefer the following ODTR as a function of the CATE:

$d^{*}(W)\equiv\mathbb{I}\Big{[}E_{P_{U,X}}[Y_{1}-Y_{0}|W]>0\Big{]}\text{ .}$

In other words, if a subject’s expected counterfactual outcome is better under
treatment versus no treatment given his or her covariate profile, then assign
treatment; otherwise, assign control. A subject’s counterfactual outcome under
the ODTR is $Y_{d^{*}}$, and the expected outcome had everyone received the
treatment assigned by the ODTR is $E_{P_{U,X}}[Y_{d^{*}}]$.

Following our applied example, $Y_{1}$, $Y_{0}$, and $Y_{d}$ are the
counterfactual outcomes for a person if he/she were given CBT, TAU, and either
CBT or TAU based on the rule $d$, respectively; here, $d^{*}$ is the rule for
assigning CBT versus TAU using subjects’ covariates that would yield the
highest probability of no re-arrest, $E_{P_{U,X}}[Y_{d^{*}}]$.

###  2.3 Identification and Statistical Parameter

We assume that our observed data were generated by sampling $n$ independent
observations $O_{i}\equiv(W_{i},A_{i},Y_{i})$, $i=1,\ldots,n$, from a data
generating system described by $\mathcal{M}^{F}$ above (e.g., the
“Interventions” study consists of 441 i.i.d. observations of $O$). The
distribution of the observed data can be written as:

$P_{0}(O)=P_{W,0}(W)g_{0}(A|W)P_{Y,0}(Y|A,W)\text{ ,}$

where $P_{W,0}$ is the true distribution of $W$; $g_{0}$ is the true
conditional distribution of $A$ given $W$, or the treatment mechanism;
$P_{Y,0}$ is the true conditional distribution of $Y$ given $A$ and $W$. The
distribution of the data $P_{0}$ is an element of the statistical model
$\mathcal{M}$, which in our RCT example is semi-parametric. Further, if the
data are generated from an RCT design, as in the “Interventions” study, then
the true $g_{0}$ is known, and the backdoor criteria (with the implied
randomization assumption [38, 40]), $Y_{d}\perp A|W\hskip 14.22636pt\forall
d\in\mathcal{D}\text{ ,}$ and the positivity assumption,
$Pr\Big{(}\min_{a\in\\{0,1\\}}g_{0}(A=a|W)>0\Big{)}=1\text{ ,}$ hold by
design; in an observational data setting the randomization assumption requires
measurement of a sufficient set of baseline covariates, and the positivity
assumption may also pose greater challenges [41].

Define $Q(a,w)\equiv E[Y|A=a,W=w]$. Under the above assumption,
$E_{P_{U,X}}[Y_{d}]$ (a parameter of the counterfactual distribution) is
identified as $E_{0}[Q_{0}(A=d,W)]$ (a parameter of the observed distribution)
for any candidate rule $d$. Thus, the ODTR is identified by

$d_{0}^{*}\in\operatorname*{arg\,max}_{d\in\mathcal{D}}E_{0}[Q_{0}(A=d,W)]\text{
.}$

In addition, the CATE is identified as $Q_{0}(1,W)-Q_{0}(0,W)$, where the
latter is sometimes referred to as the blip function $B_{0}(W)$. Then, the
true optimal rule can also be defined as a parameter of the observed data
distribution using the blip function:

$d_{0}^{*}(W)\equiv\mathbb{I}[B_{0}(W)>0]\text{ .}$

Analogous to the definition of the ODTR as a function of the CATE, in words,
the blip function essentially says that if treatment for a type of subject
$W=w$ is effective (i.e., greater than 0), treat that type of person. If not,
do not treat him/her. If all subjects were assigned treatment in this way,
then this would result in the highest expected outcome, which is the goal.

##  3 Estimation of the ODTR

We denote estimators with a subscript $n$, so that, for example, an estimator
of the true ODTR $d^{*}_{0}$ is $d^{*}_{n}$. Estimates are functions of
$P_{n}$, which is the empirical distribution that gives each observation
weight $\frac{1}{n}$; $P_{n}\in\mathcal{M}_{NP}$, and $\mathcal{M}_{NP}$ is a
non-parametric model. In what follows, we briefly describe examples of common
methods for estimating the ODTR. We first describe methods that estimate the
ODTR via an estimate of the blip function. We then describe methods that
directly estimate a rule that maximizes the mean outcome.

###  3.1 Blip-based Approaches

Blip-based approaches aim to learn the blip, which implies an ODTR. A benefit
of doing this is that one can look at the distribution of the predicted
estimates of the blip for a given sample. Having the blip distribution allows
one to identify the patients in a sample who benefit most (or least, or
little) from treatment. Additionally, estimating the blip function can allow
for estimating the ODTR under resource constraints; for example, an ODTR in
which only $k$% of the population can receive treatment [16]. Below we
illustrate two methods of estimating the ODTR by way of the blip function
(i.e., blip-based estimators of the ODTR).

##### Single stage Q-learning

A plug-in estimator naturally follows from the above definition of the optimal
rule. One can estimate $Q_{n}(A,W)$ using any regression-based approach for
estimating an outcome regression and predict at $Q_{n}(1,W)$ and $Q_{n}(0,W)$.
This provides an estimate of the blip: $B_{n}(W)=Q_{n}(1,W)-Q_{n}(0,W)$, which
implies an estimate of the optimal rule: $d_{n}^{*}=\mathbb{I}[B_{n}(W)>0]$
[3, 18, 25, 42].

##### Estimating the blip function

Consider the double-robust pseudo-outcome [43]:

$D(Q,g)=\frac{2A-1}{g(A|W)}[Y-Q(A,W)]+Q(1,W)-Q(0,W)\text{ .}$

Importantly, $E_{0}[D(Q,g)|W]=B_{0}(W)$ if $Q=Q_{0}$ or $g=g_{0}$. Using this
result, one could estimate the blip by regressing the pseudo-outcome
$D_{n}(Q_{n},g_{n})$ (which we abbreviate from here on as $D_{n}$) on $W$
using any regression-based approach. As in the previous method, this estimate
of the blip implies an estimate of the optimal rule
$d_{n}^{*}=\mathbb{I}[B_{n}(W)>0]$ [15, 34, 44].

###  3.2 Direct Estimation Approaches for Maximizing the Expected Outcome

Instead of estimating the blip function, which implies an ODTR, one could
estimate the ODTR directly by selecting a rule $d$ that maximizes the
estimated $E_{U,X}[Y_{d}]$. Below we illustrate outcome weighted learning
(OWL) – one example of a direct-estimation method for the ODTR.

##### Single stage outcome weighted learning

We briefly describe the general concept of outcome weighted learning here, but
refer to zhao2012estimating [22] and rubin2012statistical [45] for a more
thorough explanation. The optimal rule defined above as a function of $P_{0}$
could equivalently be written as an inverse probability of treatment weighted
(IPTW) estimand:

$d_{0}^{*}\in\operatorname*{arg\,max}_{d\in\mathcal{D}}E_{0}[Q_{0}(A=d,W)]=\operatorname*{arg\,max}_{d\in\mathcal{D}}E_{0}\Bigg{[}\frac{Y}{g_{0}(A|W)}\mathbb{I}[A=d]\Bigg{]}\text{
.}$

Written this way, estimating $d_{0}^{*}$ could be regarded as a classification
problem, where the weighted outcome $\frac{Y}{g(A|W)}$ helps us learn what
kind of patients should get treatment: if a certain kind of patient $W=w$ has
large weighted outcomes and they were treated according to candidate rule $d$,
future patients with that covariate profile should be treated using that rule.
Conversely, the smaller the weighted outcome among patients $W$ who were
treated according to $d$, the larger the “misclassification error” and the
less likely those kinds of patients should be treated according to $d$. This
maximization problem is equivalent to the following minimization problem:

$\displaystyle d_{0}^{*}$
$\displaystyle\in\operatorname*{arg\,min}_{d\in\mathcal{D}}E_{0}\Bigg{[}\frac{Y}{g_{0}(A|W)}\mathbb{I}[A\neq
d]\Bigg{]}\text{ .}$ (2)

Now, if patients $W=w$ who did not follow the rule $d$ have large weighted
outcomes (and thus larger “misclassification error”), those kinds of patients
should be given the opposite treatment that $d$ proposes. Note that in the RCT
setting, if one uses the known $g_{0}$ and if treatments are given with equal
probability, then this reduces to finding the rule that minimizes the mean
outcome among patients who did not follow the rule. Equation 2 could
alternatively be written as a minimization problem for, instead of a rule $d$,
a function $f$:

$\displaystyle f_{0}^{*}$
$\displaystyle\in\operatorname*{arg\,min}_{f\in\mathcal{F}}E_{0}\Bigg{[}\frac{Y}{g_{0}(A|W)}I\Big{[}A\neq\frac{sign(f(W))+1}{2}\Big{]}\Bigg{]}\text{
,}$ (3)

where $sign(x)=-1$ if $x\leq 0$ and $sign(x)=1$ if $x>0$. Under the true data
distribution $P_{0}$, $f_{0}$ is the blip function, $B_{0}$. In order to solve
this minimization problem using data, we can use a plug-in estimator of (3);
however, since it is a 0-1 function (i.e., it is discontinuous and non-
convex), one could use a convex surrogate function to approximate it, to
instead minimize:

$\displaystyle f_{n}^{*}$
$\displaystyle\in\operatorname*{arg\,min}_{f\in\mathcal{F}}\frac{1}{n}\sum_{i=1}^{n}\frac{Y_{i}}{g_{n}(A_{i}|W_{i})}\Phi(A_{i}f(W_{i}))+\lambda_{n}\left\lVert
f\right\rVert^{2}\text{ ,}$ (4)

where $\Phi(t)$ is the surrogate loss function (e.g., hinge loss, exponential
loss, logistic loss), $\left\lVert f\right\rVert$ is the norm of $f$, and
$\lambda_{n}$ is the estimated penalization parameter on $f$ to avoid
overfitting of the rule. This can also be generalized with the IPTW function
replaced by the augmented IPTW [34, 45]. Once $f_{n}^{*}$ is found as the
solution to equation (4), the estimated ODTR is:

$d_{n}^{*}=sign(f_{n}^{*}(W))\text{ .}$

###  3.3 SuperLearner to Estimate ODTR

The overarching goal of SuperLearner is to let a user-supplied library of
candidate algorithms, such as specific implementations of the general
approaches described above, “team up” to improve estimation of the ODTR. In
order to implement the ODTR SuperLearner, there are three user-supplied
decisions one must make. First, one must consider the library of candidate
algorithms to include. These could include algorithms for estimating the blip
function (which imply a rule), algorithms that search for the ODTR directly
(such as OWL estimators), static rules that determine treatment regardless of
covariates, or combinations of the above classes of algorithms. Second, in
what is sometimes referred to as the metalearning step, one can either
implement a SuperLearner that chooses one algorithm out of the library of
candidate algorithms to include (i.e., “discrete” SuperLearner), or a
SuperLearner that is a combination the candidate algorithms (i.e.,
“continuous” SuperLearner). For the latter, one again has a choice of
metalearner; we consider weighted convex combinations of candidate estimators
of the blip and combinations of estimates of the rules themselves (through a
weighted “majority vote”). Finally, one must choose the risk function used to
judge the performance of the weighted combinations of algorithms (estimated
using V-fold cross validation). Here, we consider two risk functions: the
mean-squared error (MSE) and the mean outcome under the candidate rule.

The steps for implementing the ODTR SuperLearner are as follows; they closely
follow the implementation of the canonical SuperLearner for regression [46]:

  1. 1.

Choose $J$ candidate algorithms for estimating the optimal rule $d_{n,j}(W)$
for $j=1,...,J$. Candidates can include approaches based on estimating the
blip $B_{n,j}(W)$, e.g., candidate regressions of $D_{n}$ on $W$, where the
candidate regressions might consist of a parametric linear model
(corresponding to a classic approach of fitting a parametric outcome
regression on $A$ and $W$) as well as more flexible machine learning type
approaches such as neural networks [47], multivariate adaptive regression
splines [48], or recursive partitioning and regression trees [49]. Candidate
algorithms might also include approaches for estimating the optimal rule
directly, such an OWL estimator. Finally, the static treatment rules that
treat all patients or treat no patients, regardless of their covariate values,
can also be included as candidates. Inclusion of both simple parametric model
estimators, as well as static rules such as treating all and treating none, is
important to allow for the possibility that the underlying true ODTR may in
fact be simple (or well-approximated by a simple rule), and providing less
aggressive candidates in the SuperLearner library can help protect against
overfitting in finite samples.

  2. 2.

Split the data into $V$ exhaustive and mutually exclusive folds. Let each fold
in turn serve as the validation set and the complement data as the training
set.

  3. 3.

Fit each of the $J$ candidate algorithms on the training set. Importantly,
candidate algorithms might depend on nuisance parameters, and those nuisance
parameters should be fit on the training set, as well. For example, if a
candidate algorithm regresses $D_{n}$ on $W$ to estimate the blip (which
implies an ODTR), then $Q$ and $g$ should be fit and predicted on the training
set, and then plugged into $D$ to fit that candidate algorithm on the same
training set (this is also called “nested” cross-validation, described in
detail by [35]).

  4. 4.

Predict the estimated blip or the treatment assigned under the estimated ODTR
for each observation in the validation set for each algorithm, based on the
corresponding training set fit.

  5. 5.

Choose to either implement the discrete SuperLearner, which selects one
algorithm out of the candidate algorithms, or the continuous SuperLearner,
which creates a weighted average of the candidate algorithms.

    1. (a)

Continuous SuperLearner. Create different convex combinations of the candidate
blip or treatment rule estimates that were predicted on the validation set
(i.e., convex combinations of the predictions from the previous step).
Formally, define an estimator of $B_{n,\alpha}(W)$ or $d_{n,\alpha}(W)$ as a
convex combination of the candidate algorithms (indexed by $j$); each convex
combination of algorithms is indexed by a weight vector $\alpha$. A given
convex combination of blip estimates are denoted as:

$B_{n,\alpha}(W)=\sum_{j}\alpha_{j}B_{n,j}(W),\alpha_{j}\geq 0\forall
j,\sum_{j}\alpha_{j}=1\text{ .}$

Alternatively, the predicted treatments under the candidate ODTRs can be
combined as a weighted “majority vote” of the convex combination of the
candidate rules:

$d_{n,\alpha}(W)=\mathbb{I}\Big{[}\sum_{j}\alpha_{j}d_{n,j}(W)>\frac{1}{2}\Big{]},\alpha_{j}\geq
0\forall j,\sum_{j}\alpha_{j}=1\text{ .}$

    2. (b)

Discrete SuperLearner. The discrete SuperLearner, which chooses only one
candidate algorithm, can be thought of as a special case of the continuous
SuperLearner, where algorithms are still combined using a convex combination,
but each algorithm weight $\alpha_{j}$ must be either $0$ or $1$. Such an
approach may be particularly advantageous when sample size is small:

$B_{n,\alpha}(W)=\sum_{j}\alpha_{j}B_{n,j}(W),\alpha_{j}\in\\{0,1\\}\forall
j,\sum_{j}\alpha_{j}=1$

$d_{n,\alpha}(W)=\mathbb{I}\Big{[}\sum_{j}\alpha_{j}d_{n,j}(W)>\frac{1}{2}\Big{]},\alpha_{j}\in\\{0,1\\}\forall
j,\sum_{j}\alpha_{j}=1\text{ .}$

  6. 6.

Calculate an estimated risk within each validation set for each combination of
algorithms (i.e., for each convex combination indexed by a particular value
for $\alpha$). Here, we discuss two choices of risk functions for step (6)
above. First, mean-squared error risk targeting the blip function, which we
will refer to as $R_{MSE}$:

$R_{MSE}=E[(D(Q,g)-B(W))^{2}]\text{ .}$

Because the MSE compares $D_{n}$ to the predicted blip, the candidate
estimators under the MSE risk function must output estimated blip values. Of
note, $E_{P_{U,X}}[[Y_{1}-Y_{0}-B(W)]^{2}]$ is identified as $R_{MSE_{0}(B)}$
if either $Q=Q_{0}$ or $g=g_{0}$. The second risk function, which we call
$R_{E[Y_{d}]}$, uses the expected rule-specific outcome as criterion:

$R_{E[Y_{d}]}=-E[E[Y|A=d,W]]\text{ .}$

Intuitively, SuperLearner aims to choose the combination of treatment rule
algorithms that minimizes a cross-validated empirical risk, so it makes sense
to have that risk be the negative of the expected outcome, such that
SuperLearner chooses the combination of algorithms that maximizes the expected
outcome, since that is the ultimate goal of the ODTR. Candidate estimators for
the SuperLearner that use the expected mean outcome under the rule as the risk
function can include both blip estimators that imply treatment rules as well
as direct estimators of the treatment rules. When the expected rule specific
outcome is chosen as the risk function, a further practical choice is how to
estimate this quantity; we focus here on a cross-validated targeted maximum
likelihood estimator (TMLE) due to its favorable theoretical properties
(double robustness, semi-parametric efficiency, and greater protection against
overfitting through the use of sample splitting [46]); however, one can use
any estimator of treatment specific mean outcomes to estimate this quantity
[10, 50, 11].

  7. 7.

Average the risks across the validation sets resulting in one estimated cross
validated risk for each candidate convex combination (i.e., each possible
choice of $\alpha$).

  8. 8.

Choose the estimator (i.e., the convex combination $\alpha$) that yields the
smallest cross-validated empirical risk. Call this “best” weighting of the
algorithms $\alpha_{n}$.

  9. 9.

Fit each candidate estimator $B_{n,j}(W)$ of the blip or $d_{n,j}(W)$ of the
optimal rule on the entire data set. Generate predictions for each candidate
algorithm individually, and then combine them using the weights $\alpha_{n}$
obtained in the previous step. This is the SuperLearner estimate of the ODTR,
where $d_{n,B}^{*}=\mathbb{I}[B_{n,\alpha_{n}}(W)>0]$ or
$d_{n,d}^{*}=d_{n,\alpha_{n}}(W)$ directly.

We summarize the practical implications of 3 key choices for implementing ODTR
SuperLearner here and in Table 1. In the first dimension, selection of the
candidate algorithms, for illustration we consider having a library with only
blip function estimators (called “Blip only” library) or a library with blip
estimators, direct-estimation estimators, and static treatment rules that
treat everyone or no one (called “Full” library). If one chooses to include
direct-search estimators of the ODTR or static rules (i.e., functions that do
not output a blip estimate in the process), then one is constrained to using
the vote-based metalearner and $R_{E[Y_{d}]}$ risk function, because the blip-
based metalearner and $R_{MSE}$ risk function both rely on estimates of the
blip for combining and choosing the algorithms, respectively. For the second
dimension, the choice of how to combine algorithms, we consider either the
metalearner that (a) selects only one candidate algorithm (called “Discrete”),
(b) uses a weighted average to combine predicted blip estimates and directly
plugs those into the risk (called “Blip-based”), (c) uses a weighted average
to combine predicted treatments under the candidate combinations of rules and
creates a weighted majority vote of these candidate rules as input into the
risk (called “Vote-based”). The third dimension is the choice of performance
measure, risk function – either the MSE ($R_{MSE}$) or the mean outcome under
the candidate rule ($R_{E[Y_{d}]}$). For the second and third dimensions, if
one uses the vote-based metalearner, then the $R_{MSE}$ risk cannot be used
because $R_{MSE}$ requires an estimate of the blip to choose the best
algorithm, and the vote-based metalearner does not output an estimate of the
blip.

Choice 1: Library | Blip only | Full
---|---|---
Choice 2: Metalearner | Discrete | Continuous | Discrete | Continuous
Blip-based | Vote-based | Blip-based | Vote-based
Choice 3: Risk | $R_{MSE}$ | $R_{E[Y_{d}]}$ | $R_{MSE}$ | $R_{E[Y_{d}]}$ | $R_{MSE}$ | $R_{E[Y_{d}]}$ | $R_{MSE}$ | $R_{E[Y_{d}]}$ | $R_{MSE}$ | $R_{E[Y_{d}]}$ | $R_{MSE}$ | $R_{E[Y_{d}]}$
Possible? | ✓ | ✓ | ✓ | ✓ | ✗ | ✓ | ✗ | ✓ | ✗ | ✗ | ✗ | ✓

Table 1: Summary of the possible ODTR SuperLearner configurations across the
library, metalearner, and risk dimensions. The last row (“Possible?”)
indicates whether the particular configuration is possible to implement. The
checkmarks (✓) in the following table indicate that it is possible to
construct that kind of ODTR SuperLearner algorithm; the x-marks (✗) indicate
otherwise.

##  4 Simulation: Comparisons and Considerations of SuperLearner ODTR
Estimators

We use simulations to: (1) illustrate the potential benefit to estimating the
ODTR using a SuperLearner approach, as compared to a more traditional approach
to studying treatment effect heterogeneity based on fitting an outcome
regression with interaction terms on covariates and treatment, as is often
standard practice [1, 51, 52, 53]; and (2) investigate the implications of
practical choices when implementing an ODTR SuperLearner in finite samples,
including specification of candidate algorithms in the library, choice of
metalearner, and choice of risk function.

###  4.1 Data Generating Processes

All simulations were implemented in R [54], and the code, simulated data, and
results can be found at https://github.com/lmmontoya/SL.ODTR. We examine these
comparisons using two types of data generating processes (DGPs). Each
simulation consists of 1,000 iterations of either $n=1,000$ or $n=300$, to
assess the impacts of the different configurations as a function of sample
size. Both DGPs generate the covariates, treatment, and outcome as follows:

$\displaystyle W_{1},W_{2},W_{3},W_{4}$ $\displaystyle\sim
Normal(\mu=0,\sigma^{2}=1)$ $\displaystyle A$ $\displaystyle\sim
Bernoulli(p=0.5)$ $\displaystyle Y$ $\displaystyle\sim Bernoulli(p)\text{ .}$

The probability of having a successful outcome differs for the two DGPs,
which, in this case, means that the blip functions differ as well. The first
DGP is an example of one with a complex blip function, and the second DGP is
one with a blip function that is a simpler function of one variable. The first
DGP is directly from work by Luedtke and van der Laan [34, 16, 44], and the
second is modified from the first. The probability of a successful outcome for
DGP 1 is:

$\displaystyle p$
$\displaystyle=0.5logit^{-1}(1-W_{1}^{2}+3W_{2}+5W_{3}^{2}A-4.45A)+0.5logit^{-1}(-0.5-W_{3}+2W_{1}W_{2}+3|W_{2}|A-1.5A)\text{
,}$

then the true blip function is:

$\displaystyle B_{0}(W)=$ $\displaystyle
0.5[logit^{-1}(1-W_{1}^{2}+3W_{2}+5W_{3}^{2}-4.45)+logit^{-1}(-0.5-W_{3}+2W_{1}W_{2}+3|W_{2}|-1.5)$
$\displaystyle-
logit^{-1}(1-W_{1}^{2}+3W_{2})+logit^{-1}(-0.5-W_{3}+2W_{1}W_{2})]\text{ .}$

For DGP 1, $E_{P_{U,X}}[Y_{d^{*}}]\approx 0.5626$ and the true optimal
proportion treated $E_{P_{U,X}}[d^{*}]\approx 55.0\%$.
$E_{P_{U,X}}[Y_{1}]\approx 0.4638$ and $E_{P_{U,X}}[Y_{0}]\approx 0.4643$.

DGP 2’s probability of the outcome’s success is:

$\displaystyle p$ $\displaystyle=logit^{-1}(W_{1}+0.1A+W_{1}A)\text{ .}$

Thus the true blip function is:

$\displaystyle B_{0}(W)=$ $\displaystyle
logit^{-1}(W_{1}+0.1+W_{1})-logit^{-1}(W_{1})\text{ .}$

For DGP 2, $E_{P_{U,X}}[Y_{d^{*}}]\approx 0.5595$ and
$E_{P_{U,X}}[d^{*}]\approx 54.0\%$; $E_{P_{U,X}}[Y_{1}]\approx 0.5152$ and
$E_{P_{U,X}}[Y_{0}]\approx 0.5000$.

###  4.2 ODTR Estimators

For each data generating process, we consider a number of estimators of the
ODTR. First, mirroring epidemiologic practice, we model the outcome as an
additive function of the treatment and covariates, and interactions with the
treatment and all covariates) [1, 51, 53, 52]. Such an approach translates to
using the following parametric model for the outcome regression:

$h(E[Y|A,W])=\beta_{0}+\sum_{i=1}^{p}\beta_{i}W_{i}+\left(\gamma_{0}+\sum_{i=1}^{p}\gamma_{i}W_{i}\right)A\text{
,}$

where $h(.)$ denotes a link function, and $p$ is the number of baseline
covariates in $W$. Using a linear link, the following parametric model for the
blip function is implied:

$B(W)=\gamma_{0}+\sum_{i=1}^{p}\gamma_{i}W_{i}\text{ .}$

Next, we examine the finite sample implications of the aforementioned user-
supplied choices in implementing a SuperLearner estimator of the ODTR,
providing guidance for practical data analysis. First, we examine the choice
of library. We consider the library that only combines candidate blip
estimators (“Blip only” library; i.e., a library with candidate algorithms
suited for regressing $D_{n}$ on $W$) versus a library that has blip
estimators, direct-estimation algorithms, and static treatment rules (“Full”
library). The “Blip only” libraries consist of either:

  1. (a)

Simple parametric models only (denoted “Parametric blip models”). This
consisted of univariate GLMs with each covariate.

  2. (b)

Machine learning algorithms only (denoted “ML blip models”), such as SL.glm
(generalized linear models), SL.mean (the average), SL.glm.interaction
(generalized linear models with interactions between all pairs of variables),
SL.earth (multivariate adaptive regression splines [48]), SL.nnet (neural
networks [47]), SL.svm (support vector machines [55]), and SL.rpart (recursive
partitioning and regression trees [49]) from the SuperLearner package [56]

  3. (c)

A combination of (a) and (b) above, denoted “Parametric + ML blip models”

The “Full” library includes other ODTR algorithms like direct-estimation
methods, static rules, and other blip-based methods. Specifically, the “Full”
library includes either the “ML blip models” or “Parametric + ML blip models”
from the “Blip only” libraries above, in addition to Q-learning [25], OWL
[22], residual weighted learning (RWL) [57], efficient augmentation and
relaxation learning (EARL) [58], optimal classification algorithms [59] (the
latter 4 are from the DynTxRegime package [60], with function names owl, rwl,
earl, and optimalClass, respectively), and static rules that treat all
patients and no patients, regardless of the patient covariate profiles. For
the algorithms from the DynTxRegime package, except for nuisance parameters
$Q_{n}$ and $g_{n}$, we use default parameters, and the rule as a function of
all covariates. Additionally, for the optimal classification algorithm, the
solver method is recursive partitioning for regression trees (rpart). Thus,
the possible “Full” libraries are:

  1. (d)

Algorithms from the “ML blip models” library, plus direct maximizers and
static rules, denoted “ML blip models and $E[Y_{d}]$ maximizers”

  2. (e)

All possible algorithms – that is, algorithms from the “Parametric + ML blip
models” library, plus direct maximizers and static rules, denoted “All blip
models and $E[Y_{d}]$ maximizers”

Second, we examine the performance of different metalearners for combining the
candidate ODTR algorithms. We examine the blip-based metalearner using the
“Blip only” libraries, and the discrete and vote-based metalearners using both
the “Blip only” libraries and “Full” libraries.

Third, we examine the performance of using either the MSE $R_{MSE}$ or the
expected outcome under the candidate rule $R_{E[Y_{d}]}$ as risk criteria for
choosing the optimal linear combination of candidate ODTR algorithms. In
particular, CV-TMLE is used for estimating $R_{E[Y_{d}]}$.

We fully estimate the ODTR SuperLearner by additionally estimating nuisance
parameters (as opposed to using the true nuisance parameter functions) in a
nested fashion [35] as described above. Specifically, we estimate $Q_{n}$ and
$g_{n}$ using the canonical SuperLearner [27, 56] and a correctly specified
parametric model (i.e., a logistic regression of $A$ on the intercept),
respectively. We use 10-fold cross-validation throughout.

###  4.3 Performance Metrics

We measure performance by computing the percent accuracy of the algorithm;
that is, in a sample, the proportion of times the treatment assigned by the
estimated ODTR matches the true optimal treatment (i.e., the treatment that
would have been assigned under the true ODTR) for each observation in the
sample, averaged across simulation repetitions. We also evaluate performance
metrics of the difference between the true conditional expected outcome under
the estimated rule, averaged across the sample, compared to the true mean
outcome under the true optimal rule
$E_{n}[Q_{0}(Y|A=d_{n}^{*},W)]-E_{0}[Y_{d_{0}^{*}}]$ (as an approximation to
the regret $E_{0}[Q_{0}(Y|A=d_{n}^{*},W)]-E_{0}[Y_{d_{0}^{*}}]$). We compute
the mean and variance of this difference across the simulation repetitions.
Instead of presenting the raw variance of the regret, we present a variance
relative to the regret yielded by estimating the blip, and thus the optimal
rule, using as a parametric GLM that models the blip as an additive function
of all covariates. Additionally, we compute $2.5^{th}$ and $97.5^{th}$
quantiles of the distribution of $E_{n}[Q_{0}(Y|A=d_{n}^{*},W)]$ across the
simulation repetitions.

###  4.4 Simulation Results

Figure 1 displays simulation results. Below we discuss results specific to
each DGP, configuration dimension, and sample size. In general, results within
each DGP (i.e., across sample sizes) follow generally similar patterns;
however, any differences in performance between libraries, metalearners, or
risks are more pronounced for the smaller sample size.

####  4.4.1 DGP 1 Results - “Complex” Blip Function

Above, we showed that DGP 1 yields a blip function that is a complex function
of all of the available covariates. Here, for a larger sample size, we would
expect a benefit to more aggressive approaches to estimating the ODTR, such as
including more flexible machine learning-based approaches in the library of
candidates, as well as use of more aggressive metalearners (either vote- or
blip-based) over a discrete SuperLearner due to the better ability of these
approaches to approximate the true underlying function. That said, for smaller
sample sizes, this benefit might be attenuated, or even result in worse
performance than simpler alternatives. For this DGP, at sample size of 1,000,
indeed we find a benefit to the use of both more aggressive metalearners and
larger libraries. Interestingly, however, this benefit is maintained for
sample size 300. Specifically, libraries that included data adaptive, machine
learning algorithms (as opposed to incorrectly specified parametric models
alone) more accurately and precisely approximated the rule, even for sample
size of 300. Results also show that for both sample sizes, within the discrete
metalearner, the $R_{E[Y_{d}]}$ risk performed better than $R_{MSE}$ risk, and
more saliently, the blip-based and vote-based metalearners performed better
than the discrete SuperLearner. Finally, as predicted by theory, all
SuperLearner approaches evaluated substantially outperformed a traditional
parametric regression approach at both sample sizes. Below we describe results
specific to each sample size.

#####  $n$ = 1,000

Libraries that contain machine learning algorithms (i.e., “ML blip models,”
“Parametric + ML blip models,” “ML blip models and EYd maximizers, ” and “Blip
models and EYd maximizers”) overall perform better than libraries with
parametric models only (i.e., “Parametric blip models”) and the standard GLM
(i.e., “GLM”), across all performance metrics. For example, the percent match
between the true ODTR and the estimated ODTR spans from 72.0%-77.7% for any
libraries with machine learning algorithms, whereas the percent match for
libraries with only parametric models is from 56.4% to 58.0%.

There are no stark differences within the libraries that contain machine
learning algorithms across the metalearner and risk dimensions, except when
using a discrete metalearner and $R_{MSE}$ risk. Specifically, the discrete
metalearner that uses $R_{MSE}$ has a higher average regret and relative
variance than all other algorithms that use machine learning (e.g., for the
“Parametric + ML blip models” “Blip only” library that uses a discrete
metalearner, the average regret when using $R_{MSE}$ versus $R_{E[Y_{d}]}$ is
-0.0389 versus -0.0284, respectively, and the relative variance when using
$R_{MSE}$ versus $R_{E[Y_{d}]}$ is 2.137 versus 0.7781, respectively).

#####  $n$ = 300

As expected, given the limited data available to estimate a complex underlying
function, both accuracy of treatment assignment and approximated regret (the
extent to which the expected outcome under the estimated rule fell short of
the best outcomes available) deteriorated with smaller sample sizes. That
said, even in this challenging situation of a complex true pattern of
treatment effect heterogeneity and limited data with which to discover it, the
ODTR SuperLearner would have improved the expected outcome by approximately
4.5% relative to the static rule that treats everyone, an approach that would
have been suggested based on estimation of the ATE.

Libraries with only parametric models perform worse than libraries that
contain machine learning algorithms in terms of average regret and accuracy.
For example, SuperLearner ODTRs that contain libraries with parametric blip
models match 54%-55.4% of the time with the true ODTR, while the SuperLearners
that contain libraries with machine learning algorithms match 60.9%-66.1% of
the time. These results parallel those found with sample size 1,000, except
the discrepancy between libraries with machine learning models and parametric
models is not as pronounced.

Similar to the $n=1,000$ case, among the libraries that used machine learning
algorithms, using the performance of the rule as risk is better across all
performance metrics than using MSE as risk for the discrete metalearner. As
long as machine learning methods were included in the library, performance was
similar across risk functions and choice of metalearners, with the exception
of the MSE risk combined with the discrete metalearner.

####  4.4.2 DGP 2 Results - “Simple” Blip Function

As shown above, DGP 2 has a true blip function that is a simple function of
one covariate (referred to here as a “simple” blip function). Here, the true
optimal rule is described by a simple parametric model for the blip; thus, we
expect this approach to perform well. However, in practice one is unlikely to
be sure that the truth can be well approximated by a simple rule; it is thus
of interest to evaluate what price is paid for expanding the library to
include more aggressive machine learning algorithms and metalearners. In
particular, one might expect that, for smaller sample size, adding machine
learning-based candidates and more complex metalearners risks substantial
drop-off in performance. However, specifying a library that includes simpler
parametric models, in addition to machine learning approaches, may help
mitigate this risk. In fact, for this particular DGP, we see, across
metalearners and risks, only a small price in performance from adding machine
learning algorithms to a library including only simple parametric models. In
short, having an ODTR SuperLearner library that also includes parametric
models is better than having a library that only consists of data adaptive,
machine learning algorithms. Within the libraries that did contain parametric
models, particularly for the discrete metalearner, $R_{MSE}$ risk performs
slightly better than $R_{E[Y_{d}]}$ risk; for other metalearners there is
little difference in performance in terms of risk. Performance of the
metalearners varies slightly by sample size. Below we describe results
specific to each sample size.

#####  $n$ = 1,000

In terms of accuracy, the libraries that only contain parametric models
perform the best, followed closely by libraries that contain parametric models
and machine learning models, followed by the library with only machine
learning models. This pattern is evident in the percent match with the true
ODTR: for example, within the discrete metalearner with $R_{MSE}$ risk, the
percent accuracy is 90.7% for the library with parametric models only
(“Parametric blip models” library), 88.8% for the library that contain both
parametric models and machine learning models (“Parametric + ML blip models”),
and 81.9% for the library that contains machine learning algorithms only (“ML
blip models”). This same pattern is apparent in terms of average regret; that
is, the libraries that contain parametric blip models or a combination of
parametric blip models and machine learning models have the lowest average
regret (-0.0041 to -0.0095), while the libraries that only contain machine
learning models have the highest average regret (-0.0100 to -0.0138). Modeling
the blip with a single, parametric model often used in standard
epidemiological analysis (which, in this case, is incorrectly specified)
yields an average regret of -0.0100 (higher than the libraries with a
combination of parametric models and/or machine learning algorithms, and at
the same level as having machine learning algorithms only).

Within the libraries that contain parametric models and use a discrete
metalearner, $R_{MSE}$ performs better than $R_{E[Y_{d}]}$. For example, the
mean regret and relative variance for the discrete metalearner that only used
parametric models in the library is -0.0041 and 1.0267, respectively, when
using $R_{MSE}$ risk, and -.0046 and 1.3174, respectively, when using
$R_{E[Y_{d}]}$ risk. Otherwise, there were no apparent differences in
performance by risk.

For libraries that contain parametric models and use $R_{MSE}$, the discrete
ODTR SuperLearner performs better than the blip-based ODTR SuperLearner, with
regards to all performance metrics. For example, the average regret and
relative variance for the library with only parametric models that uses
$R_{MSE}$ was -0.0041 and 1.0267, respectively, when using a discrete
metalearner versus -.0059 and 1.1855, respectively, when using the blip-based
metalearner. This pattern is also evident for the library that has both
parametric models and machine learning algorithms.

#####  $n$ = 300

As in the case where $n=1,000$, the library with only parametric models
performs best in terms of accuracy, followed by libraries with parametric
models and machine learning models, and finally libraries with only machine
learning algorithms; again, however, DGP 1 illustrates the risks of such a
strategy. Moreover, even at this small sample size, there is only a small
price to pay for adding machine learning-based learners to a library that also
includes simple parametric models. For example, for the discrete metalearner
that uses $R_{MSE}$, the percent accuracy for the library that uses only
parametric models is 78%, followed by a 75.5% accuracy when there is a
combination of parametric models and machine learning, while the library with
only machine learning models had a 61.7% match with the true ODTR. While this
pattern parallels that of the $n=1,000$ case, the dropoff in accuracy when the
library uses only parametric models versus when the library only uses machine
learning algorithms is larger in terms of accuracy in the smaller sample size
(16.3% difference) versus the larger sample size (8.8% difference).

Similar to the $n=1,000$ case, among libraries that contain parametric models
and in the discrete metalearner case, $R_{MSE}$ yields slightly better
performance results than $R_{E[Y_{d}]}$. In contrast to the $n=1,000$ case,
for libraries that contain parametric models and used $R_{MSE}$, the blip-
based metalearner performs slightly better than the discrete metalearner, with
regards to all performance metrics. For example, the average regret and
relative variance for the library with only parametric models that use
$R_{MSE}$ is -0.0188 and 1.6102, respectively, when using a blip-based
metalearner versus -0.0190 and 1.8109, respectively, when using the discrete
metalearner.

##  5 Application of ODTR SuperLearner to “Interventions” Study

In the “Interventions” study, 231 (52.2%) participants received CBT and 210
(47.8%) TAU. Out of the 441 participants, 271 (61.5%) were not re-arrested
within the year. The estimated probability of no re-arrest had everyone been
assigned CBT is 62.2%, and the estimated probability of no-arrest had everyone
been assigned TAU is 60.7%; there was no significant difference between these
two probabilities (risk difference: 1.51%, CI: [-8.03%,11.06%]). After
adjusting for covariates using TMLE to improve the precision on this ATE
estimate [61], the risk difference is, similarly, 1.53% (CI: [-7.31%,
10.37%]).

Figure 2 shows subgroup plots for each covariate – that is, the unadjusted
difference in probability of no re-arrest between those who received CBT
versus TAU, within each covariate group level. One might begin to identify
trends towards differential treatment effects; for example, participants may
have benefited more from CBT at the San Francisco site, or if they had
offended three or more times. Accurate interpretation of any such subgroup
analyses, however, requires variance estimates and hypothesis tests with
appropriate correction for multiple testing. In addition, as mentioned before,
it may be the case that the best way to assign treatment is by using
information on more than one variable at a time, and even interactions between
those variables.

Thus, we estimated the ODTR on the “Interventions” data to determine which
justice-involved adults with mental illness should receive CBT. Specifically,
we implemented the ODTR SuperLearner with a blip-only library, a continuous,
blip-based metalearner, and $R_{E[Y_{d}]}$ risk function. We chose a blip-
based library in order to generate estimates of the blip, which themselves can
be informative. The library for $d^{*}_{n}$ consisted of a combination of
simple parametric models (univariate GLMs with each covariate) and machine
learning algorithms (SL.glm, SL.mean, SL.glm.interaction, SL.earth, and
SL.rpart). As in the simulations, the outcome regression $Q_{n}$ was estimated
using the canonical SuperLearner, $g_{n}$ was estimated as an intercept-only
logistic regression, and we used 10-fold cross validation.

Interestingly, despite implementing a continuous metalearner, the ODTR
SuperLearner algorithm assigned all weight on a GLM that modeled the blip as a
linear function of only substance use. As shown in Figure 3, a plot depicting
the distribution of the predicted blip estimates for all participants, the
algorithm can be interpreted as such: if a justice-involved person with mental
illness has a low substance use score, give him/her CBT; otherwise, give
him/her TAU. Under this ODTR estimate, for this sample, 52.38% of the
participants would receive CBT.

##  6 Conclusions

We described the ODTR SuperLearner and illustrated its performance for sample
DGPs under different configurations of the algorithm and finite sample sizes.
These results build on existing work [34, 35] by fully estimating the ODTR and
including an expanded SuperLearner library with not only blip-based regression
estimators, but also direct-estimation methods and static interventions. We
highlighted the practical choices one must consider when implementing the ODTR
SuperLearner across three dimensions: (1) the ODTR SuperLearner library of
candidate algorithms, namely, whether to include parametric models, machine
learning algorithms, or both; whether to include only estimators that output a
predicted blip or include a combination of blip estimators, direct estimators,
and static treatment rules, (2) the metalearner that either chooses a single
algorithm or combines the algorithms and (3) the risk function that chooses
the “best” estimator or combination of estimators of the candidate ODTR
algorithms.

Simulation-based results illustrate the shortcomings of an approach to
treatment effect heterogeneity based on approximating the blip as an additive
function of the available covariates, or equivalently, modeling the outcome as
an additive function of the treatment and covariates, and interactions between
the treatment and all covariates, which is common practice in epidemiologic
analyses for heterogeneous treatment effects [1, 51, 53, 52]. With respect to
choice of library, we recommend specifying a library with both simple
parametric models and more aggressive data adaptive algorithms, as well as
static rules such as treat all or treat none, allowing for flexible estimation
of both simple and complex underlying rules. Inclusion of a full range of
algorithms from simple to aggressive was particularly important for small
sample sizes. In terms of the choice of metalearner, both vote- and blip-based
ensemble learners performed well; a vote-based metalearner has the advantage,
however, of allowing for the integration of a larger library of candidate
algorithms (including direct estimation approaches) and ease of integration of
static rules. Of note, in these simulations, vote- and blip-based metalearners
outperformed the discrete ODTR SuperLearner approach, even for sample size
300. However, we caution that this may not always be the case and when sample
size is small, a discrete SuperLearner approach may provide benefits – in
fact, one could include a convex metalearner as a candidate algorithm.
Finally, with respect to choice of risk function, both MSE and the expected
outcome under the rule performed well; in practice one might prefer
$R_{E[Y_{d}]}$ because it allows for the use of a larger library of candidate
algorithms.

Additionally, as an illustration of how to apply the ODTR SuperLearner to real
data, we estimated the ODTR using the “Interventions” study to determine which
types of justice-involved adults with mental illness should be assigned CBT
versus TAU, to yield the highest probability of no re-arrest. Preliminary
results show that the ODTR SuperLearner placed all weight on a simple
parametric model with only substance use; thus, the algorithm suggests that,
in this sample, participants with lower levels of substance use may benefit
more from CBT. We note that this is an example of a case in which the ODTR
SuperLearner generated a ODTR estimate that was fully interpretable – although
we used a continuous metalearner and thus the SuperLearner could have allowed
for combinations of algorithms, the SuperLearner happened to only choose one
algorithm: a GLM with substance use as a regressor. To guarantee
interpretability in the SuperLearner (for example, if working with
practitioners who may want a treatment decision rule that could be easily
written down [3, 62]), one could implement the ODTR SuperLearner with a
discrete metalearner and a simple parametric library only.

Importantly, in a companion paper, we _evaluate_ this ODTR – that is, we ask
the causal question: what would have been the probability of no re-arrest had
participants been assigned CBT according the ODTR SuperLearner (i.e., using
only substance use)? Further, is assigning CBT according to the ODTR
SuperLearner significantly better than assigning CBT to everyone or no one? In
this way, we can determine if it is of clinical significance to assign CBT
according to this rule – namely, if assigning CBT using only substance use
scores results in a statistically significant reduction of recidivism, and if
so, how much better one does with this ODTR compared to a non-individualized
rule (such as giving CBT to all). Under the appropriate causal assumptions,
one could use any of the methods for estimating treatment specific means to
interpret this estimate as the expected outcome under the true optimal rule or
the estimated optimal rule.

Future work could extend these simulations to the multiple treatment (i.e.,
more than 2 treatment levels) [35] and multiple timepoint setting (i.e.,
estimating a sequential ODTR from, for example, a SMART design [7, 25, 36]
instead of an RCT design). We also wish to apply the ODTR SuperLearner on the
full “Interventions” dataset (n = 720), once data collection is finished.

This work contributes to understanding the toolbox of methods for analyzing
the heterogeneity in how patients respond to treatment. By learning which
patients respond best to what treatment in a flexible manner, we can improve
patient outcomes – moving us closer to the goals of precision health.

| TAU ($A=0$) | CBT ($A=1$) | $p$
---|---|---|---
$n$ | 211 | 230 |
No re-arrest ($Y=1$) (%) | 128 (60.7) | 143 (62.2) | 0.820
Site = San Francisco (%) | 87 (41.2) | 104 (45.2) | 0.455
Gender = Female (%) | 38 (18.0) | 37 (16.1) | 0.682
Ethnicity = Hispanic (%) | 50 (23.7) | 42 (18.3) | 0.198
Age (mean (SD)) | 38.08 (11.05) | 37.01 (11.22) | 0.317
CSI (mean (SD)) | 32.35 (11.13) | 33.46 (11.27) | 0.300
LSI (mean (SD)) | 5.59 (1.33) | 5.50 (1.48) | 0.472
SES (mean (SD)) | 3.81 (1.89) | 3.81 (2.12) | 0.995
Prior adult convictions (%) |  |  | 0.156
Zero to two times | 74 (35.1) | 93 (40.4) |
Three or more times | 134 (63.5) | 129 (56.1) |
Missing | 3 (1.4) | 8 (3.5) |
Most serious offense (mean (SD)) | 5.29 (2.54) | 5.09 (2.52) | 0.415
Motivation (mean (SD)) | 3.22 (1.36) | 3.27 (1.37) | 0.720
Substance use (%) |  |  | 0.184
0 | 53 (25.1) | 76 (33.0) |
1 | 47 (22.3) | 55 (23.9) |
2 | 109 (51.7) | 98 (42.6) |
Missing | 2 (0.9) | 1 (0.4) |

Table 2: Distribution of “Interventions” data by treatment assignment. Figure
1: Performance of $E_{n}[Q_{0}(Y|A=d_{n}^{*},W)]$ (an approximation of the
true mean outcome under the estimated ODTR) for DGP 1 (top two) and DGP 2
(bottom two). The horizontal black line depicts $E_{P_{U,X}}[Y_{d_{0}^{*}}]$;
the horizontal red line depicts $E_{P_{U,X}}[Y_{1}]$; the horizontal blue line
depicts $E_{P_{U,X}}[Y_{0}]$. We compare the SuperLearner algorithm to an
incorrectly specified GLM (in gray, with N/A as the metalearner and a diamond
with no fill). We also compare (1) having a SuperLearner library with (a) only
algorithms that estimate the blip (i.e., “Blip only” libraries) that only have
parametric blip models (blue) or only have machine-learning blip models (red)
or both (purple) versus (b) an expanded or “Full” library with blip function
regressions estimated via machine learning only (orange-yellow) or machine
learning and parametric models (green), with methods that directly estimate
the ODTR and static rules, (2) having a metalearner (depicted on the x-axis)
either that chooses one algorithm (i.e., the “discrete” SuperLearner) or
combines blip predictions/treatment predictions (i.e., the “continuous”
SuperLearner), and (3) using the MSE risk function ($R_{MSE}$ as a square)
versus the mean outcome under the candidate rule risk function ($R_{E[Y_{d}]}$
as a triangle). Figure 2: Subgroup plots for each of the covariates for the
“Interventions” data. The x-axis for each of the plots is the different levels
of the covariates; the y-axis is the difference in proportion of people who
were not re-arrested between those who received CBT versus TAU, in that
covariate subgroup. Figure 3: Distribution of predicted blip estimates from
the ODTR SuperLearner. The frequencies are divided into three groups because
the ODTR SuperLearner allocated all coefficient weights to a GLM using
substance use, a variable with only 3 treatment levels. One can interpret the
ODTR SuperLearner for this sample as follows: CBT may reduce the probability
of re-arrest among justice-involved adults with low levels of substance use.
Estimation and inference of the value of the ODTR SuperLearner compared to,
for example, treating everyone or no one, informs us if there is, in fact, a
differential effect by substance use, and thus a benefit to assigning CBT in
this individualized way.

## References

  * [1] Issa J Dahabreh, Rodney Hayward, and David M Kent.  Using group data to treat individuals: understanding heterogeneous treatment effects in the age of precision medicine and patient-centred evidence.  International journal of epidemiology, 45(6):2184–2193, 2016.
  * [2] David M Kent, Ewout Steyerberg, and David van Klaveren.  Personalized evidence based medicine: predictive approaches to heterogeneous treatment effects.  Bmj, 363:k4245, 2018.
  * [3] Michael R Kosorok and Eric B Laber.  Precision medicine.  Annual review of statistics and its application, 6:263–286, 2019\.
  * [4] Jennifer L Skeem, Sarah Manchak, and Jillian K Peterson.  Correctional policy for offenders with mental illness: Creating a new paradigm for recidivism reduction.  Law and human behavior, 35(2):110–126, 2011.
  * [5] Jennifer L Skeem, Eliza Winter, Patrick J Kennealy, Jennifer Eno Louden, and Joseph R Tatar II.  Offenders with mental illness have criminogenic needs, too: Toward recidivism reduction.  Law and human behavior, 38(3):212, 2014.
  * [6] Mark W Lipsey, Nana A Landenberger, and Sandra J Wilson.  Effects of cognitive-behavioral programs for criminal offenders.  Campbell systematic reviews, 3(1):1–27, 2007.
  * [7] Huitian Lei, Inbal Nahum-Shani, K Lynch, David Oslin, and Susan A Murphy.  A” smart” design for building individualized treatment sequences.  Annual review of clinical psychology, 8:21–48, 2012.
  * [8] Feifang Hu and William F Rosenberger.  The theory of response-adaptive randomization in clinical trials, volume 525.  John Wiley & Sons, 2006.
  * [9] Ilya Lipkovich, Alex Dmitrienko, and Ralph B D’Agostino Sr.  Tutorial in biostatistics: data-driven subgroup identification and analysis in clinical trials.  Statistics in medicine, 36(1):136–196, 2017.
  * [10] Oliver Bembom and Mark J van der Laan.  A practical illustration of the importance of realistic individualized treatment rules in causal inference.  Electronic journal of statistics, 1:574, 2007.
  * [11] Mark J van der Laan and Maya L Petersen.  Causal effect models for realistic individualized treatment and intention to treat rules.  The international journal of biostatistics, 3(1), 2007.
  * [12] James Robins.  A new approach to causal inference in mortality studies with a sustained exposure period—application to control of the healthy worker survivor effect.  Mathematical modelling, 7(9-12):1393–1512, 1986.
  * [13] Bibhas Chakraborty and Erica EM Moodie.  Statistical methods for dynamic treatment regimes.  Springer, 2013.
  * [14] Susan A Murphy.  Optimal dynamic treatment regimes.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 65(2):331–355, 2003.
  * [15] James M Robins.  Optimal structural nested models for optimal sequential decisions.  In Proceedings of the second seattle Symposium in Biostatistics, pages 189–326. Springer, 2004.
  * [16] Alexander R Luedtke and Mark J van der Laan.  Optimal individualized treatments in resource-limited settings.  The international journal of biostatistics, 12(1):283–303, 2016\.
  * [17] Eric B Laber, Kristin A Linn, and Leonard A Stefanski.  Interactive model building for q-learning.  Biometrika, 101(4):831–847, 2014.
  * [18] Phillip J Schulte, Anastasios A Tsiatis, Eric B Laber, and Marie Davidian.  Q-and a-learning methods for estimating optimal dynamic treatment regimes.  Statistical science: a review journal of the Institute of Mathematical Statistics, 29(4):640, 2014.
  * [19] Erica EM Moodie, Bibhas Chakraborty, and Michael S Kramer.  Q-learning for estimating optimal dynamic treatment rules from observational data.  Canadian Journal of Statistics, 40(4):629–645, 2012.
  * [20] Min Qian and Susan A Murphy.  Performance guarantees for individualized treatment rules.  Annals of statistics, 39(2):1180, 2011.
  * [21] Erica EM Moodie, Thomas S Richardson, and David A Stephens.  Demystifying optimal dynamic treatment regimes.  Biometrics, 63(2):447–455, 2007.
  * [22] Yingqi Zhao, Donglin Zeng, A John Rush, and Michael R Kosorok.  Estimating individualized treatment rules using outcome weighted learning.  Journal of the American Statistical Association, 107(499):1106–1118, 2012.
  * [23] Baqun Zhang, Anastasios A Tsiatis, Marie Davidian, Min Zhang, and Eric Laber.  Estimating optimal treatment regimes from a classification perspective.  Stat, 1(1):103–114, 2012.
  * [24] Ying-Qi Zhao, Donglin Zeng, Eric B Laber, and Michael R Kosorok.  New statistical learning methods for estimating optimal dynamic treatment regimes.  Journal of the American Statistical Association, 110(510):583–598, 2015.
  * [25] Michael R Kosorok and Erica EM Moodie.  Adaptive Treatment Strategies in Practice: Planning Trials and Analyzing Data for Personalized Medicine, volume 21.  SIAM, 2015.
  * [26] Ying-Qi Zhao and Eric B Laber.  Estimation of optimal dynamic treatment regimes.  Clinical Trials, 11(4):400–407, 2014.
  * [27] Mark J van der Laan, Eric C Polley, and Alan E Hubbard.  Super learner.  Statistical applications in genetics and molecular biology, 6(1), 2007.
  * [28] Leo Breiman.  Stacked regressions.  Machine learning, 24(1):49–64, 1996.
  * [29] Erin LeDell, Mark J van der Laan, and Maya Petersen.  Auc-maximizing ensembles through metalearning.  The international journal of biostatistics, 12(1):203–218, 2016\.
  * [30] Eric C Polley and Mark J van der Laan.  Super learner in prediction.  2010\.
  * [31] Romain Pirracchio, Maya L Petersen, and Mark van der Laan.  Improving propensity score estimators’ robustness to model misspecification using super learner.  American journal of epidemiology, 181(2):108–119, 2014.
  * [32] Romain Pirracchio, Maya L Petersen, Marco Carone, Matthieu Resche Rigon, Sylvie Chevret, and Mark J van der Laan.  Mortality prediction in intensive care units with the super icu learner algorithm (sicula): a population-based study.  The Lancet Respiratory Medicine, 3(1):42–52, 2015.
  * [33] Maya L Petersen, Erin LeDell, Joshua Schwab, Varada Sarovar, Robert Gross, Nancy Reynolds, Jessica E Haberer, Kathy Goggin, Carol Golin, Julia Arnsten, et al.  Super learner analysis of electronic adherence data improves viral prediction and may provide strategies for selective hiv rna monitoring.  Journal of acquired immune deficiency syndromes (1999), 69(1):109, 2015.
  * [34] Alexander R Luedtke and Mark J van der Laan.  Super-learning of an optimal dynamic treatment rule.  The international journal of biostatistics, 12(1):305–332, 2016\.
  * [35] Jeremy Robert Coyle.  Computational Considerations for Targeted Learning.  PhD thesis, UC Berkeley, 2017.
  * [36] Daniel Almirall, Inbal Nahum-Shani, Nancy E Sherwood, and Susan A Murphy.  Introduction to smart designs for the development of adaptive interventions: with application to weight loss research.  Translational behavioral medicine, 4(3):260–274, 2014.
  * [37] Maya L Petersen and Mark J van der Laan.  Causal models and learning from data: integrating causal modeling and statistical estimation.  Epidemiology (Cambridge, Mass.), 25(3):418, 2014.
  * [38] Judea Pearl.  Causality: models, reasoning and inference, volume 29.  Springer, 2000.
  * [39] Alexander R Luedtke and Mark J van der Laan.  Statistical inference for the mean outcome under a possibly non-unique optimal treatment strategy.  Annals of statistics, 44(2):713, 2016.
  * [40] James M Robins and Miguel A Hernán.  Estimation of the causal effects of time-varying exposures.  Longitudinal data analysis, 553:599, 2009.
  * [41] Maya L Petersen, Kristin E Porter, Susan Gruber, Yue Wang, and Mark J van der Laan.  Diagnosing and responding to violations in the positivity assumption.  Statistical methods in medical research, 21(1):31–54, 2012.
  * [42] Inbal Nahum-Shani, Min Qian, Daniel Almirall, William E Pelham, Beth Gnagy, Gregory A Fabiano, James G Waxmonsky, Jihnhee Yu, and Susan A Murphy.  Q-learning: A data analysis method for constructing adaptive interventions.  Psychological methods, 17(4):478, 2012.
  * [43] Daniel Rubin and Mark J van der Laan.  A doubly robust censoring unbiased transformation.  The international journal of biostatistics, 3(1), 2007.
  * [44] Mark J van der Laan and Alexander R Luedtke.  Targeted learning of the mean outcome under an optimal dynamic treatment rule.  Journal of causal inference, 3(1):61–95, 2015.
  * [45] Daniel B Rubin and Mark J van der Laan.  Statistical issues and limitations in personalized medicine research with clinical trials.  The international journal of biostatistics, 8(1), 2012.
  * [46] Mark J van der Laan and Sherri Rose.  Targeted learning: causal inference for observational and experimental data.  Springer Science & Business Media, 2011.
  * [47] Brian D Ripley and NL Hjort.  Pattern recognition and neural networks.  Cambridge university press, 1996.
  * [48] Jerome H Friedman et al.  Multivariate adaptive regression splines.  The annals of statistics, 19(1):1–67, 1991.
  * [49] Leo Breiman.  Classification and regression trees.  Routledge, 2017.
  * [50] Mark J van der Laan and Susan Gruber.  Targeted minimum loss based estimation of an intervention specific mean outcome.  2011\.
  * [51] Ravi Varadhan, Jodi B Segal, Cynthia M Boyd, Albert W Wu, and Carlos O Weiss.  A framework for the analysis of heterogeneity of treatment effect in patient-centered outcomes research.  Journal of clinical epidemiology, 66(8):818–825, 2013.
  * [52] David van Klaveren, Yvonne Vergouwe, Vasim Farooq, Patrick W Serruys, and Ewout W Steyerberg.  Estimates of absolute treatment benefit for individual patients required careful modeling of statistical interactions.  Journal of clinical epidemiology, 68(11):1366–1374, 2015.
  * [53] Salim Yusuf, Janet Wittes, Jeffrey Probstfield, and Herman A Tyroler.  Analysis and interpretation of treatment effects in subgroups of patients in randomized clinical trials.  Jama, 266(1):93–98, 1991.
  * [54] R Core Team.  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria, 2018.
  * [55] Chih-Chung Chang and Chih-Jen Lin.  Libsvm: A library for support vector machines.  ACM transactions on intelligent systems and technology (TIST), 2(3):27, 2011.
  * [56] Eric Polley, Erin LeDell, Chris Kennedy, and Mark van der Laan.  SuperLearner: Super Learner Prediction, 2018.  R package version 2.0-24.
  * [57] Xin Zhou, Nicole Mayer-Hamblett, Umer Khan, and Michael R Kosorok.  Residual weighted learning for estimating individualized treatment rules.  Journal of the American Statistical Association, 112(517):169–187, 2017.
  * [58] Ying-Qi Zhao, Eric B Laber, Yang Ning, Sumona Saha, and Bruce E Sands.  Efficient augmentation and relaxation learning for individualized treatment rules using observational data.  Journal of Machine Learning Research, 20(48):1–23, 2019.
  * [59] Baqun Zhang, Anastasios A Tsiatis, Marie Davidian, Min Zhang, and Eric Laber.  Estimating optimal treatment regimes from a classification perspective.  Stat, 1(1):103–114, 2012.
  * [60] S. T. Holloway, E. B. Laber, K. A. Linn, B. Zhang, M. Davidian, and A. A. Tsiatis.  DynTxRegime: Methods for Estimating Optimal Dynamic Treatment Regimes, 2019.  R package version 4.1.
  * [61] Kelly L Moore and Mark J van der Laan.  Covariate adjustment in randomized trials with binary outcomes: targeted maximum likelihood estimation.  Statistics in medicine, 28(1):39–64, 2009.
  * [62] Zachary D Cohen and Robert J DeRubeis.  Treatment selection in depression.  Annual Review of Clinical Psychology, 14, 2018.