Deep Learning for Survival Outcomes

TL;DR

This paper introduces a novel deep learning framework tailored for censored survival data, enabling accurate estimation of survival probabilities and mean survival times, with demonstrated strong performance on simulated and real datasets.

Contribution

It develops a new class of deep learning algorithms that handle censored outcomes by using a censoring unbiased transformation, extending deep learning methods to survival analysis.

Findings

Strong performance in simulations

Effective estimation of survival probabilities

Accurate restricted mean survival estimates

Abstract

This manuscripts develops a new class of deep learning algorithms for outcomes that are potentially censored. To account for censoring, the unobservable loss function used in the absence of censoring is replaced by a censoring unbiased transformation. The resulting class of algorithms can be used to estimate both survival probabilities and restricted mean survival. We show how the deep learning algorithms can be implemented using software for uncensored data using a form of response transformation. Simulations and analysis of the Netherlands 70 Gene Signature Data show strong performance of the proposed algorithms.

Equations

$$\frac{1}{n}\sum _{i=1}^{n}L(h(T_i),f_{\beta }(W_i))+\eta ||\beta |{|}_c^2.$$

$$\alpha (m)=\sum _{l=1}^D\sum _{i=1}^n{A}_{i,l}L\left(h(T_i),f_{{\hat{\beta }}_{{\eta }_m}^{(l)}}(W_i)\right).$$

$$\frac{1}{n}\sum _{i=1}^n\frac{{\Delta }_{i}L(h(T_i),\beta (W_i))}{G_0({\tilde{T}}_i|W_i)}.$$

$$L_{IPCW}(\mathcal{O},\beta ;\hat{G})=\frac{1}{n}\sum _{i=1}^n\frac{{\Delta }_{i}L(h(T_i),\beta (W_i))}{\hat{G}({\tilde{T}}_i|W_i)}.$$

$$L_{IPCW}(\mathcal{O},\beta ;G_0)+\frac{1}{n}\sum _{i=1}^n\left(\frac{(1-{\Delta }_i)m_L({\tilde{T}}_i,W_i;S_0)}{G_0({\tilde{T}}_i|W_i)}-{\int }_0^{{\tilde{T}}_i}\frac{m_L(u,W_i;S_0)}{G_0(u|W_i)}d{\Lambda }_{G_0}(u|W_i)\right),$$

$$m_L(u,w;S)=E_S[L(h(T),W)|T\ge u,W=w]=-{\int }_u^{\infty }\frac{L(h(t),w)}{S(u|w)}dS(t|w).$$

$$L_{DR}(\mathcal{O},\beta ;\hat{G},\hat{S})=L_{IPCW}(\mathcal{O},\beta ;\hat{G})+L_{Aug}(\mathcal{O},\beta ;\hat{G},\hat{S}),$$

$$L_{Aug}(\mathcal{O},\beta ;\hat{G},\hat{S})=\frac{1}{n}\sum _{i=1}^n\left(\frac{(1-{\Delta }_i)m_L({\tilde{T}}_i,W_i;\hat{S})}{\hat{G}({\tilde{T}}_i|W_i)}-{\int }_0^{{\tilde{T}}_i}\frac{m_L(u,W_i;\hat{S})}{\hat{G}(u|W_i)}d{\Lambda }_{\hat{G}}(u|W_i)\right).$$

$$\frac{1}{n}\sum_{i=1}^{n}\big{(}\Delta_{i}L(h(T_{i}),\beta(W_{i}))+(1-\Delta_{i})m_{L}(C_{i},W_{i};S_{0})\big{)}.$$

$$L_{BJ}(\mathcal{O},\beta ;\hat{S})=\frac{1}{n}\sum _{i=1}^n\left({\Delta }_{i}L(h(T_i),\beta (W_i))+(1-{\Delta }_i)m_L(C_i,W_i;\hat{S})\right).$$

$$\mathcal{O}(\tau )=\{{\tilde{T}}_i(\tau )=\min (T_i,C_i,\tau ),{\Delta }_i(\tau )=I(\min (T_i,\tau )\le C_i),W_i;i=1,\dots ,n\}$$

$$\frac{1}{n}\sum _{i=1}^n\frac{{\Delta }_i(I(T_i\ge t)-\beta (W_i){)}^2}{\hat{G}({\tilde{T}}_i|W_i)}.$$

$$L_{IPCW,t}(\mathcal{O}(t),\beta ;\hat{G})=\frac{1}{n}\sum _{i=1}^n\frac{{\Delta }_i(t)(I({\tilde{T}}_i\ge t)-\beta (W_i){)}^2}{\hat{G}({\tilde{T}}_i(t)|W_i)}.$$

$$L_{DR,t}(\mathcal{O}(t),\beta ;\hat{G},\hat{S})=L_{IPCW,t}(\mathcal{O}(t),\beta ;\hat{G})+L_{Aug,t}(\mathcal{O}(t),\beta ;\hat{G},\hat{S}),$$

$$L_{Aug,t}(\mathcal{O}(t),\beta ;\hat{G},\hat{S})$$

$$=\frac{1}{n}\sum _{i=1}^n\frac{(1-{\Delta }_i(t))m_{L_2,t}(C_i,W_i;\hat{S})}{\hat{G}({\tilde{T}}_i(t)|W_i)}+{\int }_0^{{\tilde{T}}_i(t)}\frac{m_{L_2,t}(u,W_i;\hat{S})d{\hat{\Lambda }}_G(u|W_i)}{\hat{G}(u|W_i)},$$

$$L_{BJ,t}(\mathcal{O}(t),\beta ;\hat{S})=\frac{1}{n}\sum _{i=1}^n\left({\Delta }_i(t)(I({\tilde{T}}_i\ge t)-\beta (W_i){)}^2+(1-{\Delta }_i(t))m_{L_2,t}(C_i,W_i;\hat{S})\right).$$

$$A_{ki}(G)=\frac{{\Delta }_{i}h(\widetilde{T_i}{)}^k}{G({\tilde{T}}_i|W_i)},B_{ki}(G,S)=\frac{(1-{\Delta }_i)m_k({\tilde{T}}_i,W_i;S)}{G({\tilde{T}}_i|W_i)},$$

$$C_{ki}(G,S)={\int }_0^{{\tilde{T}}_i}\frac{m_k(u,W_i;S)d{\Lambda }_G(u|W_i)}{G(u|W_i)}.$$

$$m_k(t,w;S)=E_S[h^k(T)|T\ge t,W=w]=-[S(t|w){]}^{-1}{\int }_t^{\infty }[h(u){]}^{k}dS(u|w)$$

$$D(O_i;G,S)=A_{1i}(G)+B_{1i}(G,S)-C_{1i}(G,S).$$

$$L_2^{\ast }(\mathcal{O},\beta ;G,S)=n^{-1}\sum _{i=1}^n{\left(D(O_i;G,S)-\beta (W_i)\right)}^2.$$

$$L_{DR}^{(2)}(\mathcal{O},\beta ;G,S)=\frac{1}{n}\sum _{i=1}^n\left(A_{2i}+B_{2i}-C_{2i}-2(A_{1i}+B_{1i}-C_{1i})\beta (W_i)+(A_{0i}+B_{0i}-C_{0i})\beta (W_i{)}^2\right)$$

$$L_{DR}^{(2)}(\mathcal{O},\beta ;G,S)=\frac{1}{n}\sum _{i=1}^n\left((A_{2i}+B_{2i}-C_{2i})-2(A_{1i}+B_{1i}-C_{1i})\beta (W_i)+\beta (W_i{)}^2\right).$$

$$L_2^{\ast }(\mathcal{O},\beta ;G,S)=\frac{1}{n}\sum _{i=1}^n\left((A_{1i}+B_{1i}-C_{1i}{)}^2-2(A_{1i}+B_{1i}-C_{1i})\beta (W_i)+\beta (W_i{)}^2\right).$$

$$L_2^{\ast }(\mathcal{O},\beta ;G,S)+\eta ||\beta |{|}_p^2$$

$$L_{DR}^{(2)}(\mathcal{O},\beta ;G,S)+\eta ||\beta |{|}_p^2.$$

$$\sum _{l=1}^D\sum _{i=1}^n{A}_{i,l}{L}_2^{\ast }(\mathcal{O},f_{{\hat{\beta }}_{{\eta }_m}^{(l)}}(W);G,S)=\sum _{l=1}^D\sum _{i=1}^n{A}_{i,l}{L}_{DR}^{(2)}(\mathcal{O},f_{{\hat{\beta }}_{{\eta }_m}^{(l)}}(W);G,S)+K(\mathcal{O};G,S).$$

$$argmi{n}_{m\in \{1,\dots ,M\}}\sum _{l=1}^D\sum _{i=1}^n{A}_{i,l}{L}_2^{\ast }(\mathcal{O},f_{{\hat{\beta }}_{{\eta }_m}^{(l)}}(W);G,S)$$

$$=argmi{n}_{m\in \{1,\dots ,M\}}\sum _{l=1}^D\sum _{i=1}^n{A}_{i,l}{L}_{DR}^{(2)}(\mathcal{O},f_{{\hat{\beta }}_{{\eta }_m}^{(l)}}(W);G,S).$$

Full text

Deep Learning for Survival Outcomes

Jon Arni Steingrimsson

Department of Biostatistics, Brown University

Abstract: This manuscripts develops a new class of deep learning algorithms for outcomes that are potentially censored. To account for censoring, the unobservable loss function used in the absence of censoring is replaced by a censoring unbiased transformation. The resulting class of algorithms can be used to estimate both survival probabilities and restricted mean survival. We show how the deep learning algorithms can be implemented using software for uncensored data using a form of response transformation. Simulations and analysis of the Netherlands 70 Gene Signature Data show strong performance of the proposed algorithms.

Keywords: Machine Learning, Risk Estimation, Censoring Unbiased Transformations, $L_{2}$-loss, Semiparametric Theory, Doubly Robust Estimation

1 Introduction

Prediction models built using deep learning algorithms have had great success in many application areas including natural language processing (Goldberg, 2016), speech recognition (Graves et al., 2013), and image recognition (LeCun et al., 2015). Deep learning algorithms create a sequence of layers where each layer depends on an unknown vector of weights. The weight vector is estimated by minimizing a loss function often subject to some regularization.

In medical studies, the outcome of interest is commonly time to a specific event such as time until death or disease progression. Such outcomes are frequently subject to right-censoring, which occurs when a subject drops out from the study, dies from other causes, or the study ends before the participant experiences the event of interest. The main difficulty of adapting the deep learning algorithm to right censored outcomes is that the full data loss used in the absence of censoring cannot be calculated.

To overcome that challenge, Liao and Ahn (2016) and Ranganath et al. (2016) proposed a deep learning algorithm where the loss function assumes a Weibull distributed failure time. Building on previous work by Faraggi and Simon (1995), Katzman et al. (2018) proposed a deep learning algorithm to estimate the functional form of the covariates in an underlying proportional hazard model using a loss function based on the partial likelihood of a proportional hazard model. Mobadersany et al. (2018) used the algorithm to predict cancer survival based on digital pathology images and Li et al. (2019) used a similar algorithm to predict overall survival for rectal cancer patients. Finally, Luck et al. (2017) implemented a deep learning algorithm for potentially censored outcomes using a Cox model likelihood incorporating Efron’s method (Efron, 1977) to handle tied event times.

All of these algorithms share two common themes: i) they construct a prediction model where the survival time is assumed to have a parametric form or the loss function used relies on the proportional hazard assumption; and ii) they use a loss function that in the absence of censoring does not reduce to any commonly used deep learning algorithm for uncensored outcomes, creating a gap between methods used for censored and uncensored outcomes.

In the context of building a survival tree, Molinaro et al. (2004) developed an inverse probability censoring weighted (IPCW) loss that i) is an unbiased estimator for the full data risk that would be used when there is no censoring and ii) reduces to the corresponding full data loss in the absence of censoring. IPCW estimators are inefficient as they fail to utilize information in censored observations. Using semi-parametric efficiency theory for missing data, Steingrimsson et al. (2016) developed a “doubly robust” estimator for the full data risk that is both more efficient and more robust to the modeling choices made than the IPCW loss. Steingrimsson et al. (2017) proposed a class of censoring unbiased loss (CUL) functions that includes both the IPCW and the “doubly robust” losses as a special case.

This manuscript develops a class of deep learning algorithms for censored outcomes, referred to as censoring unbiased deep learning (CUDL), where the unobservable full data loss is replaced by the CUL functions. We show how the full data loss can be selected to estimate both survival probabilities and restricted mean survival. Furthermore, we show how the censoring unbiased deep learning algorithm can be implemented using software for fully observed continuous outcomes using a form of response transformation.

Section 2 reviews the deep learning algorithm when there is no censoring. Section 3 discusses loss estimation for time-to-event outcomes. Section 4 defines the censoring unbiased deep learning algorithm and shows how different choices of full data loss functions result in deep learning algorithms that estimate both survival probabilities and restricted mean survival. Section 5 shows how the CUDL algorithms with the full data loss as the $L_{2}$ loss can be implemented using software for fully observed data. Sections 6 and 7 discuss implementation of the CUDL algorithms and evaluate the performance of the deep learning algorithms using simulations and by analyzing the Netherlands Breast Cancer dataset, respectively. A Supplementary Web Appendix contains additional simulation results and proofs.

2 Deep Learning for Fully Observed Outcomes

In the absence of censoring, the dataset is assumed to consist of a positive continuous failure time $T\in\mathbb{R}^{+}$ and a covariate vector $W\in\mathbb{R}^{p}$. Assume that the outcome is possibly transformed using a monotone function $h:\mathbb{R}^{+}\rightarrow\mathbb{R}$ (e.g. $h(u)=u$ or $h(u)=log(u)$). With no censoring, the data is $\mathcal{F}=\{(W_{i},h(T_{i}));i=1,\ldots,n\}$. The dataset $\mathcal{F}$ is referred to as the fully observed data.

An important component of the deep learning algorithm is specification of a loss function (Goodfellow et al., 2016). A loss function $L(h(T),f(W))$ measures the discrepancy between a prediction $f(W)$ and an outcome $h(T)$. Examples of loss functions commonly used in connection with a continuous outcome are the $L_{2}$ loss $(h(T)-f(W))^{2}$ and the $L_{1}$ loss $|h(T)-f(W)|$.

In this manuscript we focus on the following deep learning algorithm based on feedforward networks. For a predetermined loss function, the full data deep learning algorithm is defined by the following steps:

Create a layer of hidden features. Fix the number of layers $K$ and for each layer $k\in\{1,\ldots,K\}$ pre-specify the function $f_{\beta_{k}}^{(k)}(x)$. The final function outputted by the hidden layer architecture is $f_{\beta}(W)=f_{\beta_{K}}^{(K)}\circ f_{\beta_{K-1}}^{(K-1)}\circ\ldots\circ f_{\beta_{1}}^{(1)}(W)$. Here, $\beta=(\beta_{1}^{T} ,\ldots,\beta_{K}^{T})^{T}$ is the vector of unknown weights that need to be estimated. 2. 2.

Estimate the weight vector. For a pre-specified loss function $L(h(T),\beta(W))$ and a fixed value of the scalar penalization parameter $\eta$, the vector of weights $\beta$ is estimated by minimizing the empirical penalized loss function

[TABLE]

Here, $||\beta||_{c}^{2}=\sum_{i=1}^{c}|\beta_{i}|^{2}$ and $c$ is the length of the weight vector. 3. 3.

Use cross-validation to select the penalization parameter $\eta$ from a pre-determined sequence $\eta_{1},\ldots,\eta_{M}$. Randomly split the dataset into $D$ disjoint sets $K_{1},\ldots K_{D}$. For fixed $l\in\{1,\ldots,D\}$ and $m\in\{1,\ldots,M\}$, define $\hat{\beta}_{\eta_{m}}^{(l)}(W)$ as the vector of weights estimated by minimizing (1) using the penalization parameter $\eta_{m}$ calculated using the data $\mathcal{F}/K_{l}$. Let $A_{i,l}$ be one if observation $i$ falls in dataset $K_{l}$ and zero otherwise. The cross-validation error corresponding to $\eta_{m}$ is defined as

[TABLE]

The final value of $\eta$ is $\eta_{m^{*}}$ where $m^{*}=argmin_{m\in\{1,\ldots,M\}}\alpha(m)$. 4. 4.

The final prediction model is $f_{\hat{\beta}_{\eta_{m^{*}}}}(W)$, where $\hat{\beta}_{\eta_{m^{*}}}$ is the value of $\beta$ that minimizes (1) with the penalization parameter set to $\eta_{m^{*}}$.

The population parameter that the full data deep learning algorithm estimates is the minimizer of the expected loss (risk) used in the algorithm. For the $L_{2}$ risk, the population parameter that the full data deep learning algorithm estimates is the conditional mean $E[h(T)|W]$.

3 Risk Estimation with Censored Data

In the presence of censoring, denoted $C$, the failure time is sometimes only partially observed. The data on observation $i$ is assumed to be $O_{i}=(\tilde{T}_{i}=\min(T_{i},C_{i}),\Delta_{i}=I(T_{i}\leq C_{i}),W_{i})$. Here, $I(\cdot)$ is the indicator function. We refer to $\tilde{T}$ is the observed time and $\Delta$ as the failure time indicator. We assume that the observed dataset $\mathcal{O}=\{O_{i};i=1,\ldots,n\}$ consists of $n$ independent and identically distributed observations.

Define $S_{0}(u|W)=P(T>u|W)$ and $G_{0}(u|W)=P(C>u|W)$ as the conditional survivor functions for $T$ and $C$, respectively. We assume that $C$ is continuous and that $C$ is independent of $T$ conditioned on $W$ (non-informative censoring). We also make the positivity assumption $G_{0}(\tilde{T}|W)\geq\epsilon>0$ for some $\varepsilon>0$.

The full data loss function $L(h(T),f(W))$ cannot be calculated if the failure time is censored. Replacing the unobservable full data loss with a loss function that can be calculated in the presence of censoring is the main difficulty in extending the deep learning algorithm to survival data.

An estimator is said to be an observed data estimator if it is a function of $\mathcal{O}$. Observed data estimators that are unbiased estimators of $\mathcal{R}(\beta)=E[L(h(T),\beta(W))]$ are referred to as censoring unbiased estimators for $\mathcal{R}(\beta)$. We now describe three censoring unbiased estimators for $\mathcal{R}(\beta)$ that all reduce to the full data loss when there is no censoring.

Inverse probability weighing (Robins and Rotnitzky, 1992) is a general missing data technique that performs a weighted complete case analysis where the weights are selected such that the weighted complete case estimator is an unbiased estimator for the desired full data target parameter. In the context of censored data risk estimation, the IPCW loss function is given by

[TABLE]

Using the law of iterated expectation, standard calculations show that the IPCW loss function is an unbiased estimator for the full data risk.

The censoring survival curve $G_{0}(u|W)$ is usually unknown and needs to estimated using some observed data estimator $\hat{G}(u|W)$. The empirical IPCW loss is given by

[TABLE]

The loss $L_{IPCW}(\mathcal{O},\beta;\hat{G})$ is a consistent estimator for $R(\beta)$ if the model for $C|W$ is correctly specified.

Censored observations do not contribute to the $L_{IPCW}(O,\beta;\hat{G})$ loss apart from potentially through the calculation of $\hat{G}(\cdot|\cdot)$. This leads to the IPCW loss being an inefficient estimator for the full data risk. With the aim of improving efficiency, Steingrimsson et al. (2016) used semi-parametric efficiency theory for missing data (Tsiatis, 2006; Robins et al., 1994) to develop an augmented estimator for $\mathcal{R}(\beta)$. This augmented loss function is given by

[TABLE]

where for a survival curve $S$

[TABLE]

Here, $\Lambda_{G}(u|W)=-\int_{0}^{u}dG(t|W)/G(u|W)$ is the cumulative hazard function.

The loss function (4) consists of the IPCW loss plus an augmentation term constructed to use information from censored responses. As shown in Steingrimsson et al. (2016), the loss (4) is the estimator of $R(\beta)$ with the smallest asymptotic variance among all unbiased estimators for $R(\beta)$ that can be written as the IPCW loss plus an augmentation term. This implies that (4) is an asymptotically more efficient estimator for $R(\beta)$ than $L_{IPCW}(O,\beta;G_{0})$.

Implementation of the augmented loss function (4) relies on estimating $G_{0}(\cdot|\cdot)$ and $S_{0}(\cdot|\cdot)$. Plugging observed data estimators $\hat{G}(\cdot|\cdot)$ and $\hat{S}(\cdot|\cdot)$ into (4) results in the empirical loss

[TABLE]

where

[TABLE]

The loss function $L_{DR}(\mathcal{O},\beta;\hat{G},\hat{S})$ is doubly robust in that it is a consistent estimator for $\mathcal{R}(\beta)$ if one of the models for $T|W$ or $C|W$ are correctly specified but not necessarily both. For that reason, $L_{DR}(\mathcal{O},\beta;\hat{G},\hat{S})$ is referred to as the empirical doubly robust loss.

A related class of simple and intuitive observed data estimators are the Buckley-James estimators (Buckley and James, 1979; Rubin and van der Laan, 2007). In the context of risk estimation, the Buckley-James estimator for $\mathcal{R}(\beta)$ is given by

[TABLE]

Implementation requires specifying an estimator for $S_{0}(\cdot|\cdot)$. Using a plug-in estimator for $S_{0}(\cdot|\cdot)$ results in the empirical Buckley-James loss

[TABLE]

The Buckley-James loss requires a consistent estimator for $S_{0}(\cdot|\cdot)$ in order to consistently estimate $\mathcal{R}(\beta)$. The Buckley-James estimator has the optimality property that $L_{BJ}(\mathcal{O},\beta;S_{0})$ is the function of the observed data that minimizes $E[(\mathcal{R}(\beta)-f(\mathcal{O}))^{2}]$ for any observed data function $f$ (Fan and Gijbels, 1994).

As $L_{BJ}(\mathcal{O},\beta;\hat{S})=L_{DR}(\mathcal{O},\beta;\hat{S},\hat{G}=1)$, the Buckley-James loss is equivalent to the doubly robust loss using the (incorrect) model specification $\hat{G}(t|w)=1$ for all $(t,w)$.

4 Censoring Unbiased Deep Learning

Replacing the full data loss function in the full data deep learning algorithm by any of the censoring unbiased loss functions results in a prediction model that can be calculated using the censored data $\mathcal{O}$. The algorithm obtained by replacing the full data loss by the $L_{IPCW}(O,\beta;\hat{G})$, $L_{DR}(O,\beta;\hat{G},\hat{S})$, and $L_{BJ}(O,\beta;\hat{S})$ are respectively referred to as the IPCW, doubly robust, and Buckley-James deep learning algorithms. Collectively we refer to these three deep learning algorithms as censoring unbiased deep learning (CUDL). When implemented using the full data loss as the $L_{2}$ loss we refer to the algorithms as $L_{2}$ censoring unbiased deep learning.

4.1 Algorithms Estimating Restricted Mean Survival

A popular target estimator for the full data deep learning algorithm for continuous outcomes is $E[T|W]$. Setting $h(u)=u$ and selecting the full data loss as the $L_{2}$ loss in the full data deep learning algorithm results in an estimator for $E[T|W]$. Estimating the mean with censored outcomes requires strong assumptions (ding2014estimating). A popular alternative for censored outcomes is to estimate the restricted mean survival $E[\min(T,\tau)|W]$ for some pre-specified constant $\tau$.

Selecting the full data loss as the $L_{2}$ loss and fitting the CUDL algorithms on the modified dataset

[TABLE]

results in an estimator for the restricted mean survival $E[\min(T,\tau)|W]$. We refer the CUDL algorithm estimating $E[\min(T,\tau)|W]$ as the restricted mean survival CUDL algorithm.

4.2 Algorithms Estimating $P(T\geq t|W)$

Survival probabilities of the form $P(T\geq t|W)$ for a fixed time-point $t$ are common target parameters for time-to-event data. The Brier risk $E[(I(T\geq t)-\beta(W))^{2}]$ induces $P(T\geq t|W)$ as the target parameter as it is the value of $\beta(W)$ that minimizes the Brier risk.

Applying the IPCW theory from Section 3 with the full data loss as the the Brier loss $(I(T\geq t)-\beta(W))^{2}$ gives

[TABLE]

The IPCW Brier loss fails to utilize that the value of $I(T\geq t)$ is known for censored observations with censoring times larger than $t$. Graf et al. (1999) used this information to create a censored data Brier loss and Lostritto et al. (2012) showed that the empirical censored data Brier loss has the IPCW representation

[TABLE]

Augmenting (10) results in the doubly robust Brier loss (Steingrimsson et al., 2017)

[TABLE]

where

[TABLE]

and $m_{L_{2},t}(u,W;S)=E_{S}[(I(\tilde{T}\geq t)-\beta(W))^{2}|T\geq u,W]$.

Following the developments in Section 3, the empirical Buckley-James Brier loss is given by

[TABLE]

Incorporating any of the ICPW, doubly-robust, or Buckley-James Brier loss functions into the CUDL algorithm results in a prediction model for $P(T\geq t|W)$. We refer to the CUDL algorithms with the full data loss selected as the Brier loss as the Brier CUDL algorithms.

5 Implementation of the $L_{2}$ Censoring Unbiased Deep Learning Algorithms

This section shows how a form of response transformation can be used to implement the doubly robust and Buckley-James CUDL algorithms described in Sections 4.1 and 4.2 using deep learning software implementing the full data deep learning algorithm with the $L_{2}$ loss.

Both the restricted mean survival and Brier CUDL algorithms use an $L_{2}$ full data loss of the form $(h(T)-\beta(W))^{2}$ with $h(T)=min(T,\tau)$ fit on the dataset $\mathcal{O}(\tau)$ and $h(T)=I(T\geq t)$ fit on the dataset $\mathcal{O}(t)$, respectively. Hence, both algorithms can be implemented using the form of response transformation now described.

For $k=0,1,2$, define

[TABLE]

and

[TABLE]

Here,

[TABLE]

for $k=1,2$ and $m_{0}(t,w;S)=1$ for all $(t,w)$. Define the response transformation

[TABLE]

Following Steingrimsson et al. (2017), $E[D(O;G_{0},S)]=E[D(O;G,S_{0})]=E[h(T)]$. Hence, $D(O;G,S)$ is a censoring unbiased transformation for $E[h(T)]$ if at least one of the models for $T|W$ and $C|W$ is correctly specified. Define the response transformed $L_{2}$ loss

[TABLE]

The response transformed $L_{2}$ loss is just the $L_{2}$ loss using the censoring unbiased outcome transformation $D(O;G,S)$ as the outcome. The doubly robust $L_{2}$ loss, denoted by $L^{(2)}_{DR}(\mathcal{O},\beta;\hat{G},\hat{S})$, is given by equation (6) using the full data $L_{2}$ loss $L(h(T),\beta(W))=(h(T)-\beta(W))^{2}$.

Steingrimsson et al. (2017) show that the loss functions $L^{*}_{2}(\mathcal{O},\beta;G,S)$ and $L^{(2)}_{DR}(\mathcal{O},\beta;G,S)$ are equivalent up to a term that is independent of $\beta$. An important consequence of that equivalence is Theorem 5.1. A proof of the theorem is presented in Supplementary Web Appendix S.2.

Theorem 5.1**.**

The prediction model created using the CUDL algorithm with the loss function $L^{(2)}_{DR}(\mathcal{O},\beta;G,S)$ is identical to the prediction model built using the fully observed deep learning algorithm implemented using the loss function $L^{*}_{2}(\mathcal{O},\beta;G,S)$.

Theorem 5.1 is general enough to allow for $G(\cdot|\cdot)=1$ and as $L_{DR}(\mathcal{O},\beta;G=1,S)=L_{BJ}(\mathcal{O},\beta;S)$ the result also holds for the Buckley-James loss. Requiring $m_{0}(t,w;S)=1$ for all $(t,w)$ is a key condition used to proof Theorem 5.1. The IPCW loss is a special case of the doubly robust loss with $m_{0}(t,w;S)=0$. As a consequence, the IPCW deep learning algorithm is not a special case of Theorem 5.1.

The main practical utility of Theorem 5.1 is that it allows the $L_{2}$ doubly robust and Buckley-James CUDL algorithms to be implemented using the following road-map:

Within the R software environment, the fully observed $L_{2}$ deep learning algorithm can be implemented using the Keras interface to R. Keras is a high level application programming interface that incorporates various types of backend engines such as Tensorflow to train deep learning models. There is a large literature on optimization techniques to estimate the weight vector for fully observed outcomes. The road-map described above allows easy implementation of the CUDL algorithms using optimization procedures available for fully observed outcomes and the $L_{2}$ loss. All the simulations and analysis presented in Section 6 and 7 are implemented using the response transformation with the Keras interface.

6 Simulations

6.1 Architecture and Implementation Choices

Implementation of the CUDL algorithms requires specifying the full data loss, the hidden layers, and the models for $S_{0}(\cdot|\cdot)$ and/or $G_{0}(\cdot|\cdot)$.

The CUDL algorithms used in the simulations have two layers (i.e. $K=2$). The first layer consists of a rectified linear unit activation function with $d_{1}=15$ output units. That is, $f_{\beta_{1}}^{(1)}:\mathbb{R}^{p}\rightarrow\mathbb{R}^{d_{1}}$ with $f_{\beta_{1}}^{(1)}(x)=max(x^{\prime}\beta_{1}^{*}+b_{1}^{*},0)$. Here, $\beta_{1}^{*}$ is a matrix of dimension $p\times d_{1}$ and the intercept $b_{1}^{*}$ is a vector of length $d_{1}$. The maximum in the above equation is an element-wise maximum. The vector of weights $\beta_{1}$ is given by $\beta_{1}=((\beta_{1}^{(*1)})^{T},\ldots,(\beta_{1}^{(*d_{1})})^{T},(b_{1}^{*})^{T})^{T}$, where $\beta_{1}^{(*k)}$ is the k-th column of $\beta_{1}^{*}$.

For the Brier CUDL, the second layer $f_{\beta_{2}}^{(2)}:\mathbb{R}^{d_{1}}\rightarrow[0,1]$ uses the sigmoid activation function $f_{\beta_{2}}^{(2)}(x)=(1+e^{-x^{\prime}\beta_{2}^{*}+b_{2}^{*}})^{-1}$, where $\beta_{2}^{*}$ is a vector of length $d_{1}$ and $b_{2}^{*}$ is a scalar. The second vector of weights is $\beta_{2}=((\beta_{2}^{*})^{T},b_{2}^{*})^{T}$. The sigmoid activation function ensures that the final prediction falls in the interval $[0,1]$, respecting the natural boundary of the target parameter $P(T\geq t|W)$.

The restricted mean survival CUDL algorithms use the same architecture apart from the second layer being a rectified linear unit activation function instead of the sigmoid activation function (i.e. $f_{\beta_{2}}^{(2)}(x)=max(x^{\prime}\beta_{2}^{*}+b_{2}^{*},0)$). This reflects that the target estimator is no longer a probability. For both CUDL algorithms, the dimension of the weight vector $\beta=(\beta_{1}^{T},\beta_{2}^{T})^{T}$ is $p*d_{1}+2*d_{1}+1$.

The censoring survival curve $G_{0}(t|w)$ required for implementation of the CUDL algorithms is estimated using the survival tree method of LeBlanc and Crowley (1992). It is implemented using the rpart package in R. Default tuning parameters in rpart are used apart from that the minimum number of observations in a terminal node is set to $30$. To ensure that the positivity assumption holds empirically, a “Method 2” truncation as described in Steingrimsson et al. (2016) is used within each terminal node. “Method 2” truncation within a terminal node modifies the data by setting the failure time indicator for the largest $10\%$ of observations falling in each terminal node to one. To estimate $S_{0}(t|W)$ we use the random survival forest procedure (Ishwaran et al., 2008) with default tuning parameters.

Five fold cross-validation is used to select a final penalization parameter from the sequence $(0,0.001,0.01,0.1)$. All covariates are standardized prior to fitting the CUDL algorithms. To solve the minimization problem required to estimate the weight vector we use the rmsprop minimization procedure in the keras interface. We use default values for the tuning parameters except for setting the epochs parameter to 100 and we use $20\%$ dropout (Dahl et al., 2013). R code implementing the CUDL algorithm analysis presented in Section 7 is publicly available from github.com/jonsteingrimsson/CensoringDL.

6.2 Simulation Setup and Results

We use simulations to evaluate the performance of the Buckley-James and doubly robust CUDL algorithms. We do not include the IPCW CUDL algorithm as i) the IPCW loss is both less efficient and less robust to model misspecifications than the doubly robust loss and ii) Theorem 5.1 does not hold for the IPCW CUDL algorithm, preventing implementation using the response transformation described in Section 5.

We compare the two CUDL algorithms to three other prediction models for censored data, a main effects Cox model, a penalized Cox model, and random survival forests. The main effects Cox model is implemented using the coxph function in the survival package in R. The penalized Cox model is implemented using the cv.glmnet function in R. All default values of the tuning parameters were used for the penalized Cox model. This includes selecting the penalization parameter using $10$ fold cross-validation. The random survival forest algorithm (Ishwaran et al., 2008) is implemented using the rfsrc function from the randomForestSRC package (Ishwaran and Kogalur, 2007) with all tuning parameters set to their default values.

We use the following simulation settings to compare the performance of the five algorithms.

1. Setting 1. The covariate vector is simulated from a $30$ dimensional multivariate normal distribution with mean zero and covariance matrix with element $(i,j)$ equal to $0.5^{|i-j|}$. The failure time distribution is exponential with mean $e^{0.1\sum_{j=1}^{10}W^{(j)}}$, where $W^{(j)}$ is the $j$-th component of $W$. The censoring distribution is exponential with mean $1.14$, which results in approximately $47\%$ censoring. The training set consists of $1000$ i.i.d. draws from the joint distribution of $(\tilde{T},\Delta,W)$. In this setting, the proportional hazard assumption is satisfied.
2. Setting 2. The covariate vector is simulated from a $30$ dimensional multivariate normal distribution with mean zero and covariance matrix with element $(i,j)$ equal to $0.5^{|i-j|}$. The failure time is simulated from a gamma distribution with shape parameter $0.5+0.3|\sum_{j=11}^{15}W^{(j)}|$ and scale parameter equal to $2$. The censoring times are uniformly distributed on the interval $[0,15]$, which results in approximately $18\%$ censoring. The training set consists of $1000$ i.i.d. draws from the joint distribution of $(\tilde{T},\Delta,W)$. In this setting, the proportional hazard assumption is violated.

We compare the algorithms both when predicting $P(T\geq t|W)$ for a pre-specified time-point $t$ and when predicting restricted mean survival $E[\min(T,\tau)|W]$ for a pre-specified time-point $\tau$. For both simulation settings, the timepoint $t$ is selected as the median of the marginal failure time distribution and $\tau$ is set to the $85$th quantile of the marginal observed time distribution.

Both Cox models and the random survival forest algorithm estimate the survival curve $S_{0}(t|w)$. To estimate the restricted mean survival for these algorithms we calculate $\hat{S}(t|w)$ and use the formula $\hat{E}[\min(T,\tau)|W=w]=-\int_{0}^{\infty}\min(t,\tau)d\hat{S}(t|w)$.

The doubly robust and Buckley-James CUDL algorithms are implemented using the doubly robust and Buckley-James loss functions. To estimate $P(T\geq t|W)$ the CUDL algorithms use the full data loss as the Brier loss and to estimate $E[\min(T,\tau)|W]$ the CUDL algorithms use the estimation procedure detailed in Section 4.1. For both simulation settings the length of the weight vector $\beta$ is $481$.

When estimating the survival probability $P(T\geq t|W)$, each of the algorithms is fit on the training set and the resulting model fit is used to predict $P(T\geq t|W)$ on an independent test set consisting of $1000$ observations simulated using the corresponding full data distribution. The final evaluation measure is the average $L_{2}$ distance between the predicted test set probability and the true probability $\frac{1}{1000}\sum_{i=1}^{1000}(P(T\geq t|W_{i})-\hat{P}(T\geq t|W_{i}))^{2}$. When the target parameter is the restricted mean survival, the probability $P(T\geq t|W)$ is replaced by $E[\min(T,\tau)|W]$ in the description above.

The results from a $1000$ simulations are shown in Figure 1. We see that both CUDL algorithms perform substantially better than the random survival forest algorithm for both settings and both target parameters. When compared to the Cox models, the CUDL algorithms perform substantially better when the proportional hazard assumption is violated (setting 2). For the setting where the proportional hazard assumption holds (setting 1), the CUDL algorithms show similar performance to the Cox model when estimating survival probabilities but perform slightly worse than the Cox models when estimating restricted mean survival. The doubly robust and Buckley-James CUDL algorithms show similar performance for both settings and target parameters.

Supplementary Web Appendix S.1 presents additional simulation results for estimating $P(T\geq t|W)$ using modifications of settings one and two.

• •

Figure S-1 shows comparisons for the five algorithms when the sample size is $250$, $500$, $1500$, and $3000$. The performance of all algorithms improves as sample size increases and the relative performance of the algorithms is similar to what is seen in Figure 1. In setting 2, the improvements of the CUDL algorithms compared to the Cox model becomes larger as the sample size is increased.

• •

Figure S-2 shows simulations results when the covariate dimension is increased to $100$. In the setting where the proportional hazard assumption holds, the penalized Cox model shows the best performance followed by the Buckley-James CUDL algorithm. When the proportional hazard assumption is violated, the CUDL algorithms outperform the other methods.

• •

Figures S-3 and S-4 show simulations results when the time-point $t$ is set to the $25$th and $75$th quantile of the marginal failure time distribution, respectively. For the $75$th quantile and setting 2, the relative performance of the random survival forest algorithm is better than for the $50$th quantile and the performance is comparable to both CUDL algorithms. For all other combinations of quantiles and simulation settings, the relative performance of the algorithms is similar to what is seen in Figure 1

7 Comparing Prediction Accuracy using the Netherlands 70 Gene Signature Data

The Netherlands Cancer Institute 70 gene signature dataset consists of data from $144$ lymph node positive breast cancer patients. The dataset includes five risk factors (diameter of tumor, number of positive nodes, ER status, grade, and age) and 70 measures of gene expression (Van’t Veer et al., 2002). We use the data to evaluate the prediction accuracy of the CUDL algorithms when predicting the probability of metastasis-free survival beyond a specific time-point (measured in months). Patients who were alive at the end of study, developed a second primary cancer, had recurrence of regional or local disease, or died from other causes than breast cancer are considered censored. The censoring rate is $67\%$. The dataset is publicly available from the R package penalized.

We compare the prediction accuracy of the doubly robust and Buckley-James CUDL algorithms to a penalized Cox model and the random survival forest algorithm. The main effects Cox model is not included as the algorithm failed to converge and had low prediction accuracy. All implementation choices needed to fully define the four algorithms are as described in Section 6.2, this includes using a survival tree to estimate $G_{0}(u|w)$ and the random survival forest algorithm to estimate $S_{0}(u|w)$.

The four algorithms estimate $P(T\geq t|W)$ for a sequence of fixed time-points $t$ and we compare the prediction accuracy using the censored data Brier score given by equation (10). To calculate the censored data Brier score we use five fold cross-validation where the cross-validation is done such that all the cross-validation sets have approximately equal censoring rate.

Figure 2 shows the median of the censored data Brier score as a function of $t$ across $100$ different splits into cross-validation sets. For all time-points $t$, the CUDL algorithms have lower or similar prediction accuracy compared to the penalized Cox model and the random survival forest algorithm. The doubly robust and Buckley-James algorithms show similar prediction accuracy.

8 Discussion

This manuscript developed a class of deep learning algorithms for time-to-event outcomes by replacing the unobservable full data loss used in the absence of censoring by a censoring unbiased loss function. We show how the full data loss can be selected to estimate both survival probabilities and restricted mean survival. Furthermore, we show that the doubly robust and Buckley-James deep learning algorithms can be implemented using standard software for fully observed outcomes using a form of response transformation. The performance of the algorithms is evaluated both using simulations and by analyzing a dataset on breast cancer patients.

The Brier CUDL algorithms estimate $P(T\geq t|W)$ for a fixed time-point $t$, while many algorithms such as the Cox model and the random survival forest algorithm estimate the whole survival curve $S_{0}(t|w)$. The CUDL algorithm can be used to calculate an estimator for the whole survival curve by by using the CUDL algorithm to calculate $P(T\geq t|W)$ setting $t$ to each unique failure time in the dataset and assuming that the survival curve only jumps at observed failure times. However, this procedure becomes computationally intensive for large sample sizes.

There are several interesting future research directions arising from this work. Extensions to more complex data structures such as competing risk and time-dependent covariates are of interest. Furthermore, appropriately handling missing data both when the missingness mechanism is unknown and known (e.g. case-cohort studies) is of importance.

Implementation of the doubly-robust and Buckley-James algorithm requires an estimator for $S_{0}(\cdot|\cdot)$. It would be interesting to utilize an iterative algorithm to update the estimator for $S_{0}(\cdot|\cdot)$ using the CUDL algorithm. More precisely, at the first iteration the random survival forest algorithm is used as the estimator for $S_{0}(\cdot|\cdot)$ in the CUDL algorithm. The resulting CUDL algorithm is then used to estimate $S_{0}(\cdot|\cdot)$ and an updated CUDL algorithm is calculated with the updated estimator $S_{0}(\cdot|\cdot)$. This process is repeated with an updated estimator for $S_{0}(\cdot|\cdot)$ until convergence or for a fixed amount of iterations.

Acknowledgment

The author thanks Constantine Gatsonis and and Samantha Morrison for helpful comments on an earlier draft.

Supplementary Web Appendix

References to figures, tables, theorems and equations preceded by “S-” are internal to this supplement; all other references refer to the main paper.

S.1 Additional Simulation Results

S.1.1 Impact of Sample Size on Performance

Figure S-1 shows the prediction accuracy of the five different algorithms for sample sizes of $250,500,1500$, and $3000$. The simulation settings used are the same as described in Section 6.2 apart from the sample sizes are changed.

The CUDL algorithms show the overall best performance for all sample sizes. For setting 2, the CUDL algorithms show superior performance for all sample sizes. For setting 1, the CUDL algorithms show comparable or slightly worse performance to the Cox models and superior performance compared to the random survival forest algorithm.

S.1.2 Increasing Covariate Dimension

Figure S-2 shows simulation results when the covariate dimension is increased to $100$ for the simulation settings described in Section 6.2. For both settings, the additional covariates are noise variables not affecting either the failure time or censoring distribution. The covariate vector is simulated from a $100$ dimensional multivariate normal distribution with mean zero and covariance matrix with element $(i,j)$ equal to $0.5^{|i-j|}$. The results show similar trends as seen in the simulations presented in Figure 1.

S.1.3 Simulations for other time-points

Figures S-3 and S-4 compare performance of the five prediction models in both settings used in Section 6.2 when estimating $P(T\geq t|W)$ with $t$ selected as the $25$th and $75$th quantile of the marginal failure time distribution, respectively. The results show similar trends as seen in Figure 1, except that for the $75$th quantile we see better relative performance of the Cox model in Setting 1 and better relative performance of the random survival forest algorithm in Setting 2.

S.2 Proof of Theorem 5.1

Proof of Theorem 5.1: As $L_{BJ}(\mathcal{O},\beta;S)=L_{DR}(\mathcal{O},\beta;1,S)$ it is enough to show the stated equivalence for the $L_{DR}(\mathcal{O},\beta;G,S)$ loss under regularity conditions that are flexible enough to allow $G(t|w)=1$ for all $(t,w)$.

Using the notation defined in Section 5

[TABLE]

Lemma $1$ in Strawderman (2000) gives that $A_{0i}+B_{0i}-C_{0i}=1$ for all $i=1,\ldots,n$. Hence,

[TABLE]

Expanding the square shows that the response transformed $L_{2}$ loss can be written as

[TABLE]

The above two formulas show that $L^{*}_{2}(\mathcal{O},\beta;G,S)$ and $L_{DR,2}(\mathcal{O},\beta;G,S)$ are equivalent up to a term that is independent of $\beta(W)$. An important consequence is that for a fixed penalization parameter the weight vector which minimizes

[TABLE]

is the same as the weight vector which minimizes

[TABLE]

Hence, it is enough to show that the penalization parameter selected by minimizing the cross-validated loss with the loss $L^{*}_{2}(\mathcal{O},\beta;G,S)$ is the same as the penalization parameter selected by minimizing the cross validated loss with the loss function $L^{(2)}_{DR}(\mathcal{O},\beta;G,S)$. Previous calculations show that

[TABLE]

where $K(\mathcal{O};G,S)$ does not depend on $\beta$. This shows that

[TABLE]

Hence, the cross-validation procedure for both loss functions results in the same final penalization parameter.