Deep Recurrent Survival Analysis

The event occurrence information of some samples may be lost, due to some limitation of the observation period or losing tracks during the study procedure [\citeauthoryearWang, Li, and Reddy2017], which is called data censorship. When dealing with time-to-event information, a more complex learning problem is to estimate the probability of the event occurrence at each time, especially for those samples without tracking logs after (or before) the observation time which is defined as right-censored (or left-censored) [\citeauthoryearWang, Li, and Reddy2017]. Survival analysis is a typical statistical methodology for modeling time-to-event data while handling censorship. There are two main streams of survival analysis.

The first view is based on traditional statistics scattering in three categories. (i) Non-parametric methods including Kaplan-Meier estimator [\citeauthoryearKaplan and Meier1958] and Nelson-Aalen estimator [\citeauthoryearAndersen et al.2012] are solely based on counting statistics, which is too coarse-grained to perform personalized modeling. (ii) Semi-parametric methods such as Cox proportional hazard model [\citeauthoryearCox1992] and its variants Lasso-Cox [\citeauthoryearTibshirani1997] assumes some base distribution functions with the scaling coefficients for fine-tuning the final survival rate prediction. (iii) Parametric models assume that the survival time or its logarithm result follows a particular theoretical distribution such as Exponential distribution [\citeauthoryearLee and Wang2003] and Weibull distribution [\citeauthoryearRanganath et al.2016]. These methods either base on statistical counting information or pre-assume distributional forms for the survival rate function, which generalizes not very well in real-world situations.

The second school of survival analysis takes from machine learning perspective. Survival random forest which was first proposed in [\citeauthoryearGordon and Olshen1985] derives from standard decision tree by modeling the censored data [\citeauthoryearWang et al.2016] while its idea is mainly based on counting-based statistics. Other machine learning methodologies include Bayesian models [\citeauthoryearRanganath et al.2015, \citeauthoryearRanganath et al.2016], support vector machine [\citeauthoryearKhan and Zubek2008] and multi-task learning solutions [\citeauthoryearLi et al.2016, \citeauthoryearAlaa and van der Schaar2017]. Note that, deep learning models have emerged in recent years. ? (?) first embedded neural network into Cox model to improve covariate relationship modeling. From that, many works applied deep neural networks into well-studied statistical models to improve feature extraction and survival analysis through end-to-end learning, such as [\citeauthoryearRanganath et al.2016, \citeauthoryearLuck et al.2017, \citeauthoryearLao et al.2017, \citeauthoryearKatzman et al.2018, \citeauthoryearGrob et al.2018]. Almost all the above models assume particular distribution forms which also suffers from the generalization problem in practice.

? (?) utilized neural network directly to predict the survival rate for each sample and ? (?) extended it to a Bayesian network method. In [\citeauthoryearLee et al.2018] the authors proposed a feed forward deep model to directly predict the probability density values at each time point and sum them for estimating the survival rate. However, in that paper, the gradient signal is quite sparse for the prediction outputs from the neural network. Moreover, to our best knowledge, none of the related literatures considers the sequential patterns within the feature space over time. We propose a recurrent neural network model predicting the conditional probability of event at each time and estimate the survival rate through the probability chain rule, which captures the sequential dependency patterns between neighboring time slices and back-propagate the gradient more efficiently.

Due to its adequate model capability and the support of big data, deep learning, a.k.a. deep neural network, has drawn great attention ranging from computer vision [\citeauthoryearKrizhevsky, Sutskever, and Hinton2012] and speech recognition [\citeauthoryearGraves, Mohamed, and Hinton2013] to natural language processing [\citeauthoryearBahdanau, Cho, and Bengio2015] during the recent decades. Among them, recurrent neural network (RNN) whose idea firstly emerged two decades ago and its variants like long short-term memory (LSTM) [\citeauthoryearHochreiter and Schmidhuber1997] employ memory structures to model the conditional probability which captures dynamic sequential patterns. In this paper we borrow the idea of RNN and well design the modeling methodology for survival function regression.

We define $z$ as the variable of the true occurrence time for the event of interest if it has been tracked. We just simplify the occurrence of the event of interest as event, and define the probabilistic density function (P.D.F.) of the true event time as $p(z)$, which means the probability that the event truely occurs at time $z$.

Now that we have the P.D.F. of the event time, we can derive the survival rate at each time $t$ as the C.D.F. as

which is the probability of the observing object surviving, i.e., event not occurring, until the observed time $t$. Then the straightforward definition of the event rate, i.e., the probability of event occurring before the observing time $t$, as that

The data of the survival analysis logs are represented as a set of $N$ triples $\{\boldsymbol{x},z,t\}_{i=1}^{N}$, where $t>0$ is the observation time for the given sample. Here $z$ is left unknown (and we marked $z$ as null) for the censored samples without the observation of the true event time. $\boldsymbol{x}$ is the feature vector of the observation which encodes different information under various scenarios.

Our goal is to model the distribution of the true event time $p(z)$ over all the historical time-to-event logs with handling the censored data of which the true event time is unknown. So the main problem of survival analysis is to estimate the probability distribution $p(z|\boldsymbol{x})$ of the event time with regard to the sample feature $\boldsymbol{x}$, for each sample. Formally speaking, the derived model is a “mapping” function $T$ which learns the patterns within the data and predicts the event time distribution over the time space as

First of all, we present the definition of the conditional hazard rate over continuous time as

which models the instant occurrence probability of the event at time $t$ given that the event has not occurred before. Note that the concept of hazard rate has been commonly utilized in many survival analysis literatures [\citeauthoryearCox1992, \citeauthoryearFaraggi and Simon1995, \citeauthoryearLuck et al.2017].

In the discrete context, a set of $L$ time slices $0<t_{1}<t_{2}<\ldots<t_{L}$ is obtained which arises from the finite precision of time determinations. Analogously we may also consider the grouping of continuous time as $l=1,2,\ldots,L$ and uniformly divide disjoint intervals $V_{l}=(t_{l-1},t_{l}]$ where $t_{0}=0$ and $t_{l}$ is the last observation interval boundary for the given sample, i.e., the tracked observation time in the logs. $V_{L}$ is the last time interval in the whole data space. This setting is appropriately suited in our task and has been widely used in clinical research [\citeauthoryearLi et al.2016, \citeauthoryearLee et al.2018], information systems [\citeauthoryearJing and Smola2017, \citeauthoryearGrob et al.2018] and other related fields [\citeauthoryearWu, Yeh, and Chen2015, \citeauthoryearWang et al.2016].

As such, our event rate function and survival rate function over discrete time space is

where the input to the two functions is the observed time $t_{l}$ from the log. And the discrete event time probability function at the $l$-th time interval is

The discrete conditional hazard rate $h_{l}$, defined as the conditional probability as

which approximates the continuous conditional hazard rate function $h(t_{l})$ in Eq. (4) as the intervals $V_{l}$ become infinitesimal.

Till now, we have presented the discrete time model and discuss the death (i.e., event) and survival probability over the discrete time space. We here propose our DRSA model based on recurrent neural network with the parameter $\boldsymbol{\theta}$, which captures the sequential patterns for conditional probability $h^{i}_{l}$ at every time interval $V_{l}$ for the $i^{\text{th}}$ sample.

The detailed structure of DRSA network is illustrated in Figure 1. At each time interval $V_{l}$, the $l$-th RNN cell predicts the instant hazard rate $h^{i}_{l}$ given the sample feature $\boldsymbol{x}^{i}$ and the current time $t_{l}$ conditioned upon the previous events as

where $f_{\boldsymbol{\theta}}$ is the RNN function taking $(\boldsymbol{x}^{i},t_{l})$ as input and $h^{i}_{l}$ as output. $\boldsymbol{r}_{l-1}$ is the hidden vector calculated from the RNN cell at the last time step which contains the information about the conditional. It is quite natural for using the recurrent cell to model the conditional probability over time [\citeauthoryearBahdanau, Cho, and Bengio2015]. In our paper we implement the RNN function as a standard LSTM unit [\citeauthoryearHochreiter and Schmidhuber1997], which has been widely used in sequence data modeling. The details of the implementation of our RNN architecture can be referred in our supplemental materials and our reproductive code published in the experiment part.

From Eqs. (5), (7) and (8), we can easily derive the survival rate function $S(t)$ through probability chain rule, and the corresponding event rate function $W(t)$ at the time $t$ for the $i^{\text{th}}$ individual sample as

where $l^{i}$ is the time interval index for the $i^{\text{th}}$ sample at $t^{i}$.

Moreover, taking Eqs. (6) and (7) into consideration, the probability of time $z$ lying in the interval of $V_{l^{i}}$ for the $i^{\text{th}}$ sample is

By means of probability chain rule, it connects all the outputs of the conditional hazard rate $h$ at each individual time to the final prediction, i.e., the probability $p(z)$ of the true event time $z$ and the survival rate $S(t)$ at each time $t$. This feed-forward calculation guarantees that the gradient signal from the loss function can be transmitted through back-propagation more effectively comparing with [\citeauthoryearLee et al.2018], which will be discussed below.

Since there is no ground truth of either the event time distribution or survival rate, here we maximize the log-likelihood over the empirical data distribution to learn our deep model. Specifically, we take three objectives as our losses.

The first loss is to minimize the negative log-likelihood of the true event time $z=z^{i}$ over the uncensored logs as

where $l^{i}$ is the index of the interval of the true event time $z^{i}\in V_{l^{i}}$.

The second loss is to minimize the negative partial log-likelihood of the event rate over the uncensored logs as

This loss adds more supervisions onto the predictions over the time range $(z^{i},t^{i})$ for the uncensored data than those [\citeauthoryearKatzman et al.2018, \citeauthoryearLee et al.2018, \citeauthoryearTibshirani1997] merely supervise on the true event time $z^{i}$.

Though the censored logs do not contain any information about the true event time, we would only know that the true event time $z$ is greater than our logged observing time $t^{i}$ then. Here we incorporate the partial log-likelihood embedded in the censored logs as the third loss to correct the learning bias of our model as

In this section, we unscramble some intrinsic properties of our deep survival model.

First of all, we analyze the model effectiveness of DRSA. In [\citeauthoryearLee et al.2018], the proposed deep model directly predicts the event probability $p(t)$ and combines to estimate the survival rate as $S(t)=\sum_{t^{\prime}\leq t}p(t^{\prime})$ while ignoring the sequential patterns. As a result, the gradient signal would only have effect on the prediction at time $t$ individually. On the contrary, from Figure 1, we can see that our DRSA model is obviously more effective since the supervision would be directly back-propagated through the chain rule calculation to all the units with strict mathematical derivation, which guarantees to transmit the gradient more efficiently and effectively. We also explicitly model the sequential patterns by conditional hazard rate prediction and we will illustrate the advantage of that in the experiments.

Then we take the view of censorship prediction of our methodology. As is known that there is a censoring status as an indicator of survival at the given time, for each sample as

In the tracking logs, each sample $(\boldsymbol{x}^{i},z^{i},t^{i})$ is uncensored where $c^{i}=0$. While for the censored logs losing tracking at the observation time, the true event time $z^{i}$ is unknown but the tracker only has the idea that $z^{i}>t^{i}$, thus $c^{i}=1$.

Moreover, for the uncensored data, it is natural to “push down” the probability of survival $S(t^{i})$. And for the censored data, it needs to “pull up” $S(t^{i})$ since we “observe event not occurred” at time $t^{i}$. However, using only $L_{z}$ to supervise the prediction of $p(z^{i})$ at time $z^{i}$ in Eq. (12) is insufficient. So that we incorporate the two partial likelihood losses $L_{\text{uncensored}}$ and $L_{\text{censored}}$ in Eqs. (13) and (14).

Therefore, taking Eqs. (13) and (14) altogether, we may find that the combination of $L_{\text{uncensored}}$ and $L_{\text{censored}}$ describes the classification of survival status at time $t^{i}$ of each sample as

which is the cross entropy loss for predicting the survival status at time $t^{i}$ given $\boldsymbol{x}^{i}$ over all the data $\mathbb{D}_{\text{full}}=\mathbb{D}_{\text{uncensored}}\bigcup\mathbb{D}_{\text{censored}}$.

Combining all the objective functions and our goal is to minimize the negative log-likelihood over all the data samples including both uncensored and censored data as

where $\boldsymbol{\theta}$ is the model parameter in Eq. (8) and the hyperparameter $\alpha$ controls the loss value balance between them. Specifically, $\alpha$ controls the magnitudes of the two losses at the same level to stabilize the model training.

We also analyze the model efficiency in the supplemental material of this paper and the time complexity of model inference is the same as the traditional RNN model which has proven practical efficiency in the industrial applications [\citeauthoryearZhang et al.2014].

We evaluate all the compared models in three real-world tasks. We also published the processed full datasets²² 2 We have put sampled data in the published code. The three processed full datasets link: https://goo.gl/nUFND4..

CLINIC

is a dataset for tracking the patient clinic status [\citeauthoryearKnaus et al.1995]. Here the goal of survival analysis is to estimate the time till the event (death), and predict the probability of the event with waning effects of baseline physiologic variables over time.

MUSIC

is a user lifetime analysis dataset [\citeauthoryearJing and Smola2017] that contains roughly 1,000 users with entire listening history from 2004 to 2009 on last.fm, a famous online music service. Here the event is the user visit to the music service and the goal is to predict the time elapsed from the last visit of one user to her next visit.

BIDDING

is a real-time bidding dataset in the computational advertising field [\citeauthoryearRen et al.2018, \citeauthoryearWang et al.2016]. In this scenario, the time is correspondent to the bid price of the bidder and the event is just winning of the auction. The feature contains the auction request information. Many researchers [\citeauthoryearWu, Yeh, and Chen2015, \citeauthoryearWang et al.2016] utilized survival analysis for unbiased winning probability estimation of a single auction while handling the losing (censored) logs without knowing the true winning price.

The statistics of the three datasets are provided in Table 1. We split the CLINIC and MUSIC datasets to training and test sets with ratio of 4:1 and 6:1, respectively. For feature engineering, all the datasets have been one-hot encoded for both categorical and numerical features. The original BIDDING data have already been feature engineered and processed as training and test datasets. Note that, the true time of the event of all the testing data have been preserved for the performance evaluation. In these datasets, since all the time is integer value, we bucketize the discrete time interval as interval size $s_{\text{intv}}=1$ and the maximal time interval number $L$ is equal to the largest integer time in each dataset. The discussion about various interval sizes has been included in the supplemental materials.

The first evaluation metric is the time-dependent concordance index (C-index), which is the most common evaluation used in survival analysis [\citeauthoryearLi et al.2016, \citeauthoryearLuck et al.2017, \citeauthoryearLee et al.2018] and reflects a measure of how well a model predicts the ordering of sample event times. That is, given the observing time $t$, two samples $d^{1}=(\boldsymbol{x}^{1},z^{1})$ with large event time $z^{1}>t$ and $d^{2}=(\boldsymbol{x}^{2},z^{2})$ with small event time $z^{2}\leq t$ should be ordered as $d^{1}\prec d^{2}$ where $d^{1}$ is placed before $d^{2}$. This evaluation is quite similar to the area under ROC curve (AUC) metric in the binary classification tasks [\citeauthoryearWang, Li, and Reddy2017]. From the ranking view of event probability estimation at time $t$, C-index assesses the ordering performance among all the uncensored and censored pairs at $t$ among the test data.

We also use average negative log probability (ANLP) to evaluate the regression performance among different forecasting models. ANLP is to assess the likelihood of the co-occurrence of the test sample with the corresponding true event time, which is correspondent to the event time likelihood loss $L_{z}$ in Eq. (12). Here we compute ANLP as

where $p(z|\boldsymbol{x})$ is the learned time-to-event probability function of each model.

Finally, we conduct the significance test to verify the statistical significance of the performance improvement of our model w.r.t. the baseline models. Specifically, we deploy a MannWhitney U test [\citeauthoryearMason and Graham2002] under C-index metric, and a t-test [\citeauthoryearBhattacharya and Habtzghi2002] under ANLP metric.

We compare our model with two traditional statistic methods and five machine learning methods including two deep learning models.

• •

  KM is Kaplan-Meier estimator [\citeauthoryearKaplan and Meier1958] which is a statistic-based non-parametric method counting on the event probability at each time over the whole set of samples.

• •

  Lasso-Cox is a semi-parametric method [\citeauthoryearZhang and Lu2007] based on Cox proportional hazard model [\citeauthoryearCox1992] with $l1$-regularization.

• •

  Gamma is a parametric gamma distribution-based regression model [\citeauthoryearZhu et al.2017a]. The event time of each sample is modeled by a unique gamma distribution with respect to its features.

• •

  STM is a survival random tree [\citeauthoryearWang et al.2016] model which splits the data into small segments using between-node heterogeneity and utilizes Kaplan-Meier estimator to derive the survival analysis for each segment.

• •

  MTLSA is the recently proposed multi-task learning with survival analysis model [\citeauthoryearLi et al.2016]. It transforms the original survival analysis problem into a series of binary classification subproblems, and uses a multi-task learning method to model the event probability at different time.

• •

  DeepSurv is a Cox proportional hazard model with deep neural network [\citeauthoryearKatzman et al.2018] for feature extraction upon the sample covariates.

• •

  DeepHit is a deep neural network model [\citeauthoryearLee et al.2018] which predicts the probability $p(z)$ of event over the whole time space with the input $\boldsymbol{x}$. This method achieved state-of-the-art performance in survival analysis problem.

• •

  DRSA is our proposed model which has been described above. The implementation details can be referred to supplemental materials and the published code.

We present the evaluation results according to the category which the compared models belong to. KM and Lasso-Cox model are statistic-based methods, while Gamma, STM and MTLSA are machine learning based models. The rest models are deep neural network models with end-to-end learning paradigm.