Targeted Learning Ensembles for Optimal Individualized Treatment Rules with Time-to-Event Outcomes

TL;DR

This paper develops targeted ensemble methods for estimating optimal individualized treatment rules to maximize time-to-event outcomes, using semiparametric efficiency theory and super learning in clinical trial data.

Contribution

It introduces two novel estimators for optimal treatment rules based on different loss functions, with theoretical guarantees and application to cancer treatment data.

Findings

Proposed estimators are doubly robust and achieve oracle inequalities.

Ensemble methods improve the estimation of individualized treatment rules.

Application demonstrates effectiveness in a breast cancer clinical trial.

Abstract

We consider estimation of an optimal individualized treatment rule from observational and randomized studies when a high-dimensional vector of baseline variables is available. Our optimality criterion is with respect to delaying expected time to occurrence of an event of interest (e.g., death or relapse of cancer). We leverage semiparametric efficiency theory to construct estimators with desirable properties such as double robustness. We propose two estimators of the optimal rule, which arise from considering two loss functions aimed at (i) directly estimating the conditional treatment effect (also know as the blip function), and (ii) recasting the problem as a weighted classification problem that uses the 0-1 loss function. Our estimated rules are super learning ensembles that minimize the cross-validated risk of a linear combination in a user-supplied library of candidate estimators. We prove oracle inequalities bounding the finite sample excess risk of the estimator. The bounds depend on the excess risk of the oracle selector and a doubly robust term related to estimation of the nuisance parameters. We discuss some important implications of these oracle inequalities such as the convergence rates of the value of our estimator to that of the oracle selector. We illustrate our methods in the analysis of a phase III randomized study testing the efficacy of a new therapy for the treatment of breast cancer.

Tables

Table 1. Table 1: Coefficients of the super learner ensemble for estimation of g R subscript 𝑔 𝑅 g_{R} and h ℎ h . RF is random forests, XGB is extreme gradient boosting, MLP is multilayer perceptron, GLM is logistic regression, MARS is multivariate adaptive splines, and LASSO is L 1 subscript 𝐿 1 L_{1} regularized logistic regression.

		RF	XGB	MLP	GLM	MARS	LASSO
${\hat{g}}_R$	$A=1$	0.100	0.056	0.000	0.000	0.312	0.532
	$A=0$	0.000	0.301	0.000	0.000	0.294	0.405
$\hat{h}$	$A=1$	0.023	0.050	0.000	0.000	0.237	0.691
	$A=0$	0.154	0.086	0.098	0.060	0.123	0.480

Table 2. Table 2: Coefficients of each candidate in each ensemble (standardized to sum one). The rows represent the candidates, the columns the ensemble.

	SL-Reg	SL-Class-0-1	SL-Class-Hinge	SL-Class-Log
D-Class-RF	0.000	0.005	0.000	0.000
D-Class-XGB	0.792	0.001	0.945	0.869
D-Class-GLM	0.000	0.031	0.000	0.000
D-Reg	0.007	0.017	0.001	0.006
B-Reg	0.201	0.947	0.054	0.125

Table 3. Table 3: P-values of pair-wise comparisons of the value of each rule estimated in the validation data.

SL-Class-Log	0.026
SL-Class-Hinge	0.026	0.159
SL-Reg	0.025	0.160	0.443
D-Class-RF	0.001	0.001	0.001	0.001
D-Class-XGB	0.001	0.009	0.007	0.007	0.003
D-Class-GLM	0.001	0.001	0.001	0.001	0.011	0.001
D-Reg	0.001	0.004	0.002	0.002	0.158	0.021	0.056
B-Reg	0.120	0.025	0.023	0.023	0.001	0.001	0.001	0.001
Always-Treat	0.001	0.001	0.001	0.001		0.003	0.011	0.158	0.001
Never-Treat	0.001	0.001	0.001	0.001	0.001	0.001	0.036	0.002	0.001	0.001
	SL-Class-0-1	SL-Class-Log	SL-Class-Hinge	SL-Reg	D-Class-RF	D-Class-XGB	D-Class-GLM	D-Reg	B-Reg	Always-Treat

Equations

$$O=(W,A,R_0,L_1,R_1,L_2\dots ,R_{K-1},L_K),$$

$$h(m,a,w)=P_0(L_m=1\mid I_m=1,A=a,W=w),$$

$$g_R(m,a,w)=P_0(R_m=1\mid J_m=1,A=a,W=w).$$

$${\theta }_{\text{mbox}c}(z)=E\{\min (T_1,\tau )-\min (T_0,\tau )\mid Z=z\}.$$

$${\theta }_{\text{mbox}c}(z)=\sum _{t=1}^{\tau -1}\{S_{\text{mbox}c}(t,1,z)-S_{\text{mbox}c}(t,0,z)\},$$

$$E\{\min (T_d,\tau )\}=E\{d(Z)[\min (T_1,\tau )-\min (T_0,\tau )]\}+E\{\min (T_0,\tau )\}.$$

$$V_{\text{mbox}c}(d)=E\{d(Z)[\min (T_1,\tau )-\min (T_0,\tau )]\}=E\{d(Z){\theta }_{\text{mbox}c}(Z)\}.$$

$$d_{\text{mbox}c}(z)=\underset{d\in \mathcal{D}}{\mathrm{argmax}}{V}_{\text{mbox}c}(d)=\text{mathds}1\{{\theta }_{\text{mbox}c}(z)>0\},$$

$$S(t,a,w)=P(T>t\mid A=a,W=w),G(t,a,w)=P(C\ge t\mid A=a,W=w).$$

$$S(t,a,w)$$

$$\theta (z)=\sum _{t=1}^{\tau -1}E\left\{\prod _{m=1}^t\{1-h(m,1,W)\}-\prod _{m=1}^t\{1-h(m,0,W)\}\mid Z=z\right\}.$$

$$D_{\eta}(O)=\sum_{m=1}^{\tau-1}\big{[}I_{m}Z(m,A,W)\left\{L_{m}-h(m,A,W)\right\}+S(m,1,W)-S(m,0,W)\big{]},$$

$$Z_a(m,A,W)=-\sum _{t=m}^{\tau -1}\frac{\text{mathds}1\{A=a\}}{g_A(a,W)G(m,a,W)}\frac{S(t,a,W)}{S(m,a,W)}.$$

$${\hat{R}}_{\hat{\eta }}({\hat{\theta }}_j)=\frac{1}{K}\sum _{i=1}^n\frac{1}{|{\mathcal{V}}_{k(i)}|}{L}_{{\hat{\eta }}_{k(i)}}(O_i,{\hat{\theta }}_{j,k(i)}).$$

$${\hat{\theta }}_{\alpha }(z)=\sum _{j=1}^J{\alpha }_j{\hat{\theta }}_j(z),{\alpha }_j\ge 0,\sum _{j=1}^J{\alpha }_j=1.$$

$$\hat{\alpha }=\underset{\alpha }{\mathrm{argmin}}\sum _{i=1}^n\frac{1}{|{\mathcal{V}}_{k(i)}|}{\left\{D_{{\hat{\eta }}_{k(i)}}(O_i)-\sum _{j=1}^J{\alpha }_j{\hat{\theta }}_{j,k(i)}(V_i)\right\}}^2\text{ subject to }{\alpha }_j\ge 0,\sum _{j=1}^J{\alpha }_j=1.$$

$${\tilde{R}}_{{\eta }_1}({\hat{\theta }}_{\alpha })=\frac{1}{K}\sum _{k=1}^K\int {\left\{D_{{\eta }_1}(o)-{\hat{\theta }}_{\alpha ,k}(z)\right\}}^{2}d{P}_0(o).$$

$$\tilde{\alpha }=\underset{\alpha }{\mathrm{argmin}}{\tilde{R}}_{{\eta }_1}({\hat{\theta }}_{\alpha })\text{ subject to }{\alpha }_j\ge 0,\sum _{j=1}^J{\alpha }_j=1,$$

$$E^2({\hat{\theta }}_{\hat{\alpha }})=\mathbb{E}\{{\tilde{R}}_{{\eta }_1}({\hat{\theta }}_{\hat{\alpha }})-{\tilde{R}}_{{\eta }_1}({\theta }_0)\}=\mathbb{E}{P}_0({\hat{\theta }}_{\alpha }-{\theta }_0{)}^2,$$

$$B_1(\hat{\eta },{\eta }_0)$$

$$B_2(\hat{\eta },{\eta }_0)$$

$$\begin{gathered} E({\hat{\theta }}_{\mathrm{sl}})\le (1+2\delta {)}^{1/2}E({\hat{\theta }}_{\mathrm{or}})+C_1{\left\{(1+\log n)/n\right\}}^{1/2}+ \\ {C}_2{B}_1(\hat{\eta },{\eta }_0)+C_3{\left(\log n/n\right)}^{1/4}{\left\{B_2(\hat{\eta },{\eta }_0)\right\}}^{1/2} \end{gathered}$$

$$B_1(\hat{\eta },{\eta }_0)=O\left((\log n/n{)}^{1/2}\right),B_2(\hat{\eta },{\eta }_0)=O\left((\log n/n{)}^{1/2}\right).$$

$$\mathbb{E}\{V_0(\tilde{d})-V_0(\hat{d})\}\lesssim \{E^2({\hat{\theta }}_{\mathrm{sl}})-E^2({\hat{\theta }}_{\mathrm{or}}){\}}^{(1+\lambda )/(2+\lambda )},$$

$$\mathbb{E}\{V_0(\tilde{d})-V_0(\hat{d})\}=O\left((\log n/n{)}^{(1+\lambda )/(2+\lambda )}\right).$$

$$f_0\in \underset{f}{\mathrm{argmax}}\int d_f(z)D_{\eta }(o)d{P}_0(o).$$

$$L_{\eta }(o;f)=|D_{\eta }(o)|\text{mathds}1\left[\text{mathds}1\{D_{\eta }(o)>0\}\ne d_f(z)\right].$$

$${\hat{f}}_{\alpha }(z)=\sum _{j=1}^J{\alpha }_j{\hat{f}}_j(z),{\alpha }_j\ge 0.$$

$$\hat{\alpha }\in \underset{\alpha }{\mathrm{argmin}}\sum _{i=1}^n\frac{1}{|{\mathcal{V}}_{k(i)}|}{L}_{{\hat{\eta }}_k(i)}\left(O_i,{\hat{f}}_{\alpha ,k(i)}\right)\text{ subject to }{\alpha }_j\ge 0,$$

$${\tilde{R}}_{{\eta }_1}(\hat{f})=\frac{1}{K}\sum _{k=1}^K\int L_{{\eta }_1}(o,{\hat{f}}_k)d{P}_0(o).$$

$$E(\hat{f})=\mathbb{E}\{{\tilde{R}}_{{\eta }_1}(\hat{f})-{\tilde{R}}_{{\eta }_1}(f_0)\}=V_0(f_0)-\mathbb{E}{V}_0(\hat{f}).$$

Full text

1 Introduction

Individualized treatment rules play a fundamental role in the precision medicine model for healthcare, whereby medical decisions are targeted to the individual based on their expected clinical response, instead of the traditional one-size-fits-all approach. Mathematically, a treatment rule is a function that maps an individual’s pre-treatment covariates into an optimal treatment choice. In this paper, we are concerned with learning the optimal rules from data collected as part of an observational or randomized study, where optimality is defined as the maximum delay in the expected time of occurrence of an undesirable event (e.g., death or relapse).

Recent advances in biomedical imaging and gene expression technology produce large amounts of data that can be used to tailor treatment to very specific patient characteristics. Methods to estimate the optimal rule when it is defined with respect to a single time-point outcome include the work of Qian and Murphy (2011); Zhao et al. (2012); Song et al. (2015); Rubin et al. (2012); McKeague and Qian (2014), among others. Methods to solve the problem using survival outcomes subject to informative censoring have been proposed by Zhao et al. (2011); Goldberg and Kosorok (2012). The latter methods use Q-learning, relying on sequential support vector regressions, to estimate the optimal sequential treatment rule that optimizes a survival outcome under right-censoring. Geng et al. (2015) also tackle estimation of the optimal rule in a survival setting using $\ell_{1}$ regularization for the outcome regression under the strong assumption that censoring is independent of covariates and the outcome, but their decision functions are restricted to linear functions. Zhao et al. (2015) generalize the weighted classification approach of Zhao et al. (2012) to allow for informative censoring and doubly robust loss functions, but their decision functions are restricted to support vector machines. Bai et al. (2016) present methods for estimating optimal rules with a survival outcome subject to informative censoring. They consider two strategies based on estimation of the blip function and based on a classification perspective. Their methods are restricted to decision functions that can be indexed by a Euclidean vector and parametric nuisance estimators, and are therefore of limited applicability to high-dimensional data. All the above methods are potential candidates in the library of estimators that constitute our ensembles.

In this article, we propose two methods to construct an ensemble of decision functions for the optimal rule. Our ensembles are linear combinations of estimators in a user-supplied library, where the coefficients in the linear combination are chosen to minimize the cross-validated risk. We propose to use a doubly robust loss function with roots in efficient estimation theory for marginal causal effects (Moore and van der Laan, 2009; Díaz et al., 2015). In our context, double robustness means that the estimated rules will have certain optimality properties under consistent estimation of at least one of two nuisance parameters: (a) the hazard of the outcome at each time point conditional on covariates and treatment, and (b) the hazard of censoring and the treatment mechanism.

The library of candidate estimators may contain any of the algorithms discussed in the previous paragraphs. In light of the no free lunch theorems of Wolpert (2002) for supervised learning, for any given dataset, our ensembles are expected to have better or equal generalization error than any of the individual candidates in the library. We provide a formal proof of this claim in the form of an oracle inequality, which bounds the excess risk of our estimator in terms of the excess risk of the oracle estimator, defined as the combination of estimators that would be chosen in a hypothetical world in which an infinite validation sample is available and at least one of the nuisance parameters is known. Our methods are developed under the assumption that censoring is at random (Rubin, 1987), which means that censoring is random within strata of treatment and baseline variables. We also assume that treatment is randomized within strata of the covariates, either by nature or by experimentation.

The finite sample bounds we present are inspired by developments in the targeted learning literature, which establish the optimality of cross-validation in estimator selection for high-dimensional parameters (van der Laan and Dudoit, 2003). Related to our work, Luedtke and van der Laan (2016) consider super learning ensembles for estimation of optimal DTRs in two time points. They present oracle inequalities for super learning of the optimal rule using a loss function indexed by the treatment mechanism, which is assumed known. We generalize their results in the following ways: (i) we provide oracle inequalities under a doubly robust loss function indexed by two nuisance parameters, when neither of the nuisance parameters is known, (ii) we show that the oracle inequalities inherit the double robustness property of the loss function, and (iii) we present comparable oracle inequalities for the 0-1 loss function. In addition, we discuss how these oracle inequalities are related to the convergence of the value of the rule under a margin assumption describing the behavior of the blip function in the boundary of the decision threshold.

2 Data and Notation

Assume individuals are monitored at $K$ time points $t=\{1,\dots,K\}$. Let $T$ denote a time-to-event outcome taking values in $\{1,\dots,K\}\cup\{\infty\}$, where $T=\infty$ represents no event occurring in the follow-up period. Let $C\in\{0,\dots,K\}$ denote the censoring time defined as the time at which the individual is last observed in the study, and let $C=K$, represent administrative censoring. Let $A\in\{0,1\}$ denote study arm assignment, and let $W$ denote a vector of baseline variables, which may include gene expression as well as demographic, comorbidity, and other clinical data. Denote $\mathds{1}(\cdot)$ the indicator variable taking value $1$ if the argument is true and [math] otherwise. The observed data vector for each participant is $O=(W,A,\Delta,\tilde{T})$, where $\tilde{T}=\min(C,T)$, and $\Delta=\mathds{1}\{T\leq C\}$ is the indicator that the participant’s event time is observed (uncensored). For a random variable $X$, we let $X$ take values on a set $\mathbf{O}$.

We assume the observed data vector for each participant $i$, denoted $O_{i}=(W_{i},A_{i},\Delta_{i},\tilde{T}_{i})$, is an independent, identically distributed draw from the unknown joint distribution $P_{0}$ on $(W,A,\Delta,\tilde{T})$. The empirical distribution of $O_{1},\ldots,O_{n}$ is denoted with $P_{n}$. We assume $P_{0}\in\mathcal{M}$, where $\mathcal{M}$ is the nonparametric model defined as all continuous densities on $O$ with respect to a dominating measure $\nu$. We use $P$ to denote a generic distribution $P\in\mathcal{M}$, and $E_{0}(\cdot)$ to denote expectation with respect to $P_{0}$, and $\mathbb{E}(\cdot)$ is used to denote expectation over draws of $O_{1},\ldots,O_{n}$. For a function $f(o)$, we denote $Pf=\int f(o)dP(o)$, and $||f||^{2}=P_{0}f^{2}$. We use $a\lesssim b$ to denote that $a$ is smaller or equal than $b$ up to a universal constant.

We can equivalently encode a single participant’s data vector $O$ using the following longitudinal data structure:

[TABLE]

where $R_{t}=\mathds{1}\{\tilde{T}=t,\Delta=0\}$ and $L_{t}=\mathds{1}\{\tilde{T}=t,\Delta=1\}$, for $t\in\{0,\ldots,K\}$. For a random variable $X$, we denote its history through time $t$ as $\bar{X}_{t}=(X_{0},\ldots,X_{t})$. For a given scalar $x$, the expression $\bar{X}_{t}=x$ denotes element-wise equality.

Define the following indicator variables for each $t\geq 1$: $I_{t}=\mathds{1}\{\bar{R}_{t-1}=0,\bar{L}_{t-1}=0\}$, $J_{t}=\mathds{1}\{\bar{R}_{t-1}=0,\bar{L}_{t}=0\}.$ The variable $I_{t}$ is the indicator based on the data through time $t-1$ that a participant is at risk of the event being observed at time $t$. Analogously, $J_{t}$ is the indicator based on the outcome data through time $t$ and censoring data before time $t$ that a participant is at risk of censoring at time $t$. We define $J_{0}=1$.

Define the discrete hazard function for survival at time $m\in\{1,\dots,K\}$:

[TABLE]

among the population at risk at time $m$ within strata of study arm and baseline variables. Similarly, for the censoring variable $C$, define the censoring hazard at time $m\in\{0,\dots,K\}$:

[TABLE]

We use the notation $g_{A}(a,w)=P_{0}(A=a\mid W=w)$, $g=(g_{A},g_{R})$, and $\eta=(h,g_{A},g_{R})$. Let $p_{W}$ denote the marginal distribution of the baseline variables $W$. We add the subscript [math] to $p_{W},g,h$ to denote the corresponding quantities under $P_{0}$.

3 Treatment Effect, Identification, and Optimal Individualized

Treatment Rules

3.1 Potential Outcomes and Causal Parameter

Define the potential outcomes $T_{a}:a\in\{0,1\}$ as the event times that would have been observed had study arm assignment $A=a$ and censoring time $C=K$ been externally set with probability one. For a restriction time $\tau\in\{1,\ldots,K\}$ of interest, we define the restricted survival time under treatment arm $A=a$ as $\min(T_{a},\tau)$. For a transformation $Z$ of $W$, the treatment effect within strata of the covariates $Z$ may be defined in terms of the so-called full-data blip function (see e.g., Robins, 1997) of the restricted mean survival time:

[TABLE]

The transformation $Z$ may represent a subset of covariates (e.g., gene expression), or the whole vector $W$. We define the marginal treatment effect as $\theta_{\mbox{\footnotesize c,m}}=E\{\min(T_{1},\tau)-\min(T_{0},\tau)\}$.

The subscript $c$ denotes a causal parameter, that is, a parameter of the distribution of the potential outcomes $T_{1}$ and $T_{0}$. It can be shown (see Díaz et al., 2015) that $E\{\min(T_{a},\tau)\mid Z=z\}=\sum_{t=0}^{\tau-1}S_{\mbox{\footnotesize c}}(t,a,z),$ where $S_{\mbox{\footnotesize c}}(t,a,z)=P(T_{a}>t\mid Z=z)$ is the survival probability corresponding to the potential outcome under assignment to arm $A=a$ within strata $Z=z$. As a result, $\theta_{\mbox{\footnotesize c}}(z)$ may be expressed as

[TABLE]

since $S_{\mbox{\footnotesize c}}(0,a,z)=1$ for $a\in\{0,1\}$ and for all $z$.

An individualized treatment rule $d$ is a function that maps the covariate values $z$ of a given participant to a personalized treatment decision in $\{0,1\}$. The potential time to event under a rule $d$ is defined as $T_{d}=d(z)T_{1}+\{1-d(z)\}T_{0}$. Accordingly, the restricted mean survival time under a treatment rule that assigns treatment according to $d(z)$ is equal to

[TABLE]

Because the last term does not depend on $d(z)$, we define the value of the rule $d$ as

[TABLE]

The above equation provides the basis for the definition of an optimal rule as

[TABLE]

where $\mathcal{D}=\{d:\mathbf{Z}\to\{0,1\}\}$ is the space of functions that map the range of $Z$ into a treatment decision in $\{0,1\}$. We define optimality of an rule with respect to the restricted mean survival time, though other effect measures could also be used.

3.2 Identification of Parameters in Terms of Observed Data

Generating Distribution $P_{0}$

In this section we show how the blip function $\theta_{\mbox{\footnotesize c}}(z)$, the value function $V_{\mbox{\footnotesize c}}(d)$, and the optimal rule $d_{\mbox{c}}(z)$, which are defined above in terms of the distribution of potential outcomes, can be equivalently expressed as functions $\theta_{0}(z)$, $V_{0}(d)$, and $d_{0}(z)$ of the observed data distribution $P_{0}(W,A,\Delta,\tilde{T})$, under the assumptions C.1-C.4 below. This is useful since the potential outcomes are not always observed, in contrast to the observed data vector $(W,A,\Delta,\tilde{T})$ for each participant, whose distribution we can make direct statistical inferences about.

Define the following assumptions:

C.1****Consistency.

$T=\mathds{1}(A=0)T_{0}+\mathds{1}(A=1)T_{1}$**

C.2****Randomization.

$A$* is independent of $T_{a}$ conditional on $W$, for each $a\in\{0,1\}$*

C.3****Random censoring.

$C$* is independent of $T_{a}$ conditional on $(A,W)$, for each $a\in\{0,1\}$*

C.4****Strong positivity.

$P_{0}(g_{A,0}(a,W)>\epsilon)=1$* and $P_{0}(g_{R,0}(t,a,W)<1-\epsilon)=1$ for each $a\in\{0,1\}$ and $t\in\{0,\ldots,\tau-1\}$ and some $\epsilon>0$.*

We make assumptions C.1-C.4 throughout the manuscript. Denote the survival and censoring function for $T$ at time $t\in\{1,\dots,\tau-1\}$ conditioned on study arm $a$ and baseline variables $w$ by

[TABLE]

Under assumptions C.1-C.4, we have T\mbox{\perp!!!\perp}C\mid A,W and therefore $S(t,a,w)$ and $G(t,a,w)$ have the following product formula representations:

[TABLE]

The potential outcome survival function $S_{\mbox{\footnotesize c}}(t,a,z)$ can be equivalently represented in terms of the observed data distribution as $S(t,a,z)=E\{S(t,a,W)\mid Z=z\}$. It follows from (2) that the causal parameter $\theta_{\mbox{\footnotesize c}}(z)$ is equal to the following observed-data blip function:

[TABLE]

Thus, the value $V_{\mbox{\footnotesize c}}(d)$ of a rule $d$ is equal to $V(d)=E\{d(Z)\theta(Z)\}$, and a corresponding optimal treatment rule is equal to $d_{0}(z)=\mathds{1}\{\theta_{0}(z)>0\}$, where we denote the corresponding true quantities (i.e., quantities computed w.r.t. $P_{0}$) as $\theta_{0}(z)$, $V_{0}(d)$, and $d_{0}(z)$.

In addition to assumptions C.1-C.4 above, we sometimes make the following margin assumption, which is common in the classification literature for plug-in estimators:

C.5****Margin assumption.

There exists a constant $\lambda\geq 0$ such that $P_{0}(0<\theta_{0}(Z)\leq t)\lesssim t^{\lambda}$ for all $t>0$.

The case $\lambda=0$ is trivial and implies no assumption, whereas $\lambda=\infty$ corresponds to the strongest assumption since it implies that $\theta_{0}(Z)$ is bounded away from zero. This assumption characterizes the behavior of the decision function in the boundary, and has been shown crucial to establish the convergence of certain classifiers (e.g., Audibert et al., 2007; Luedtke and Chambaz, 2017).

4 Plug-in Estimation of the Blip Function and the Optimal

Rule

In this section we discuss various estimators for $\theta_{0}(z)$, which can be mapped to a plug-in estimators through $d_{0}(z)=\mathds{1}\{\theta_{0}(z)>0\}$. Our general strategy relies on the concept of censoring unbiased transformation, given in Definition 1 below. This concept was first introduced by Fan and Gijbels (1994) and is further discussed in Rubin and van der Laan (2007), among others.

Definition 1** (Unbiased transformation).**

$D:\mathcal{O}\to\mathbb{R}$* is referred to as an unbiased transformation for $\theta_{0}(z)$ if $E_{0}\left\{D(O)\mid Z=z\right\}=\theta_{0}(z)$.*

The above definition motivates the construction of estimators of $\theta_{0}(z)$ by regressing the transformation $D(O)$ on the covariates $Z$. A common complication in this step is that most unbiased transformations typically depend on unknown nuisance parameters which must be estimated prior to carrying out the analysis. In this work, we focus on the doubly robust censoring unbiased transformation $D_{\eta}$ defined in Lemma 1 below. In addition to being a doubly robust unbiased transformation for $\theta_{0}(z)$ (i.e., providing robustness to inconsistent estimation of one out of two nuisance parameters), $D_{\eta}$ is an efficient estimating function in the non-parametric model in the sense that it may be used to construct efficient estimators of the marginal treatment effect $\theta_{\mbox{\footnotesize c,m}}$ (see e.g., Díaz et al., 2015).

Lemma 1** **(Doubly robust censoring unbiased

transformation).

Define

[TABLE]

where $Z(m,A,W)=Z_{1}(m,A,W)-Z_{0}(m,A,W)$, and

[TABLE]

Assume $\eta=(h,g_{A},g_{R})$ is such that $h=h_{0}$ or $(g_{A},g_{R})=(g_{A,0},g_{R,0})$. Then $D_{\eta}$ is an unbiased transformation for $\theta_{0}(z)$, that is, $E_{0}\left\{D_{\eta}(O)\mid Z=z\right\}=\theta_{0}(z)$ .

As a consequence of the previous lemma, the expected value of the quadratic loss function $L_{\eta}(O;\theta)=\{D_{\eta}(O)-\theta(Z)\}^{2}$ is minimized at $\theta_{0}$ if $\eta=(h,g)$ is such that either $h=h_{0}$ or $g=g_{0}$.

For a loss function $L_{\eta}$, we denote its expected value as $R_{0,\eta}(\theta)=E_{0}\{L_{\eta}(O;\theta)\}$ and refer to it as the risk. We now discuss the construction of super learning ensembles of candidate estimators for $\theta$ that target minimization of the quadratic risk. Consider a collection of estimation algorithms for estimating $\theta_{0}$, hereby called a library, ${\cal L}=\{\hat{\theta}_{j}:j=1,\ldots,J\}$. For an estimator $\hat{\eta}$ of $\eta_{0}$, in light of the discussion of the previous section, this library may be constructed by considering any predictive algorithm that minimizes the quadratic risk for prediction of the doubly robust unbiased transformation $D_{\hat{\eta}}(O)$. The literature in machine and statistical learning provides us with a wealth of algorithms that may be used in this step. Examples include algorithms based on regression trees (e.g., random forests, Bayesian regression trees), algorithms based on smoothing (e.g., generalized additive models, local polynomial regression, multivariate adaptive regression splines), and others (e.g., support vector machines, neural networks).

Consider the following cross-validation set up. Let ${\cal V}_{1},\ldots,{\cal V}_{K}$ denote a random partition of the index set $\{1,\ldots,n\}$ into $K$ validation sets of approximately the same size. That is, ${\cal V}_{k}\subset\{1,\ldots,n\}$; $\bigcup_{k=1}^{K}{\cal V}_{k}=\{1,\ldots,n\}$; and ${\cal V}_{k}\cap{\cal V}_{k^{\prime}}=\emptyset$. In addition, for each $k$, the associated training sample is ${\cal T}_{k}=\{1,\ldots,n\}\setminus{\cal V}_{k}$. Denote $\hat{\eta}_{k}$ the estimator of $\eta_{0}$ trained only using data in $\mathcal{T}_{k}$. Likewise, denote by $\hat{\theta}_{j,k}$ the estimator of $\theta_{0}$ obtained by training the $j$-th predictive algorithm in $\mathcal{L}$ using only data in the sample ${\cal T}_{k}$ (e.g., regressing $D_{\hat{\eta}_{k}}(O_{i})$ on $V_{i}$ for $i\in{\cal T}_{k}$). We use $k(i)$ to denote the index of the validation set that contains observation $i$. The cross-validated prediction risk of $\hat{\theta}_{j}$ is defined as

[TABLE]

In this paper we consider an ensemble learner given by a convex combination

[TABLE]

The weights $\alpha_{j}$ are chosen to minimize the cross-validated risk of the above combination, that is:

[TABLE]

The above expression is a weighted ordinary least squares problem with constraints on the coefficients, and may therefore be solved using standard off-the-shelf regression or optimization software. We denote this super learner with $\hat{\theta}_{\operatorname{sl}}=\hat{\theta}_{\hat{\alpha}}$.

The optimality of general cross-validation selection procedures is discussed in van der Laan & S. Dudoit & A.W. van der Vaart (2006); van der Vaart et al. (2006). Optimality here is defined in terms of asymptotic equivalence with the oracle risk, which we define as the risk computed when (i) one of the components of the nuisance parameter $\eta_{0}$ is known, and (ii) a validation sample of infinite size is available to assess the performance of the estimator. Specifically,

Definition 2** (Oracle risk and oracle selector).**

Let $\eta_{1}=(g_{1},h_{1})$, where either $g_{1}=g_{0}$, or $h_{1}=h_{0}$. The oracle risk of a candidate $\hat{\theta}_{\alpha}$ is defined as

[TABLE]

The oracle selector is equal to

[TABLE]

and the corresponding oracle blip function is denoted with $\hat{\theta}_{\operatorname{or}}=\hat{\theta}_{\tilde{\alpha}}$.

The risk $\tilde{R}_{\eta_{1}}(\theta_{0})=\int L_{\eta_{1}}(o;\theta_{0})dP_{0}(o)$ is the optimal risk (with respect to the loss function $L_{\eta_{1}}$, which in light of Lemma 1 is a valid loss function) achieved by the true $\theta_{0}$. The following theorem provides a bound on the excess risk of the estimator $\hat{\theta}_{\operatorname{sl}}$ and the excess risk of $\hat{\theta}_{\operatorname{or}}$. The excess risk for a selector $\hat{\alpha}$ is defined as the difference between the oracle risk of the selector $\hat{\alpha}$ and the optimal risk, i.e.,

[TABLE]

where we remind the reader that the expectation is taken over draws of $O_{1},\ldots,O_{n}$. We denote this excess risk as $\operatorname{\mathcal{E}}^{2}(\hat{\theta}_{\hat{\alpha}})$, below we refer to its square root as $\operatorname{\mathcal{E}}(\hat{\theta}_{\hat{\alpha}})$. We show that the above excess risk is bounded by two terms: one depending on the excess risk of the oracle selector $\operatorname{\mathcal{E}}(\hat{\theta}_{\operatorname{or}})$, and another one depending on doubly robust terms associated to estimation of $\eta_{0}$.

Theorem 1** **(Oracle inequality for the super learner of the blip

function).

Let $\eta_{1}=(g_{1},h_{1})$ denote the element-wise $L_{2}(P_{0})$ limit of $\hat{\eta}$ as $n\to\infty$, and assume that either $g_{1}=g_{0}$ or $h_{1}=h_{0}$. Define

[TABLE]

Then, for $\delta>0$

[TABLE]

*for constants $C_{1}$, $C_{2}$, and $C_{3}$. *

Note that the terms $B_{1}(\hat{\eta},\eta_{0})$ and $B_{2}(\hat{\eta},\eta_{0})$ converge to zero if either $\hat{g}$ or $\hat{h}$ converge to $g_{0}$ or $h_{0}$, respectively, in $L_{2}(P_{0})$ norm. This implies that the doubly robust property of $D_{\eta}(o)$ is transferred to the oracle inequality. To the best of our knowledge this result had not been previously shown in the literature.

The super learner $\hat{\theta}_{\operatorname{sl}}$ may be used to construct a plug-in estimator of the optimal rule as $\hat{d}(z)=\mathds{1}\{\hat{\theta}_{\operatorname{sl}}(z)>0\}$. The following remark discusses the convergence rates of the value of the selected rule $V_{0}(\hat{d})$ to the value of the oracle rule $V_{0}(\tilde{d})$.

Remark 1** (Convergence rates to the oracle value).**

Assume

[TABLE]

Lemma 5.3 of Audibert et al. (2007), along with Jensen’s inequality, show that under assumption C.5, we have

[TABLE]

where $\tilde{d}(z)=\mathds{1}\{\hat{\theta}_{\operatorname{or}}(z)>0\}$ is the oracle rule. This yields the following convergence rate:

[TABLE]

An example of a case yielding the above rate is a randomized study ($g_{A,0}(w)=q\in(0,1)$) with no censoring ($P_{0}(\Delta=1)=1$.) In this case, a logistic regression fit of $A$ on $W$ containing at least an intercept would yield an estimator satisfying $||\hat{g}_{A}-g_{A,0}||^{2}=O_{P}(n^{-1})$. Plugging in the true value $g_{R,0}(t,a,w)=0$ for $\hat{g}_{R}(t,a,w)$, and assuming $\hat{h}$ is inconsistently estimated yields $B_{1}(\hat{\eta},\eta_{0})=O(n^{-1/2})$ and $B_{2}(\hat{\eta},\eta_{0})=O(n^{-1/2})$. Under no margin assumption ($\lambda=0$) we get a convergence rate of $(\log n/n)^{1/2}$. Under a strong margin assumption in which $\theta_{0}(Z)$ is bounded away from zero ($\lambda=\infty$) we get a rate of $\log n/n$.

The above convergence result establishes the convergence of the value of our estimator $\hat{d}$ to the value of the oracle $\tilde{d}$. This is different from the typical result in the classification literature, which establishes convergence to the optimal value $V_{0}(d_{0})$. The latter result often involves fast learning rates (sometimes faster than $n^{-1}$) and requires restricting the class of blip functions considered to Hölder (Audibert et al., 2007) or Donsker (Luedtke and Chambaz, 2017) classes, a restriction we do not impose.

5 Super Learner Ensembles for the Optimal Rule from a

Classification Perspective

5.1 Estimators Using the 0-1 Loss Function

In this section we discuss a classification approach that aims at directly estimating the optimal rule $d_{0}(z)$. Our approach here differs from the previous section in that we do not attempt to estimate the blip function. Instead, we introduce the concept of a decision function, defined as $f:\mathbf{Z}\to\mathbb{R}$, and which yields a treatment rule $d_{f}(z)=\mathds{1}\{f(z)>0\}$. In a slight abuse of notation we use $V(f)$ to refer to the value of the rule $d_{f}$. Any function $f_{0}$ such that $\operatorname{sign}\{f_{0}(z)\theta_{0}(z)\}=1$ has optimal value $V_{0}(d_{0})$. This provides intuition on the benefits of directly optimizing the value of the loss function instead of the risk of the blip function: an inconsistent estimator of the blip function may provide an optimal rule, as long as its sign is correct. For a given rule $d_{f}$, in light of Lemma 1, we have that $V_{0}(f)=E_{0}\{d_{f}(Z)D_{\eta}(O)\}$ if $\eta$ is such that either $h=h_{0}$, or $g=g_{0}$. Thus, a decision function that optimizes the value of the rule $d_{f}$ may be found as

[TABLE]

For a binary value $b\in\{0,1\}$ and any $X$ we have $bX=\mathds{1}\{X>0\}|X|-|X|\mathds{1}\left[\mathds{1}\{X>0\}\neq b\right]$. Thus, the optimization problem may be recast as $f_{0}\in\mathcal{F}_{0}$, where $\mathcal{F}_{0}=\operatorname*{arg\,min}_{f}\int L_{\eta}(o;f)dP_{0}$ and

[TABLE]

Expression (8) is a weighted classification loss function in which we aim to classify the binary outcome $\mathds{1}\{D_{\eta}(O)>0\}$ based on data $Z$, using the 0-1 loss function with weights given by $|D_{\eta}(O)|$. The objective is to classify an individual who benefits from treatment arm $A=1$ (i.e., an individual with $D_{\eta}(O)>0$) as requiring treatment (i.e., $d_{f}(Z)=1$), while penalizing for the loss $|D_{\eta}(O)|$ incurred if the individual were misclassified.

In what follows we consider a library of algorithms for estimation of the decision function $\mathcal{L}=\{\hat{f}_{j}(z):j,\ldots,J\}$. In light of the discussion of the previous sections, the most natural choice for a decision function is the blip function $\hat{\theta}(z)$. However, we do not restrict our setup to functions with a blip interpretation. In addition to estimators of the blip function $\theta_{0}(z)$, the library may contain other decision functions such as the support vector machines proposed by Zhao et al. (2015) and the parametric decision functions of Bai et al. (2016).

We construct an ensemble of the decision functions as

[TABLE]

In this way, we generate an ensemble optimal rule as $\hat{d}_{\alpha}(z)=\mathds{1}\{\hat{f}_{\alpha}(z)>0\}$. As in the previous section, we define the super learner selector as

[TABLE]

where $\hat{f}_{\alpha,k(i)}$ represents (9) with $\hat{f}_{j}(z)$ replaced by $\hat{f}_{j,k(i)}(z)$: the $j$-th decision function estimated using the training sample $\mathcal{T}_{k(i)}$. The super learner of the decision function is defined as $\hat{f}_{\operatorname{sl}}(z)=\hat{f}_{\hat{\alpha}}(z)$, and the corresponding optimal rule is defined as $\hat{d}_{\operatorname{sl}}=\mathds{1}\{\hat{f}_{\operatorname{sl}}(z)>0\}$.

For $\eta_{1}=(g_{1},h_{1})$ such that either $g_{1}=g_{0}$ or $h_{1}=h_{0}$, the oracle risk of the decision function is defined as

[TABLE]

The oracle selector of $\alpha$ is thus defined as $\tilde{\alpha}\in\operatorname*{arg\,min}_{\alpha}\tilde{R}_{\eta_{1}}(\hat{f}_{\alpha})$, and we denote $\hat{f}_{\operatorname{or}}=\hat{f}_{\tilde{\alpha}}$. The excess risk of an estimator $\hat{f}$ is equal to

[TABLE]

In Theorem 2 below, we provide bounds on $\operatorname{\mathcal{E}}(\hat{f}_{\operatorname{sl}})$ in terms of the excess risk of the oracle selector $\operatorname{\mathcal{E}}(\hat{f}_{\operatorname{or}})$ and the bias terms $B_{1}(\hat{\eta},\eta_{0})$ and $B_{2}(\hat{\eta},\eta_{0})$ defined in Theorem 1.

Theorem 2** **(Oracle inequality for the super learner of the optimal

rule).

Assume the conditions of Theorem 1. In addition, assume that $\hat{\alpha}$ is computed in a grid of size $Mn^{q}$ for some $M>0$, $q>0$. Then,

[TABLE]

If condition C.5 holds with $\lambda=\infty$, then, for $\delta>0$

[TABLE]

for constants $C_{1}$, $C_{2}$, and $C_{3}$, where $B_{1}$ and $B_{2}$ are defined as in Theorem 1.

Remark 2**.**

Assume $B_{1}(\hat{\eta},\eta_{0})$ converges as in Remark 1. An immediate consequence of the above result is that under no margin assumption, we have $\mathbb{E}\{V_{0}(\hat{d}_{\operatorname{or}})-V_{0}(\hat{d}_{\operatorname{sl}})\}=O\left((\log n/n)^{1/2}\right)$. Under the strong margin assumption C.5 with $\lambda=\infty$ we have $\mathbb{E}\{V_{0}(\hat{d}_{\operatorname{or}})-V_{0}(\hat{d}_{\operatorname{sl}})\}=O(\log n/n)$. These rates are identical to the rates obtained in Remark 1 for the plug-in estimator. The question of whether analogous convergence rates may be obtained for other values of $\lambda$ under condition C.5 remains an open problem.

In comparison to Theorem 1, Theorem 2 has the additional assumption that the optimization of the loss function is carried out in a grid polynomial size in $n$. Inspection of the proofs of the theorems in the Supplementary Material reveals the reason for the additional assumption: the 0-1 loss function is non-smooth and the Lipschitz condition used in the proof of Theorem 1 does not apply. As demonstrated in our data application, this assumption is likely to have little practical consequences, but it is unclear to us whether it can be removed.

5.2 Using a Surrogate Loss Function for the 0-1 Loss

It is well known in the statistical learning literature that minimizing (8) is generally difficult due to the discontinuity and non-convexity of the 0-1 loss. A common approach to mitigate the issues arising from the discontinuity and non-convexity of the 0-1 loss function is to use surrogates loss functions, such as the logistic loss $\phi(x)=\log(1+\exp(-x))$ or the hinge loss $\phi(x)=\max(1-x,0)$. We have the following result, which teaches us that any decision function $d_{f_{0}}(z)$ based on a decision function $f_{0}\in\mathcal{F}_{\phi,0}$ has the same performance as the optimal rule $d_{0}(z)$.

Lemma 2**.**

Assume $\eta$ is such that either $h=h_{0}$, or $g=g_{0}$. Define

[TABLE]

where the surrogate loss $L_{\phi,\eta}$ is defined as

[TABLE]

Then we have $\mathcal{F}_{\phi,0}\subseteq\mathcal{F}_{0}$, where $\mathcal{F}_{0}=\operatorname*{arg\,min}_{f}\int L_{\eta}(o;f)dP_{0}$.

6 Estimating the Optimal Treatment for Breast

Cancer Patients in our Motivating Application

Different types of human breast cancer tumors have been shown to have heterogeneous response to treatments (Perou et al., 2000; Sotiriou and Pusztai, 2009). Amplification of ERBB2 gene and associated overexpression of human epidermal growth factor receptor (HER2) encoded by this gene occur in 25-30$\%$ of breast cancers (Slamon et al., 2001). HER2-positive breast cancer is an aggressive form of the disease and the prognosis for such patients is generally poor (Slamon, 1987; Seshadri et al., 1993). The clinical efficacy of adjuvant trastuzumab, a recombinant monoclonal antibody, in early stage HER2-positive patients was demonstrated by several large clinical trials (Perez et al., 2011; Romond et al., 2005). Despite significant improvement in disease-free and overall survival of patients treated with trastuzumab, about 20-25$\%$ patients relapse within 3-5 years (Perez et al., 2011). In this paper we use data from the North Central Cancer Treatment Group N9831 study, a phase III randomized clinical trial testing the addition of trastuzumab to chemotherapy in stage I-III HER2-positive breast cancer.

The total number of patients enrolled in the NCCTG N9831 trial was 3,505. Samples from 1,390 patients, for whom there was available tissue, were used to quantify mRNA from a custom codeset of 730 genes created by experts. The available baseline variables may be thus be categorized in three classes: demographic (e.g, race, age, ethnicity), clinical (e.g., tumor grade, tumors size, nodal status, hormone receptor status), and gene expression. Among the 1,390 patients, 483 received chemotherapy alone (control arm) and 907 patients received chemotherapy plus trastuzumab (treatment arm).

The clinical challenge is to identify genetic and demographic profiles for patients with HER2-positive breast cancer who are unlikely to benefit from adjuvant trastuzumab.

In order to estimate and assess the performance of the estimated rule using different datasets, we split our data into training and validation datasets, of sizes 1000 and 390, respectively.

6.1 Estimators of $h$ and $g_{R}$

According to our theoretical results, the optimality of the estimated treatment rules hinges upon consistent estimation of at least one of the nuisance parameters $h$ and $g_{R}$. As a result, it is crucial to employ flexible methods capable of unveiling complex patterns which are not visible to the human eye. As demonstrated below in Section 6.3, simple parsimonious models such as the Cox proportional hazards or logistic regression fail to detect these complex relations in the data.

In order to accurately estimate the nuisance parameters, we use an ensemble learner known as the super learner for prediction (van der Laan et al., 2007). We train the ensemble separately using data from each treatment arm, in order to fully account for treatment-covariate interactions. Like our rule ensembles, super learning predictors build a combination of candidate predictors that minimize a cross-validated user-supplied risk function. Since $g_{R}$ and $h$ are conditional probabilities, we focus on logistic regression ensembles and the negative log-likelihood loss function, using the R implementation in the SuperLearner package (Polley et al., 2016). The candidate estimators included in the ensembles are listed in Table 1, along with the coefficients of each predictor in the ensemble, when trained in the complete dataset using 5-fold cross-validation.

For random forests, extreme gradient boosting, and multilayer perceptron, the tuning parameters are tuned using data splitting with the aid of the R caret package (Kuhn et al., 2016). To avoid $p>n$, logistic regression and multivariate adaptive splies are estimated with a variable screening algorithm which computes univariate t-statistics and keeps only the 50 variables with a larger value.

6.2 Candidate Estimators for the Optimal Treatment Rule

According to our discussion in Sections 4 and 5, there are at least three types of estimators for the optimal rule $d_{0}$. The first type is a simple substitution estimator, obtained through inspection of equation (4), which consists in regressing the blip function $\hat{B}(W)=\sum_{t=1}^{\tau-1}\{\hat{S}(t,1,W)-\hat{S}(t,0,W)\}$ on $Z$, where $\hat{S}$ is the estimator of the survival function corresponding to the estimator $\hat{h}$ described in Section 6.1. The second type is obtained through regression of the unbiased transformation $D_{\hat{\eta}}(O)$ on $Z$. The third type of estimation methods is obtained based on equation (8), and is obtained by classifying the binary outcome $\mathds{1}\{D_{\hat{\eta}}(O)>0\}$ as a function of $Z$, with weights given by $|D_{\hat{\eta}}(O)|$. Here, $\hat{\eta}=(\hat{g}_{R},\hat{h})$, where the components of $\hat{\eta}$ are as described in Section 6.1. Any regression or supervised classification technique available in the statistical learning literature may be used as a candidate for solving these problems.

In our application, we focus on the following candidates for estimating $d_{0}$:

[TABLE]

According to our discussion in Sections 4 and 5, we also train four super learning ensembles of the above candidate estimators, using different loss functions:

[TABLE]

The coefficients of each candidate estimator in each ensemble are presented in Table 2. These coefficients were computed using the Subplex (Rowan, 1990) routine implemented in the NLopt nonlinear-optimization R package. For improved robustness, the 0-1 loss was optimized using 1000 different random starting values.

6.3 Assessing the Performance of The Estimated Treatment Rule

Once each rule is estimated using only data in the training dataset, its value $V(\hat{d})$ is estimated on the validation dataset. To that effect, we use the targeted minimum loss based estimator of the restricted mean survival time proposed by Díaz et al. (2015) (See also Moore and van der Laan, 2011).

Figure 1 presents the estimated restricted mean survival time obtained with each estimated rule, along with 95% confidence intervals. For comparison, we also present the value of two static rules of interest: never treat and always treat. As is clear from the figure, the best algorithm in our application is regression of the blip function. All super learning ensembles yield a similar value, demonstrating the oracle property of the super learner. Treating patients according to the optimal rule yields a restricted mean survival of 157.1 (s.d. 3.1) months. In comparison with the always treat rule, which yields 151.2 (s.d., 3.3) months, the optimal rule improves mean patient survival by 6 months.

According to Table 2, only the super learning ensemble based on the 0-1 loss assigns a large weight to the best algorithm. In fact, its restricted mean survival time (see Figure 1) is identical to that of the optimal rule. The other ensembles assign more weight to the second best algorithm, weighted classification using extreme gradient boosting, and have slightly smaller restricted mean survival time. This is in agreement with our theoretical findings that the best performance is obtained using the 0-1 loss function.

It is also worth noting that three of the estimated rules (weighted classification using random forests and logistic regression, and regression of the function $D_{\hat{\eta}}$) yield a restricted mean survival time smaller or equal than the restricted mean survival time of the static rule always treat.

Table 3 in the Supplementary Materials shows the p-value for the pair-wise comparisons of the value of each estimated rule. A few interesting points to note are:

1. (i)

The ensemble using the 0-1 loss function outperforms the other ensembles, the difference in values is small but significant at 5% level. 2. (ii)

The value of the optimal rule, which is obtained through regression of the blip function (see Figure 1), is significantly different from all other rules, except the ensemble using the 0-1 loss function. This illustrates the theoretical property of super-learning stating that the risk of the ensemble converges to the risk of the best candidate in the library. 3. (iii)

Weighted classification using logistic regression, which is often advocated because it yields parsimonious rules (e.g., Zhang et al., 2015), has a value significantly lower than the static rule always treat.

In our application, we have decided to use data splitting to train and assess the performance of the estimated rules. Though correct, this approach may be unnecessary, since the value of the rule may be assessed using the training dataset, under certain conditions derived by Luedtke et al. (2016).

7 Discussion

We present two methods for constructing an ensemble individualized treatment rule. The methods are based on a plug-in estimator optimizing the prediction error of the blip function, and a weighted classification approach which directly estimates the decision function. Though we found no theoretical differences between the two approaches in terms of their asymptotic properties, the classification ensemble using the 0-1 loss function yielded better treatment rules than the other approaches in our illustrative application. The superiority of the classification approach has been recognized before (e.g., Zhao et al., 2012), and is a consequence of the fact that it emphasizes optimizing the decision rule rather than prediction accuracy emphasized by the blip approach.

We consider a survival time measured in a discrete time scale. Most clinical research studies measure time to event in a discrete scale. In our motivating application, time to relapse of cancer or death was measured in days. We foresee no technical difficulties in extending our approach to consider a continuous time to event. This can be achieved by replacing discrete time hazards by their continuous counterpart, as well as replacing certain sums over time by the appropriate martingale integrals (see Bai et al., 2016) in the definition of the censoring unbiased transformation $D_{\eta}(O)$. A potential practical limitation is that the software and literature for data-adaptive machine learning estimation of continuous time hazards (required for the nuisance parameters) is scarce in comparison to that of binary classification, which may be used for estimation of discrete time hazards. Among the few methods that can be used for this problem are (semi)-parametric models such as Cox regression and accelerated failure time models. Available data adaptive approaches include survival random forests and regularized Cox regression. If time is measured on a continuous scale, implementation of our methods requires discretization. The specific choice of the discretization intervals may be guided by what is clinically relevant. For example, in cancer research, the clinically relevant scale would typically be a week or a month. In the absence of clinical criteria to guide the choice of discretization level, a concern is that too coarse of a discretization may lead to relevant information loss. A question for future research is how to optimally set the level of discretization in order to trade off information loss versus estimator precision. Another area for future research is to consider discretization levels that get finer with sample size.

We present doubly robust oracle inequalities and convergence rates assuming (i) an exposure of interest that occurs at baseline, and (ii) censoring which is confounded with the time to event only by baseline variables. We conjecture that our general results apply to the more general case of a dynamic treatment regime with a time-varying treatment and time-varying confounders. Such results will be the subject of future research.

In our definition of oracle risk and oracle selector we have used a nuisance parameter $\eta_{1}=(g_{1},h_{1})$ satisfying either $h_{1}=h_{0}$, or $g_{1}=g_{0}$. Therefore these oracle quantities change depending on which of the two nuisance parameters is correctly specified. According to efficient estimation theory in semi-parametric models, we expect the case $\eta_{1}=\eta_{0}$ to yield oracle quantities with minimal variability. In the single misspecification case in which $h_{1}=h_{0}$ or $g_{1}=g_{0}$ but not both, it is unclear to us whether misspecification of one of the models yields better results than the other. Lastly, our setup includes as particular case the inverse probability weighted loss function, which may be obtained by using a constant estimator $\hat{h}(t,a,w)=1$, as well as the g-computation loss function, which is obtained by using $\hat{g}_{A}(a,w)=1$ and $\hat{g}_{R}(t,a,w)=0$.

8 Supplementary Material

8.1 Motivating application

8.2 Proofs of Theorems and Lemmas

8.2.1 Lemma 1

*Proof * For simplicity, consider the treatment-time-specific function

[TABLE]

and note that $D_{\eta}=\sum_{m=1}^{\tau-1}(D_{m,1,\eta}-D_{m,0,\eta})$. For a function $f(t,a,w)$ we denote $Pf(t)=\int f(t,a,w)dP(w)$. Conditioning first on $W$ in the above display yields

[TABLE]

Thus, we have

[TABLE]

Plugging in $g=g_{0}$ or $h=h_{0}$ yields the result. ∎

8.2.2 Theorem 1

*Proof * We start by assuming the minimization of the risk in the definition of $\hat{\alpha}$ and $\tilde{\alpha}$ is carried out in a grid $\mathcal{B}_{n}\subset\mathcal{B}=\{\alpha\in\mathbb{R}^{J}:\alpha_{j}\geq 0,\sum_{j=1}^{J}\alpha_{j}=1\}$ of polynomial size in $n$ (that is $|\mathcal{B}_{n}|\lesssim n^{q}$) for some $1\leq q<\infty$, but do away with this assumption at the end of the proof. Let $\hat{\beta}$ ad $\tilde{\beta}$ denote the cross-validated and oracle selectors when the risk minimization is performed in $\mathcal{B}_{n}$ rather than $\mathcal{B}$. We use $P_{n,k}$ to denote the empirical distribution corresponding to the validation set $\mathcal{V}_{k}$, as well as $E_{K}(X)=K^{-1}\sum_{k=1}^{K}X_{k}$ to denote an average across validation splits. We denote $PL(\theta)=\int L(o;\theta)dP(o)$. Let

[TABLE]

Define the centered loss function

[TABLE]

In this proof we denote with $R$ the corresponding centered risks, i.e., denote

[TABLE]

the corresponding cross-validated and oracle risks. For notational convenience we denote $R(\beta)=R(\hat{\theta}_{\beta})$. Note that $\tilde{R}_{\eta^{\star}}(\beta)=\operatorname{\mathcal{E}}(\hat{\theta}_{\beta})$. For $\delta>0$ we have

[TABLE]

where the second inequality is a consequence of the definition of $\hat{\beta}$ as the minimizer of $\hat{R}_{\hat{\eta}}(\beta)$, and the last equality is the result of adding and subtracting some terms. Denote (11) with $T$, (12) with $H$, and (13) with $Q(\tilde{\alpha})$.

Note that the assumptions of the theorem imply that $P_{0}\{|D_{\eta^{\star}}(O)|\leq M\}=1$ for some constant $M$. This, together with Lemma 5 below, allow the application of Lemma 3 in van der Laan and Dudoit (2003) (see also pages 143-145 of Dudoit and van der Laan, 2005) to show that

[TABLE]

It remains to analyze $Q(\hat{\beta})$ and $Q(\tilde{\beta})$. First, we write $Q_{\beta}=Q_{1}(\beta)+Q_{2}(\beta)$, where

[TABLE]

For $Q_{1}(\beta)$, note that $R_{\eta^{\star}}(\beta)=E_{K}P_{0}(\theta_{0}-\theta_{\beta,k})^{2}$. This yields, for $\beta\in(\hat{\beta},\tilde{\beta})$,

[TABLE]

Conditioning on $W$ first, from the definition of $\eta^{\star}$, Lemma 1 shows that the second term in the right hand side is zero. Conditioning on $W$ first along with the proof of Lemma 1 and the Cauchy-Schwartz inequality also yields

[TABLE]

where the last inequality follows from Lemma 5 in Appendix 8.2.5 and the definition of $\tilde{\beta}$ as the minimizer of $\tilde{R}_{\eta^{\star}}(\beta)$, and the second to last inequality follows from Cauchy-Schwartz applied to the norm defined by the inner product $<f_{k},g_{k}>=\mathbb{E}\,E_{K}P_{0}f_{k}g_{k}$. For $Q_{2}(\beta)$, note that $(P_{n,k}-P_{0})(L_{\hat{\eta}_{k}}-L_{\eta^{\star}})(\theta_{\beta,k})$ is an empirical processes with index set $\mathcal{B}_{n}$, where the latter set is finite. We will apply the following inequality for empirical processes with finite index set:

[TABLE]

where $F$ is an envelope of $\mathcal{F}$. This result is a direct consequence of Lemma 19.38 of van der Vaart (1998). Note that the all functions in $\mathcal{F}_{k}=\{(L^{0}_{\hat{\eta}_{k}}-L^{0}_{\eta^{\star}})(\theta_{\beta,k}):\beta\in\mathcal{B}_{n}\}$ satisfy

[TABLE]

where the second inequality follows from Lemma 5. Thus, the envelope $F_{k}$ of $\mathcal{F}_{k}$ is bounded by the same quantity. This, together with (14) shows

[TABLE]

This proves

[TABLE]

which is equivalent to $x^{2}-bx\leq c$ for

[TABLE]

The quadratic formula $x\leq(b+\sqrt{b^{2}+4c})/2$ implies $x\leq b+\sqrt{c}$, which yields

[TABLE]

From our definitions and assumptions, the function $f(\beta)=R_{\eta_{0}}(\hat{\theta}_{\beta})$ satisfies the Lipschitz condition

[TABLE]

where $||\cdot||_{\infty}$ denotes the supremum norm and $||\cdot||_{2}$ the Euclidean norm. Thus $f(\hat{\beta})-f(\hat{\alpha})$ and $f(\tilde{\beta})-f(\tilde{\alpha})$ are both bounded by $n^{-q}$, which allows us to replace $(\hat{\beta},\tilde{\beta})$ by $(\hat{\alpha},\tilde{\alpha})$ in (15), completing the proof of the theorem. ∎

8.2.3 Theorem 2

*Proof *

For convenience in the calculations we use the loss function

[TABLE]

which is equivalent to the one used in the Theorem. Let $\eta^{\star}$, $L_{\eta}^{0}(O;f)$, $\hat{R}_{\eta}(\beta)$, and $\tilde{R}_{\eta}(\beta)$ be defined as in the proof of Theorem 1. We have

[TABLE]

Define

[TABLE]

Since, by definition, $\hat{R}_{\hat{\eta}}(\hat{\beta})\leq\hat{R}_{\hat{\eta}}(\tilde{\beta})$, we have

[TABLE]

van der Laan and Dudoit (2003), page 26, show that

[TABLE]

In the proof of Theorem 2, we show that

[TABLE]

completing the proof of first claim of the theorem.

Assume now condition C.5 holds with $\alpha=\infty$ such that $\inf_{z\in\mathbf{Z}}|\theta_{0}(z)|>0$ The proof in this case has the same steps as the proof of Theorem 1 and we will only provide a sketch. The conditions of the Theorem allow application of Lemma 4 below to obtain

[TABLE]

For $Q_{1}(\alpha)$ and $Q_{2}(\alpha)$ we get

[TABLE]

For $Q_{2}(\alpha)$, note that $(P_{n,k}-P_{0})(L_{\hat{\eta}_{k}}-L_{\eta^{\star}})(f_{\alpha,k})$ is an empirical processes with index set $\mathcal{A}_{n}$, where the latter set is the finite set with $Mn^{q}$ points in which $\hat{\alpha}$ is computed. We will apply inequality (14). Note that the all functions in $\mathcal{F}_{k}=\{(L^{0}_{\hat{\eta}_{k}}-L^{0}_{\eta^{\star}})(f_{\alpha,k}):\alpha\in\mathcal{A}_{n}\}$ satisfy

[TABLE]

where the second inequality follows from Lemma 5. Thus, the envelope $F_{k}$ of $\mathcal{F}_{k}$ is bounded by the same quantity. This, together with (14) shows

[TABLE]

This completes the proof. ∎

8.2.4 Lemma 2

*Proof * This is a direct application of Theorems 1 (part 3) and 2 of Bartlett et al. (2006). See also Theorem 5 of Luedtke and van der Laan (2016). ∎

8.2.5 Lemmas

Lemma 3**.**

Consider the assumptions of Theorem 1. Let $Z=L_{\eta_{0}}(O;\theta)-L_{\eta_{0}}(O;\theta_{0})$. We have

[TABLE]

*Proof * First, note that

[TABLE]

In light of Lemma 1 we have

[TABLE]

Note that $P_{0}\{|2D_{\eta_{0}}-\theta(Z)-\theta_{0}(Z)|\leq 4\max(M,C_{1})\}=1$. Thus

[TABLE]

which completes the proof of the lemma. ∎

Lemma 4**.**

Consider the assumptions of Theorem 2. Let

[TABLE]

Let $Z=L_{\eta_{1}}(O;\theta)-L_{\eta_{1}}(O;\theta_{0})$. We have

[TABLE]

*Proof * We have

[TABLE]

∎

Lemma 5**.**

For each $\hat{\eta}=(\hat{g},\hat{h})\to\eta_{1}=(g_{1},h_{1})$ such that either $g_{1}=g_{0}$ or $h_{1}=h_{0}$ define

[TABLE]

We have

[TABLE]

with $B^{2}$ defined in Theorem 1.

*Proof * First let $g_{1}=g_{0}$ and $h_{1}\neq h_{0}$. Then $\eta^{\star}=(g_{0},\hat{h})$ and straightforward algebra shows

[TABLE]

Analogously, for $g_{1}\neq g_{0}$ and $h_{1}=h_{0}$ we have

[TABLE]

Now, for $g_{1}=g_{0}$ and $h_{1}=h_{0}$ we get

[TABLE]

Putting these results together proves the lemma. ∎

Lemma 6**.**

For two sequences $a_{1},\ldots,a_{m}$ and $b_{1},\ldots,b_{m}$ we have

[TABLE]

*Proof * Replace $(b_{t}-a_{t})$ by $(1-a_{t})-(1-b_{t})$ in the right hand side and expand the sum to notice it is a telescoping sum. ∎