Year: 2019

Abstract:
The goal of lifetime clustering is to develop an inductive model that maps subjects into K clusters according to their underlying (unobserved) lifetime distribution. We introduce a neural-network based lifetime clustering model that can find cluster assignments by directly maximizing the divergence between the empirical lifetime distributions of the clusters. Accordingly, we define a novel clustering loss function over the lifetime distributions (of entire clusters) based on a tight upper bound of the two-sample Kuiper test p-value. The resultant model is robust to the modeling issues associated with the unobservability of termination signals, and does not assume proportional hazards. Our results in real and synthetic datasets show significantly better lifetime clusters (as evaluated by C-index, Brier Score, Logrank score and adjusted Rand index) as compared to competing approaches.

heading: 1 INTRODUCTION
text: Survival analysis is widely used to model the relationship between subject covariates and the time until a particular terminal event of interest (e.g., death, or quitting of social media) that marks the end of all activities (or measurements) of that subject (known as the lifetime of the subject). For instance, a subject’s logins to a social network are her activities and the time until she quits the social network permanently is her lifetime.
The lifetime of a subject can be unobserved for two possible reasons: (a) the terminal event was right-censored, i.e., the subject did not have a terminal event within the finite data-collection period, or (b) the terminal events are inherently unobservable. Right-censoring happens for instance when a patient is still alive at the time of data-collection. Unobservability happens for instance, in social networks, when a subject simply stops using the service but does not provide a clear termination signal by deleting her account. In such a scenario, the terminal events remain unobserved for most if not all subjects, even if the subjects quit the service within the data-collection period.
Numerous survival methods have been proposed (Witten and Tibshirani, 2010a; Hothorn et al., 2006; Ishwaran et al., 2008) to predict the lifetime of a subject given her covariates and her activities/measurements for a brief initial period of time, while also accounting for right-censoring. More recent deep learning models for lifetime prediction (Lee et al., 2018; Ren et al., 2018; Chapfuwa et al., 2018) have achieved much success due to their flexibility to model complex relationships, and by avoiding limiting assumptions like parametric lifetime distributions (Ranganath et al., 2016) and proportional hazards (Katzman et al., 2018). In scenarios where terminal events are never observed (unobservability), it is a standard practice to introduce artificial termination signals through a predefined “timeout” for the period of inactivity, i.e., a social network user inactive for m months has her last observed activity declared a terminal event. Such a specification is typically arbitrary and can adversely affect the analysis depending on the “timeout” value used.
Notwithstanding the fact that lifetimes are hard to predict without termination signals in the training data, we are generally interested in clustering subjects based on their underlying lifetime distribution to improve decision-making. Applications include identifying disease subtypes (Gan et al., 2018), understanding the implications of distinct manufacturing processes on machine parts, and qualitatively analyzing different survival groups in a social network. Although accurately predicting time-toterminal-event for an individual is important in a variety of applications, lifetime clustering plays a complementary role and provides a more holistic picture.
Lifetime clustering remains a relatively unexplored topic despite being an important tool. Although traditional unsupervised clustering methods such as k-means and hierarchical clustering are popular
for this task (Bhattacharjee et al., 2001; Sørlie et al., 2001; Bullinger et al., 2004), they may produce clusters that are entirely uncorrelated with lifetimes (Gaynor and Bair, 2013). Semi-supervised clustering (Bair and Tibshirani, 2004) and supervised sparse clustering (Witten and Tibshirani, 2010b) employ a two-stage lifetime clustering process: (i) identify covariates associated with lifetime using Cox scores (Cox, 1992), and (ii) treat these covariates differently while performing k-means clustering. They assume proportional hazards (i.e., constant hazard ratios over time) and require the presence of termination signals. Furthermore, a decoupled two-stage process such as the above is not guaranteed to obtain clusters with maximally distinct lifetime distributions; rather, we require an end-to-end learning framework that prescribes a loss function specifically over the lifetime distributions of different clusters.
In this work we tackle the important task of inductive lifetime clustering without assuming proportional hazards, while also smoothly handling the unobservability of termination signals.
Contributions. (i) We introduce DeepCLife, an inductive neural-network based lifetime clustering model that finds cluster assignments by maximizing the divergence between empirical (nonparametric) lifetime distributions of different clusters. Whereas the subjects of different clusters have distinct lifetime distributions, within a cluster all subjects share the same lifetime distribution even if they have different lifetimes. Our model is robust to the modeling issues associated with the unobservability of termination signals and does not assume proportional hazards. (ii) We define a novel clustering loss function over empirical lifetime distributions (of entire clusters) based on the Kuiper two-sample test. We provide a tight upper bound of the Kuiper p-value with easy-to-compute gradients facilitating its use as a loss function, which until now was prohibitively expensive due to the test’s infinite sum. (iii) Finally, our results on real and synthetic datasets show that the proposed lifetime clustering approach produces significantly better clusters with distinct lifetime distributions (as evaluated by Logrank score, C-index, Brier Score and Rand index) as compared to competing approaches.

heading: 2 FORMAL FRAMEWORK
text: In this section, we formally define the statistical framework underlying the clustering approach introduced later in the paper. Notation remark: We use superscript (u) to refer to variables indexed by a subject u and use subscript k to refer to random variables indexed by a random subject of cluster k. See Table 1 for a complete list of all the variables and their meanings.
We assume that there are distinct underlying clusters with different event processes describing the activity events of a random subject (e.g., logins to a social network, measurements of a patient) in the respective clusters. To make these notions more formal, we introduce our event process.
Definition 1 (Abstract Event Process). Consider the k-th cluster. The Random Marked Point Process (RMPP) for the activity events is Φk = {Xk, {(Ak,i,Mk,i, Yk,i)}i∈N,Θk} over discrete times t = 0, 1, . . ., where Xk is the random variable representing the covariates of a random subject in cluster k, Θk is the time to the zeroth activity event (joining), Mk,i represent covariates of event i, Yk,i is the inter-event time between the i-th and the (i − 1)-st activity events (e.g., logins), and Ak,i = 1 indicates an event with a termination signal (death), otherwiseAk,i = 0. All these variables may be arbitrarily dependent. This definition is model-free, i.e., we will not prescribe a model for Φk.
We observe the RMPP over a time window [0, tm]. Using the above formalism, we define the true lifetime of a subject in cluster k as the sum of all inter-event times until termination signal Ak,i = 1 is observed.
Definition 2 (True lifetime). The random variable that defines the true lifetime until the terminal event of a subject in cluster k is Tk := maxi (∑ i′≤i Yk,i′ ∏ i′′<i(1−Ak,i′′) ) .
The true lifetime Tk is unobserved if: (a) the terminal event Ak,i = 1 does not occur within the observation time period [0, tm], or (b) the termination signals are inherently unobservable. Now, the true lifetime distribution of a subject in cluster k is defined as the probability that the subject has at least one more activity event after time t, and is given by Sk(t) := P [Tk > t] = 1− Fk(t), t ∈ N∪{0}, where Fk(t) is the underlying cumulative distribution function (CDF) of Tk. Note again that all the subjects of a cluster share the same lifetime distribution.
Training data. Our training data D consists of subjects from K underlying (hidden) clusters with distinct lifetime distributions Sk for k ∈ {1, . . . ,K}. For a subject u ∈ D, we observe the following quantities {X(u), {(M (u)i , Y (u) i )} Q(u)(tm) i=1 ,Θ (u)}, where X(u), Y (u)i ,M (u) i and Θ
(u) are analogously defined as in Definition 1, but for a given subject u. Q(u)(t) is the number of observed activity events of u after her joining and before time t. The training data may or may not contain the termination signals, {A(u)i } Q(u)(tm) i=1 , where A (u) i = 1 indicates that event i was a terminal event for subject u (e.g., death). Termination signals are typically available in healthcare applications, whereas in the case of social networks, we usually do not observe the termination signal (i.e., account deletion) for any subject.
We define the observed lifetime of a subject u as H(u) := ∑Q(u)(tm) i=1 Y (u) i , i.e., the sum of all the inter-event times within the observation period [0, tm]. In the absence of termination signals, we additionally define inactive period as the time elapsed since the last observed event of subject u, given by χ(u) := tm −Θ(u) −H(u). We use H(u) and χ(u) only during training. Figure 1 shows the true lifetime, the observed lifetime, and the inactive period for four users of two different clusters in the training data. The events are observed (denoted by solid circles) only till the time of measurement tm, whereas the rest of the events are right-censored (denoted by solid diamonds). The termination signal (denoted by dotted diamonds) may be unobserved even if it occurs before tm (eg., u3).
Test data. For a new test subject, we can only observe her covariates and her activity events for a brief initial period of time τ since her joining. Formally, we observe the following quantities
{X(u), {(M (u)i , Y (u) i )} Q(u)(Θ(u)+τ) i=1 ,Θ
(u)}. The goal of inductive lifetime clustering is to be able to assign a cluster for a new unseen subject within τ time of her joining. Note that since H(u) and χ(u) are tied to tm, they cannot be computed for test subjects.
Lastly, we formally define our clustering problem. Definition 3 (Clustering problem). Consider a dataset D with N subjects. Our goal is to find a mapping κ : ( X(u), {(M (u)i , Y (u) i )} Q(u)(Θ(u)+τ) i=1 ) → {1, . . . ,K}, that inductively maps subject
covariates and observed activity events for a brief initial period of time τ into clusters, such that the divergence ∆ between the empirical lifetime distributions of these clusters is maximized, i.e.,
κ? = arg max κ∈K min i,j∈{1...K},
i 6=j
∆(Ŝi(κ), Ŝj(κ)) , (1)
where K is a set of all mappings, Ŝk(κ) is the empirical lifetime distribution of subjects in D mapped to cluster k through κ, and ∆ is an empirical distribution divergence measure.
κ∗ optimized in this fashion guarantees that subjects in different clusters have different lifetime distributions 1.
1See Appendix A.5 for why we wish to maximize the minimum divergence across all pairs instead of the sum of divergences across all pairs.

heading: 3 THE DEEPCLIFE MODEL
text: In this section, we propose a practical lifetime clustering approach using neural networks that optimizes Equation (1). LetD be the training data as defined in Section 2. Since we want to maximize divergence between empirical lifetime distributions, we assume discrete times (relative to subject joining), t ∈ {0, 1, . . . , tmax}, where tmax = maxu∈DH(u) is the maximum observed lifetime of any subject u ∈ D. Note that it is sufficient to define the empirical distribution till tmax, since we have not observed any subject with lifetime greater than tmax.
3.1 CLUSTER ASSIGNMENTS : α(u)k (W1)
We define a neural network g that takes user covariates and the event data for a subject u during a brief initial period τ after her joining as input, and outputs her cluster assignment probabilities, α
(u) k (W1) for all k ∈ {1, . . .K},
~α(u)(W1) = g ( X(u), {M (u)i , Y (u) i } Q(u)(Θ(u)+τ) i=1 ;W1 ) , (2)
where W1 are the weights of the neural network. The final layer of g is a softmax layer with K units. Figure 2a depicts g as a feedforward neural network with L− 1 hidden layers, although our model is not restricted to a feedforward architecture. In our experiments, we compute summary statistics over the observed events {M (u)i , Y (u) i }i in order to make it compatible with the feedforward architecture.
3.2 PROBABILITY OF TERMINATION: β(u)(W2)
During the training of our model, we require termination signals or a probabilistic estimation of the termination signals in order to write the likelihood of the model. Given the model parameter W2, we define β(u)(W2) as the probability that the last observed event of u was terminal, i.e., u will have no future activity events after the last observed event.
If the termination signals A(u)i are observed in the training data D (e.g., healthcare), clearly β(u)(W2) := A (u)
Q(u)(tm) (i.e., W2 is ignored).
When such termination signals are unobservable (e.g., social network), existing survival methods commonly use a timeout window of predefined size Wfixed over the inactive period χ(u) to specify the probability of termination as β(u)(Wfixed) := 1[χ(u)>Wfixed]. However, such specification is arbitrary and precludes any learning of the window size parameter Wfixed. Instead, we model the latent termination probabilities β(u)(W2) using a smooth non-decreasing function of χ(u), i.e., higher the period of inactivity, higher the probability that the last observed event was terminal.
If the termination signals A(u)i are unobservable in the training data D (e.g., social network), we use β(u)(W2) := 1− e−W2·χ (u) where W2 ≥ 0.
3.3 EMPIRICAL LIFETIME DISTRIBUTION OF CLUSTER k : Ŝk(W1,W2;D)
Given the training data D and model parameters W1 and W2, we can obtain the soft cluster assignments α(u)k (W1) and the probability of termination β
(u)(W2) for all subjects u ∈ D and clusters k ∈ {1 . . .K} as shown in Section 3.1 and Section 3.2. In this subsection, we obtain the empirical lifetime distribution of all the clusters k = 1 . . .K, using the Kaplan-Meier estimates (Kaplan and Meier, 1958). We do not assume a parametric form for the lifetime distribution, and rather use empirical distributions in our optimization to allow lifetime curves of any shape (Figure 2b). Kaplan-Meier estimates are a maximum likelihood estimate of the lifetime distribution of a set of subjects assuming (a) hard memberships (each subject entirely belongs to the set) and (b) the presence of termination signals. We modify the estimates to account for partial memberships and probability of termination instead.
Proposition 1. Given the training data D, a cluster k, the cluster assignment probabilities {α(u)k (W1)}u∈D, and the probabilities of termination {β(u)(W2)}u∈D, the maximum likelihood estimate of the empirical lifetime distribution of cluster k is given by,
Ŝk(W1,W2;D)[t] = t∏
j=0
sk(W1;D)[j]− dk(W1,W2;D)[j] sk(W1;D)[j] , (3)
for all t ∈ {0, 1, . . . , tmax}, where, sk(W1;D)[j] = ∑ u∈D 1[H
(u) ≥ j] · α(u)k (W1), is the expected number of subjects in cluster k who are at risk (of termination) at time j, and, dk(W1,W2;D)[j] =∑ u∈D 1[H
(u) = j] · β(u)(W2) · α(u)k (W1), is the expected number of subjects that are predicted to have had a terminal event at time j.
The proof is presented in the supplementary material.
3.4 EMPIRICAL DISTRIBUTION DIVERGENCE LOSS : ∆(Ŝa, Ŝb)
We rewrite the objective function of lifetime clustering from Definition 3 with respect to the model parameters W1 and W2 as follows,
W ∗1 ,W ∗ 2 = arg max
W1,W2
min i,j∈{1...K},
i 6=j
∆ ( Ŝi, Ŝj ) , (4)
where Ŝi is the empirical distribution, a shorthand for the vector Ŝi(W1,W2;D), and ∆ is a divergence measure between two empirical distributions.
We note the following essential requirements for the divergence measure ∆: (a) ∆ defined over empirical distributions must take into account sample sizes, (b) ∆ should have easy-to-compute gradients since it is used as an objective function to train neural networks, and (c) ∆ should not assume proportional hazards, and should allow for crossing lifetime curves (see Figure 2b).
Divergence measures such as Kullback-Leibler (Kullback and Leibler, 1951) and MMD (Gretton et al., 2012) fulfill (b, c) but not (a), and will result in sample anomalies as depicted in Figure 2c (also seen in our experiments: Figure 3d). Logrank test (Mantel, 1966), commonly used for comparing lifetime distributions, fulfills (a, b), but has low statistical power when the proportional hazards assumption is not met (Peto and Peto, 1972; Bland and Altman, 2004) (e.g., Figure 2b). Finally, p-value from two-sample tests such as the Kolmogorov-Smirnov (K-S) test (Massey Jr, 1951) fulfill (a, c), but not (b) as they require the computation of an infinite sum, resulting in an impractical objective function unless heuristic approximations are made.
We propose to use the Kuiper test (Kuiper, 1960), a two-sample test closely related to the K-S test with increased statistical power in distinguishing distribution tails (Tygert, 2010). Specifically, we define ∆(Ŝa, Ŝb):=− log(KD(Ŝa, Ŝb)), where KD is the p-value from the Kuiper test between Ŝa and Ŝb. Our choice of the Kuiper test is because it is amenable to upper and lower bounds as shown next, thus avoiding the prohibitive heuristic approximations of infinite sums.
Proposition 2. (Bounds for Kuiper p-value.) Given two empirical lifetime distributions Ŝa and Ŝb with discrete support and sample sizes na and nb respectively, define the maximum positive and negative separations between them,
D̂+a,b = sup t∈{0,...tmax} (Ŝa[t]− Ŝb[t]) , D̂−a,b = sup t∈{0,...tmax} (Ŝb[t]− Ŝa[t]) .
The Kuiper test p-value Kuiper (1960) gives the probability that Λ, the empirical deviation for na and nb observations under the null hypothesis Sa = Sb 2, exceeds the observed value V = D̂+a,b + D̂ − a,b:
KD(Ŝa, Ŝb) ≡ P [Λ > V ] = 2 ∞∑ j=1 (4j2λ2a,b − 1)e−2j 2λ2a,b , (5)
λa,b = ( √ Ma,b + 0.155 +
0.24√ Ma,b )V and Ma,b = nanbna+nb is the effective sample size. Then, the
upper bound 3 for the Kuiper p-value is, KD(Ŝa, Ŝb) ≤min ( 1, 2 · 1[r(lo)a,b ≥ 1] · ( w(r (lo) a,b , λa,b)− w(1, λa,b) + v(r (lo) a,b , λa,b) ) + v(r
(up) a,b , λa,b)− w(r (up) a,b , λa,b)
) , (6)
where v(r, λ)=(4r2λ2−1)e−2r2λ2 , w(r, λ)=−re−2r2λ2 , r(lo)a,b = b 1√2λa,b c, and, r (up) a,b = d 1√2λa,b e .
The proof is presented in the supplementary material.
Figures 2d and 2e show empirical results of the upper and lower bounds as a function of V = D+ + D− against the expensive heuristic computation where the infinite sum in Equation (5) is approximated with 104 terms (Press et al., 1996). We optimize Equation (4) with ∆(Ŝa, Ŝb) := − log(KD(up)(Ŝa, Ŝb)) where KD(up) denotes the Kuiper p-value upper bound in Equation (6). ∆ defined this way satisfies all three requirements we described in the beginning of Section 3.4: takes sample sizes into account, has closed form expression with easy-to-compute gradients, and does not assume proportional hazards.
2The test is typically defined using CDFs but it is equivalent to use lifetime distributions (CCDFs) instead. 3Lower bound is presented in the supplementary material.
Implementation. We implement our clustering procedure as a feedforward neural network in Pytorch (Paszke et al., 2017) and use ADAM (Kingma and Ba, 2014) to optimize Equation (4). Each iteration of the optimization takes as input a batch of subjects, generates a single value for the set loss, calculates the gradients, and updates the parameters W1 and W2. We show in the supplementary material that the time and space complexity per iteration of the proposed approach are O(KB +K2tmax) and O(Ktmax) respectively, where K is the number of clusters and B is the batch size. For large values of K, we achieve tractability by sampling K random pairs of clusters every iteration instead of all ( K 2 ) possible pairs.

heading: 4 RELATED WORK
text: Majority of the work in survival analysis has dealt with the task of predicting the time to an observable terminal event (e.g., death), especially when the number of features is much larger than the number of subjects (Witten and Tibshirani, 2010a; pre; Hothorn et al., 2006; Shivaswamy et al., 2007). Recently, many deep learning approaches (Luck et al., 2017; Katzman et al., 2018; Lee et al., 2018; Ren et al., 2018) have been proposed for predicting the lifetime distribution of a subject given her covariates, while effectively handling censored data that typically arise in survival tasks. DeepHit (Lee et al., 2018) introduced a novel architecture and a ranking loss function in addition to the log-likelihood loss for lifetime prediction in the presence of multiple competing risks. Using a log-likelihood loss similar to DeepHit, Ren et al. (2018) propose a recurrent architecture to predict the survival distribution that captures sequential dependencies between neighbouring time points. Finally, Chapfuwa et al. (2018) introduce an adversarial learning framework to model the lifetime given the subject covariates. In contrast to these works on predicting lifetimes, our task is to cluster the subjects based on their underlying lifetime distributions.
There are relatively fewer works that perform lifetime clustering. Many unsupervised approaches have been proposed to identify cancer subtypes in gene expression data but do not consider the lifetime (Eisen et al., 1998; Alizadeh et al., 2000; Bhattacharjee et al., 2001; Sørlie et al., 2001; Bullinger et al., 2004), and may produce clusters that are entirely independent of the lifetimes. Semi-supervised clustering (Bair and Tibshirani, 2004) and supervised sparse clustering (Witten and Tibshirani, 2010b) use Cox scores (Cox, 1992) to identify features associated with the lifetime and treat these features differently while using k-means to perform the final clustering. Unlike these lifetime clustering methods, DeepCLife does not assume proportional hazards, and can smoothly handle the absence of termination signals. Our supplementary material has a more in-depth discussion of related work.

heading: 5 RESULTS
text: Baselines. We perform experiments on one synthetic dataset and two real-world datasets – Friendster social network and MIMIC III healthcare dataset (Johnson et al., 2016).
We compare the following lifetime clustering approaches: (a) SSC-Bair, a semi-supervised clustering method (Bair and Tibshirani, 2004) that performs k-means clustering on selected covariates that have high Cox scores (Cox, 1992); (b) SSC-Gaynor or supervised sparse clustering (Gaynor and Bair, 2013), a modification of sparse clustering (Witten and Tibshirani, 2010b) that weights the covariates based on their Cox scores; (c) DeepHit+GMM, a Gaussian mixture model applied over last layer embeddings learnt by DeepHit (Lee et al., 2018); (d) DeepCLife-MMD, the DeepCLife model with Maximum Mean Discrepancy (Gretton et al., 2012) as the divergence measure; (e) DeepCLifeKuiperUB, the DeepCLife model with the proposed divergence measure based on the Kuiper p-value upper bound (Equation (6)).
Termination signals for evaluation and baselines. (Timeout.) The competing methods and most evaluation metrics for survival applications require clear termination signals. In scenarios where we do not observe termination signals (e.g., Friendster experiments), we specify termination signals artificially when training the baselines using a pre-defined “timeout”, i.e., β(u)(Wfixed) = 1[χ(u) > Wfixed]. During evaluation, we use the same Wfixed to specify the termination signals and compute the metrics. This helps the competing methods since they are trained and evaluated using termination signals with the same Wfixed, whereas our approach is not trained with these termination signals.
Metrics. Given the cluster assignments κ(u′) ∈ {1, . . . ,K} and the termination signals (possibly using a “timeout” window) for all the users u′ in the test data, we can obtain the empirical lifetime distribution of all the clusters Ŝk, ∀k ∈ {1, . . . ,K} using the Kaplan-Meier estimates (over the test data alone). Then, the empirical lifetime distribution of a user u′ is given by that of her assigned cluster, i.e., Ŝ(u
′) := Ŝκ(u′). We use the following metrics4 for evaluating the clusters obtained from the methods.
• Logrank Score ↑. Logrank test (Mantel, 1966) statistic is a non-parametric test that outputs high values when it is unlikely for the K groups to have the same lifetime distribution. • Adjusted Rand index ↑. The Adjusted Rand index (Hubert and Arabie, 1985) is a measurement of cluster agreement, compared to the ground truth clustering (if available). ARI is 0.0 for random cluster assignments and 1.0 for perfect assignments. • C-Index ↑. Concordance index (Harrell et al., 1982) is a commonly used metric that calculates the fraction of pairs of subjects for which the model predicts the correct order of survival while also incorporating censoring. We use the expected lifetime obtained from the lifetime distribution Ŝ(u
′) as the predicted lifetime of user u′. C-index is 1.0 for perfect predictions and 0.5 for random predictions. • Integrated Brier Score ↓. Integrated Brier score (Brier, 1950; Graf et al., 1999) computes mean squared difference between the survival probabilities and the actual outcome over [0, tmax]. It ranges from 0.25 for random predictions to 0 for perfect predictions.
Although cluster evaluation metrics like Logrank score are more suitable for the lifetime clustering task, we also use predictive measures like C-index and Brier Score to further validate the clusters.
Evaluation. We evaluate the models using 5-fold cross validation. We use the ith fold for testing and sample N (tr) subjects from the remaining 4 folds for training. We use 20% of the N (tr) subjects as validation for early stopping and hyperparameter tuning for the different approaches.

heading: 5.1 EXPERIMENTS
text: Synthetic experiment. We test our method with a synthetic dataset6 for which we have the true clusters as ground truth. We generate 3 clusters C1, C2 and C3 with different lifetime distributions such that C2 and C3 have proportional hazards, but lifetime distribution of C1 crosses the other two curves (shown in Figure 2b). We choose an arbitrary time of measurement, tm=150, to imitate right censoring. We sample 104 subjects for each cluster, their features drawn from mixture of Gaussians. The lifetime T (u) of a subject u is randomly sampled from the lifetime distribution of her ground truth cluster. Table 2 shows performance of the methods on three synthetic datasets: D{C1,C2}, D{C1,C3},
4↑ indicates higher is better, ↓ indicates lower is better. Detailed description in the supplementary material. 4The values in bold are statistically significant with p-value < 0.01 using paired t-test over validation folds. 5More details about all the datasets and preprocessing steps can be found in the supplementary material.
D{C1,C2,C3}. DeepCLife-KuiperUB and DeepCLife-MMD were able to recover perfect ground truth cluster assignments on D{C1,C2} and D{C1,C3}, whereas the baseline methods performed invariably worse. On the harder dataset D{C1,C2,C3}, DeepCLife-KuiperUB recovered the ground truth clusters almost twice as better than any other method.
Friendster experiment. Friendster dataset5 consists of 15 million users in the Friendster online social network along with the comments sent and received by the users. We consider a subset with 1.1 million users who had participated in at least one comment. We use each user’s profile information (like age, gender, location, etc.) as covariates, and define activity events as the comments sent or received by the user. We compute summary statistics of the activity events within τ = 5 months from joining such as number of comments sent/received, number of people interacted with, mean inter-event times, etc. The task is to cluster new test users using their covariates and the summary statistics computed over initial τ = 5 months of activity. Note that we do not observe termination signals (i.e., account deletion) for any subject in the data. We use an arbitrary window of Wfixed = 10 months over the inactivity period to obtain termination signals (≈ 65% assumed to have quit) for the competing methods and for computing the evaluation metrics; DeepCLife does not require such arbitrary specification but learns a smooth timeout window during training.
Table 3 shows results for K=2, 4 clusters. We note that the proposed method obtains higher C-index values and Brier scores compared to the baselines even without termination signals. DeepCLife-
KuiperUB achieves a significant improvement in Logrank scores compared to the baselines because its loss specifically maximizes for differences in empirical distributions, while also taking sample sizes into account. DeepCLife-MMD on the other hand does not account for sample sizes, and hence performs worse. The empirical lifetime distributions of the clusters obtained from different methods for K=3 are shown in Figure 3. The clusters obtained from the baselines (a-c) are not substantially different from each other. Although DeepCLife-MMD obtains clusters with distinct lifetime distributions, it suffers from sample anomalies i.e., outputs clusters with very few or no subjects (e.g., Ŝ2 in Figure 3d). DeepCLife-KuiperUB outputs clusters that have significantly different lifetime distributions with the best Logrank scores. Corresponding plots of empirical lifetime distributions for K=2, 4, 5 are presented in the supplementary material. For K=4, we observe that DeepCLife-KuiperUB finds crossing yet distinct lifetime distributions (see Figure 7e in the supplementary material). We present a heuristic to select optimal K based on the validation performance of the clusters in Appendix A.4.
Qualitative Analysis: In the clusters found by DeepCLife-KuiperUB in Friendster, we see that a user in a low-risk cluster has on average 7.76 friends, sends 5.06 comments with an average response time of 20 days. On the other hand, a user in a high-risk cluster has just 1.56 friends on average and sends far fewer comments, around 1.07, but with a fast response time of 1.32 days. Interestingly, users that stayed longer in the system had lower activity rate in the beginning.
MIMIC III experiment. MIMIC III dataset5 (Johnson et al., 2016) consists of around 46500 patients admitted to the Intensive Care Unit (ICU). The task is to cluster the patients based on their mortality within 30 days of admission to the ICU (Purushotham et al., 2017). Unlike Friendster experiments, we can observe terminal events. Lifetime of a patient is right-censored if she was discharged within the 30-day period (≈ 84% right-censored). We use only the initial τ=24 hours of patient measurements (e.g. heart rate, respiratory rate, etc.) to perform the clustering. We compute summary statistics for each type of measurement, for example, number of heart rate measurements taken, mean heart rate, etc. to make the inputs compatible with the feedforward architecture in Figure 2a. Results for K=2 in Table 4b show that DeepCLife-KuiperUB achieves significantly better C-index and Logrank scores, followed by DeepHit+GMM. The corresponding empirical lifetime distributions of the clusters are shown in Figure 4a.

heading: 6 CONCLUSION
text: In this work we introduced Kuiper-based nonparametric loss function to maximize the divergence between empirical distributions, and a corresponding upper bound with easy-to-compute gradients. The loss function is then used to train a feedforward neural network to inductively map subjects into K lifetime-based clusters without requiring termination signals. We show that this approach produces clusters with better C-index values and Logrank scores than competing methods.

heading: A APPENDIX

heading: A.1 PROOF OF PROPOSITION 1
text: Proposition (Proposition 1). Given the training data D, a cluster k, the cluster assignment probabilities {α(u)k (W1)}u∈D, and the probabilities of termination {β(u)(W2)}u∈D, the maximum likelihood estimate of the empirical lifetime distribution of cluster k is given by,
Ŝk(W1,W2;D)[t] = t∏
j=0
sk(W1;D)[j]− dk(W1,W2;D)[j] sk(W1;D)[j] ,
for all t ∈ {0, 1, . . . , tmax}, where, sk(W1;D)[j] = ∑ u∈D 1[H
(u) ≥ j] · α(u)k (W1), is the expected number of subjects in cluster k who are at risk (of termination) at time j, and, dk(W1,W2;D)[j] =∑ u∈D 1[H
(u) = j] · β(u)(W2) · α(u)k (W1) , is the expected number of subjects that are predicted to have had a termination event at time j.
Proof. For a subject u, let e(u) ∼ Bernoulli(β(u)(W2)) be the sampled termination signal, i.e., e(u) = 1 indicates that the subject’s observed lifetime is her true lifetime, H(u) = T (u). Let ~κ(u) ∼ Categorical(α(u)1 (W1), . . . , α (u) K (W1)) be the sampled cluster assignment of the subject.
Finally, let the empirical lifetime PMF of a cluster k ∈ {1, 2, . . .K} be defined as P (Tk = t) := λk,t for t ∈ {0, 1, . . . tmax}, where Tk is the true lifetime of a random subject in cluster k. The likelihood of the training data D for a particular cluster k given the sampled termination signal {e(u)}u∈D and the sampled cluster assignments {~κ(u)}u∈D can be written using Kaplan and Meier (1958) as,
Lk = [∏ u∈D P (Tk = H (u))e (u) P (Tk > H (u))1−e (u) ]κ(u)k .
Next, we take expectation over the cluster assignments and termination signals to get,
Lk = E{~κ(u)}uE{e(u)}u [∏ u∈D P (Tk = H (u))e (u) P (Tk > H (u))1−e (u) ]κ(u)k . (7)
Using mean-field approximation and proceeding as Kaplan and Meier (1958) to write Equation (7) as a product over time instead of subjects, we get,
Lk = tmax∏ j=0 λ dk[j] k,j (1− λk,j)sk[j]−dk[j]  , (8) where, we have used P (Tk = j) = λk,j , sk[j] is a shorthand for sk(W1,W2;D)[j] =∑
u∈D 1[H (u) ≥ j] · α(u)k (W1), the expected number of subjects at risk at time j, and, dk[j] is a shorthand for dk(W1,W2;D)[j] = ∑ u∈D 1[H
(u) = j] · β(u)(W2) · α(u)k (W1), the expected number of subjects that are predicted to have had a termination event at time j.
Finally, maximizing the likelihood in Equation (8) with respect to λk,t for all k ∈ {1, . . .K} and t ∈ {0, 1, . . . tmax}, and computing the lifetime distributions (CCDF) Ŝk from the corresponding lifetime PMFs finishes the proof.

heading: A.2 PROOF OF PROPOSITION 2
text: Proposition (Proposition 2). Given two empirical lifetime distributions Ŝa and Ŝb with discrete support and sample sizes na and nb respectively, define the maximum positive and negative separations between them as,
D̂+a,b = sup t∈{0,...tmax} (Ŝa[t]− Ŝb[t]) ,
D̂−a,b = sup t∈{0,...tmax} (Ŝb[t]− Ŝa[t]) .
The Kuiper test p-value Kuiper (1960) gives the probability that Λ, the empirical deviation for na and nb observations under the null hypothesis Sa = Sb , exceeds the observed value V = D̂+a,b + D̂ − a,b:
KD(Ŝa, Ŝb) ≡ P [Λ > V ] = 2 ∞∑ j=1 (4j2λ2a,b − 1)e−2j 2λ2a,b , (9)
where λa,b = (√
Ma,b + 0.155 + 0.24√ Ma,b ) V and Ma,b = nanbna+nb is the effective sample size.
Then, the lower and upper bounds for the Kuiper p-value are, KD(Ŝa, Ŝb)
2 ≥ 1[r(lo)a,b ≥ 1] ·
( w(r
(lo) a,b − 1, λa,b)
+ v(r (lo) a,b , λa,b)
) + v(r
(up) a,b , λa,b) + w(r (up) a,b + 1, λa,b) ,
and KD(Ŝa, Ŝb)
2 ≤ min
( 1
2 ,1[r
(lo) a,b ≥ 1] ·
( w(r
(lo) a,b , λa,b)
− w(1, λa,b) + v(r(lo)a,b , λa,b) )
+ v(r (up) a,b , λa,b)− w(r (up) a,b , λa,b)
) , (10)
where v(r, λ) = (4r2λ2 − 1)e−2r2λ2 , with w(r, λ) = −re−2r2λ2 , and r(lo)a,b = b 1√2λa,b c, r (up) a,b = d 1√ 2λa,b e .
Proof. To prove these bounds, we first regard that the indefinite integral w(r, λ) = ∫ v(r, λ)dr = −re−2r2λ2 + C. Then, for a given λ, v(r, λ) is monotonically increasing w.r.t r in the interval (0, 1√
2λ ) and monotonically decreasing w.r.t r in the interval ( 1√ 2λ ,∞).
Let r(lo) = b 1√ 2λ c, r(up) = d 1√ 2λ e and J[l,m](λ) = ∑m r=l v(r, λ). Here, we assume that r
(lo) 6= r(up), i.e., 1√
2λ is not an integer. The bounds can be easily modified for when r(lo) = r(up). Now,
v(r, λ) is decreasing w.r.t r in (r(up),∞), and we have,∫ ∞ r(up)+1 v(r, λ)dr ≤ J[r(up)+1,∞](λ) ≤ ∫ ∞ r(up) v(r, λ)dr , and adding v(r(up), λ) throughout we get,
v(r(up), λ) + ∫ ∞ r(up)+1 v(r, λ)dr
≤ J[r(up),∞](λ) ≤ v(r(up), λ) + ∫ ∞ r(up) v(r, λ)dr . (11)
Similarly, since v(r, λ) is increasing in (0, r(lo)), we have,∫ r(lo)−1 0 v(r, λ)dr ≤ J[1,r(lo)−1](λ) ≤ ∫ r(lo) 1 v(r, λ)dr , and adding v(r(lo), λ) throughout we get,∫ r(lo)−1 0 v(r, λ)dr + v(r(lo), λ)
≤ J[1,r(lo)](λ) ≤ ∫ r(lo)
1
v(r, λ)dr + v(r(lo), λ) , (12)
where we assume that r(lo) ≥ 1, otherwise, J[1,r(lo)](λ) is trivially zero.
Combining equations equation 11 and equation 12, we get, 1[r(lo) ≥ 1] · (∫ r(lo)−1
0
v(r, λ)dr + v(r(lo), λ)
)
+v(r(up), λ) + ∫ ∞ r(up)+1 v(r, λ)dr
≤ 1[r(lo) ≥ 1] · J[1,r(lo)](λ) + J[r(up),∞](λ) ≤ 1[r(lo) ≥ 1] · (∫ r(lo)
1
v(r, λ)dr + v(r(lo), λ)
)
+ v(r(up), λ) + ∫ ∞ r(up) v(r, λ)dr . (13)
We use the fact that KD(Ŝa, Ŝb) = 2 · J[1,∞](λa,b), J[1,∞](λ) = 1[r(lo) ≥ 1] · J[1,r(lo)](λ) + J[r(up),∞](λ), and w(r, λ) = ∫ v(r, λ)dr, and that KD(Ŝa, Ŝb) ≤ 1 to complete the proof.

heading: A.3 ADDITIONAL RESULTS

heading: A.3.1 FRIENDSTER EXPERIMENT
text: Results on the Friendster experiments for K = 3, 5 are shown in Table 4. We observe that DeepCLife outperforms all baselines by large margins, particularly in the Logrank score. Figures 5 to 8 show empirical lifetime distributions of the clusters found by the methods on the Friendster dataset for K = 2, 3, 4, 5 respectively. Baseline methods SSC-Bair, SSC-Gaynor and DeepHit+GMM employ a two-stage clustering process and do not guarantee clusters with maximally different lifetime distributions. DeepCLife-MMD does not account for sample sizes of empirical distributions and suffers from sample anomalies, i.e., outputs clusters containing almost no subjects. DeepCLifeKuiperUB obtains clusters with the best Logrank scores. Interestingly, DeepCLife finds clusters with crossing yet distinct lifetime distributions for K = 4 (see Figure 7e).

heading: A.3.2 MIMIC III EXPERIMENT
text: Figure 9 shows the empirical lifetime distributions of the clusters found by DeepCLife-KuiperUB in MIMIC III dataset for K = 2 . . . 6.

heading: A.4 OPTIMAL NUMBER OF CLUSTERS : K
text: Similar to standard clustering algorithms, choice of K is subjective and depends on the downstream application of the clusters. Here we discuss a heuristic to use Logrank scores as a useful guide to choose the optimal number of clusters. For the Friendster experiments, Figure 10 shows the Logrank
Table 4: (Friendster) C-index (%), Integrated Brier Score (%) and Logrank score with standard errors in parentheses for different methods3 and K=3, 5 clusters with number of training examples N (tr)=105.
K = 3 K = 5
Method C-index ↑ (%) I.B.S ↓ (%) Logrank Score ↑ C-index ↑ (%) I.B.S. ↓ (%) Logrank Score ↑ SSC-Bair 64.23 (0.24) 22.19 (0.02) 5277.82 ( 43.63) 75.80 (0.16) 20.26 (0.02) 43605.16 ( 91.39) SSC-Gaynor 64.13 (0.21) 22.18 (0.02) 5400.75 ( 39.59) 75.66 (0.20) 20.48 (0.01) 41824.28 ( 151.85) DeepHit+GMM 74.81 (0.11) 20.97 (0.06) 33880.55 ( 426.74) 76.17 (0.27) 20.54 (0.11) 39254.67 (1721.56) DeepCLife-MMD 72.75 (1.24) 18.94 (0.65) 50173.94 (4809.80) 63.08 (2.69) 21.55 (0.54) 14308.90 (9149.41) DeepCLifeKuiperUB 78.48 (0.32) 19.70 (0.55) 54191.28 ( 436.72) 76.51 (1.25) 18.99 (0.16) 56130.21 (4036.02)
0 10 20 30 40 50 Time
0.0
0.2
0.4
0.6
0.8
1.0
Pr ob
ab ilit
y (C
CD F)
S1 (n1 = 123839) S2 (n2 = 67794) S3 (n3 = 20660) S4 (n4 = 8121) S5 (n5 = 6448)
(a) SSC-Bair
0 10 20 30 40 50 Time (in months)
0.0
0.2
0.4
0.6
0.8
1.0
Pr ob
ab ilit
y (C
CD F)
S1 (n1 = 118101) S2 (n2 = 8023) S3 (n3 = 10755) S4 (n4 = 59562) S5 (n5 = 30421)
(b) SSC-Gaynor
0 10 20 30 40 50 Time (in months)
0.0
0.2
0.4
0.6
0.8
1.0
Pr ob
ab ilit
y (C
CD F)
S1 (n1 = 29583) S2 (n2 = 70151) S3 (n3 = 41885) S4 (n4 = 24213) S5 (n5 = 61030)
(c) DeepHit+GMM
0 10 20 30 40 50 Time (in months)
0.0
0.2
0.4
0.6
0.8
1.0
Pr ob
ab ilit
y (C
CD F)
S1 (n1 = 662) S2 (n2 = 66) S3 (n3 = 32809) S4 (n4 = 2781) S5 (n5 = 190544)
(d) DeepCLife-MMD
0 10 20 30 40 50 Time (in months)
0.0
0.2
0.4
0.6
0.8
1.0
Pr ob
ab ilit
y (C
CD F)
S1 (n1 = 31439) S2 (n2 = 41574) S3 (n3 = 29174) S4 (n4 = 39560) S5 (n5 = 85115)
(e) DeepCLife-KuiperUB
Figure 8: (Friendster) Empirical lifetime distributions of clusters obtained from different methods for K=5 (legend shows cluster sizes).
scores on a held-out validation set for clusters obtained from DeepCLife-KuiperUB with K = 2 . . . 7. We observe that K = 4 gives the best mean Logrank score and K = 5 is close behind. However, for K > 5, the Logrank scores drop drastically indicating that there are no more clusters that have distinct lifetime distributions from the ones already found.
Similarly for the MIMIC III experiments, Figure 11 shows the Logrank scores on a held-out validation set for clusters obtained from DeepCLife-KuiperUB with K = 2 . . . 8. We see that choosing K = 6 gives the highest Logrank score. However, we also observe that from K = 2 to K = 5 the scores do not increase significantly; a practitioner might decide to use K = 2 for this reason.

heading: A.5 CHOICE OF MINIMUM VS SUM IN EQUATIONS (1) AND (4)
text: We define the clustering problem in Definition 3 as finding a mapping κ such that,
κ? = arg max κ∈K min i,j∈{1...K},
i 6=j
∆(Ŝi(κ), Ŝj(κ)) ,
whereK is a set of all mappings, Ŝk(κ) is the empirical lifetime distribution of subjects in the training data D mapped to cluster k through κ, and ∆ is an empirical distribution divergence measure. In this section, we discuss the issues that arise when one chooses to optimize the sum of divergences over all pairs of clusters rather than the minimum.
Let K = 4 and consider balanced clusters with lifetime distributions Ŝ1 = Ŝ2 6= Ŝ3 = Ŝ4. Note that such a clustering is not desirable since there are virtually only 2 clusters found (as Ŝ1 and Ŝ2 are the same, etc). But for such an assignment, the sum of divergences is relatively high∑ i 6=j ∆(Ŝi, Ŝj) = 4 · ∆(Ŝ1, Ŝ3). The sum could be further maximized by simply increasing
divergence between Ŝ1 and Ŝ3 while keeping Ŝ2 = Ŝ1 and Ŝ3 = Ŝ4. On the other hand, since mini 6=j ∆(Ŝi, Ŝj) = 0, the model should maximize the divergence between Ŝ1 and Ŝ2 first.
Figure 12 shows the clusters found by the proposed approach while maximizing the minimum divergence vs maximizing the sum of divergences in the Friendster dataset with K = 4. We observe that the clusters obtained using the sum of divergences as the loss exhibit the property discussed above (i.e., clusters 2 and 3 in Figure 12 are almost the same).

heading: A.6 ALGORITHM
text: For ease of use, we present Algorithm 1 to compute the proposed divergence measure log(KD(up)(Ŝa, Ŝb)).

heading: A.7 TIME/SPACE COMPLEXITY OF DEEPCLIFE-KUIPERUB
text: Here we discuss the time and space complexity of the proposed algorithm. We only discuss the complexity of computation of the proposed loss function for one iteration assuming a minibatch of size B, number of clusters K and discrete times with 0, 1, . . . tmax. First, we can compute the empirical lifetime distributions of all the K clusters in O(KB) time using O(Ktmax) space. Computing upper bound of the Kuiper p-value for two lifetime distributions using Algorithm 1 requires O(tmax) time and O(tmax) space. Since we compute the upper bound for all ( K 2 ) pairs, the time complexity then becomes O(K2tmax) whereas the space complexity remains O(tmax). The total time and space complexity of the proposed algorithm are O(KB +K2tmax) and O(Ktmax) respectively.
Algorithm 1 Kuiper Divergence Upper Bound (Log)
Require: Two empirical lifetime distributions Ŝa and Ŝb with sample sizes na and nb respectively
1: function KUIPERUB(Ŝa, Ŝb, na, nb) 2: D+ ← max(Ŝa − Ŝb) . Ŝa, Ŝb are vectors of length tmax+1 3: D− ← max(Ŝb − Ŝa) 4: V ← D+ +D− . Kuiper statistic 5: M ← nanbna+nb . Effective sample size 6: λ← ( √ M + 0.155 + 0.24√
M )V
7: r(lo) ← b 1√ 2λ c 8: r(up) ← d 1√ 2λ e
9: if r(lo) ≥ 1 then 10: KD-UB←W(r(lo), λ)−W(1, λ) + V(r(lo), λ) + V(r(up), λ)−W(r(up), λ) 11: else 12: KD-UB← V(r(up), λ)−W(r(up), λ) 13: return log KD-UB
14: function V(r, λ) 15: return (4r2λ2 − 1)e−2r2λ2
16: function W(r, λ) 17: return −re−2r2λ2

heading: A.7.1 TRACTABILITY THROUGH SAMPLING
text: Each iteration of optimization involves computing the Kuiper upper bound for all ( K 2 ) pairs of clusters, hence can be prohibitively expensive for large values of K. Rather, we propose to sample randomly p = O(K) pairs without replacement from the ( K 2 ) possible pairs of clusters. The time and space complexity reduces to O(KB +Ktmax) and O(Ktmax) respectively.

heading: A.8 ADDITIONAL DETAILS ABOUT DATASETS
text: Synthetic dataset. We also test our method with a synthetic dataset for which we have the true clustering as ground truth. We generate 3 clusters C1, C2 and C3 with different lifetime distributions, as depicted in Figure 13. Note that while the lifetime distributions C2 and C3 follow proportional hazards assumption, the lifetime distribution of C1 violates this assumption by crossing the other survival curves.
We simulate 104 subjects for each cluster k, and generate 20 random features for each sampled subject. The lifetime T (u) of a subject u in a cluster was randomly sampled from the corresponding
lifetime distribution of the cluster. We also choose an arbitrary time of measurement, tm = 150, to imitate right censoring.
We generate covariates X(u) of each subject u as two sets of 10 attributes. Each attribute is simulated from a mixture of three Gaussians with means µ(1)k,i , µ (2) k,i , and µ (3) k,i , for i = 1 . . . 20 and k = 1, 2, 3. For the first set of 10 features, these means are uniformly selected from the range [0, 30] and a variance of 1 is used. For the second set of features, the three means are uniformly selected from the range [0, 30] and a variance of 10 is used. Now that the modes of each Gaussian of each feature are defined, we generate the features as follows: for each feature i of user u of cluster k, we uniformly choose one of the three modes m ∈ {1, 2, 3} . Once the mode is selected, we generate a value from N (µ(m)k,i , σ2) (where σ2 is 1 for i ∈ [1, 10] and 10 for i ∈ [11, 20]). We decided to have these two sets of features as the first one (with lower variance) represents more relevant features which can help identify the clusters, while the second set (with higher variance) represents less meaningful features that may be similar for all clusters.
We create three separate datasets D{C1,C2}, D{C1,C3} and D{C1,C2,C3} such that DC is a union of all the clusters in the set C. We report the C-index and Adjusted Rand index for the proposed method and the baselines on these 3 simulated datasets in Table 2.
Friendster dataset. Friendster dataset consists of around 15 million users with 335 million friendship links in the Friendster online social network. Each user has profile information such as age, gender, marital status, occupation, and interests. Additionally, there are user comments on each other’s profile pages with timestamps that indicate activity in the site.
In our experiments, we only use data from March 2002 to March 2008, as after March 2008 Friendster’s monthly active users have been significantly affected with the introduction of “new Facebook wall” (Ribeiro and Faloutsos, 2015). From this, we only consider a subset of 1.1 million users who had participated in at least one comment, and had specified their basic profile information like age and gender. We will make our processed data available to the public at location (anonymized). We use each user’s profile information (like age, gender, relationship status, occupation and location) as features. We convert the nominal attributes to one-hot encoded features. We compute the summary statistics of the user’s activity over the initial τ=5 months such as the number of comments sent and received, number of individuals interacted with, etc. In total, we construct 60 numeric features that are used for each of the models in our experiments.
MIMIC III dataset. MIMIC III dataset consists of patients compiled from two ICU databases: CareVue and Metavision. The dataset consist of 50,000 patients with a total of 330 million different measurements recorded. We mainly use three tables in the dataset: PATIENTS, ADMISSIONS, CHART_EVENTS. PATIENTS table consists of demographic information about every patient such as gender, date of birth, and ethnicity. ADMISSIONS table consists of the admission and discharge times of the patients. A patient can have multiple admissions; in our experiments, we only consider their longest stay. Finally, CHART_EVENTS table records all the measurements of the patients taken during their stay such as daily weight, heart rate and respiratory rate. The measurements of the patients are shifted by an unknown offset (consistent across measurements for a single patient) for anonymity, hence all times are relative to the patient.
Each entry in the CHART_EVENTS table has a column for item identifier that specifies the measured item (like heart rate). Since the dataset is compiled from different databases, same item can have multiple identifiers. We use the following time series for features: ‘Peak Insp. Pressure’, ‘Plateau Pressure’, ‘Respiratory Rate’, ‘Heart Rate’, ‘Mean Airway Pressure’, ‘Arterial Base Excess’, ‘BUN’, ‘Creatinine’, ‘Magnesium’, ‘WBC’ and ‘Hemoglobin’. We observe each of these time series for τ = 24 hours from the patient’s admission to the ICU, and compute summary statistics such as number of observations, mean, variance, mean inter-arrival times, etc. Complete processing script is provided in the supplementary material.

heading: A.9 METRICS
text: Here we describe the metrics used for evaluating the clusters obtained from the methods.
Concordance Index. Concordance index or C-index (Harrell et al., 1982) is a commonly used metric in survival applications (Alaa and van der Schaar, 2017; Luck et al., 2017) to quantify a model’s ability to discriminate between subjects with different lifetimes. It calculates the fraction of pairs of subjects for which the model predicts the correct order of survival while also incorporating censoring. Concordance index can be seen as a generalization of AUROC (Area Under ROC curve) and can be interpreted in a similar fashion, i.e., a C-index score of 1 denotes perfect predictions as compared to a C-index score of 0.5 for random predictions.
Integrated Brier Score. Brier score (Brier, 1950; Graf et al., 1999) is a quadratic scoring rule, computed as follows for survival applications,
B(u) = 1
tmax + 1 tmax∑ t=0 (1[T (u) > t]− Ŝ(u)[t])2 , (14)
where for a given subject u, T (u) is her true lifetime and Ŝ(u) is her empirical lifetime CCDF. Brier score is 0 for perfect predictions and 0.25 for random predictions.
Logrank Test Score. Logrank test Mantel (1966) is a non-parametric hypothesis test that is used to compare survival distributions. The null hypothesis for the test is defined as H0 : f1(t) = . . . = fK(t), i.e., the K groups have the same survival distribution. It is most appropriate when the observations are censored and the censoring is independent of the events. High values of the logrank statistic denote that it is unlikely that the K groups have the same survival distribution.
Adjusted Rand index. The Rand index (Rand, 1971) is a measurement of cluster agreement, compared to the ground truth clustering (if available). The adjusted version (Hubert and Arabie, 1985) is a corrected-for-chance version of the Rand index, which assigns a score of 0 for the expected Rand index, obtained by random cluster assignment. The maximum value of the Adjusted Rand index is 1. We use this metric only for the toy dataset where we have the ground truth cluster assignments available.

heading: A.10 NEURAL NETWORK HYPERPARAMETERS

heading: B RELATED WORK
text: Traditional survival analysis models. The Cox regression model (Cox, 1992) is a widely used method in survival analysis to estimate the hazard function λ(u)(t) = dF
(u)(t) S(u)(t) , where dF (u) is the
probability density of F (u). The hazard function is then estimated using the covariates, X(u), of a subject u. The hazard function has the form, λ(t|X(u)) = λ0(t) · e{β
TX(u)}, where λ0(t) is a base hazard function common for all subjects, and β are the regression coefficients. The model assumes that the ratio of hazard functions of any two subjects is constant over time. This assumption is violated frequently in real-world datasets (Li et al., 2015a). A near-extreme case when this assumption does not hold is shown in Figure 2b, where the survival curves of two groups of subjects cross each other.
Majority of the work in survival analysis has dealt with the task of predicting the survival outcome especially when the number of features is much higher than the number of subjects (Witten and
Tibshirani, 2010a; pre; Hothorn et al., 2006; Shivaswamy et al., 2007). A number of approaches have also been proposed to perform feature selection in survival data (Ishwaran et al., 2010; Lagani and Tsamardinos, 2010). In the social network scenario, Sun et al. (2012) predicts the relationship building time, that is, the time until a particular link is formed in the network. Alaa and van der Schaar (2017) proposed a nonparametric Bayesian approach for survival analysis in the case of more than one competing events (multiple diseases). They not only assume the presence of termination signals but also the type of event that caused the termination.
Survival (or lifetime) clustering. There have been relatively fewer works that perform lifetime clustering. Many unsupervised approaches have been proposed to identify cancer subtypes in gene expression data without considering the survival outcome (Eisen et al., 1998; Alizadeh et al., 2000). Traditional semi-supervised clustering methods (Aggarwal et al., 2004; Basu et al., 2002; 2004; Nigam et al., 1998) do not perform well in this scenario since they do not provide a way to handle the issues with right censoring. Bair and Tibshirani (2004) proposed a semi-supervised method for clustering survival data in which they assign Cox scores (Cox, 1992) for each feature in their dataset and considered only the features with scores above a predetermined threshold. Then, an unsupervised clustering algorithm, like k-means, is used to group the individuals using only the selected features. Such an approach can miss out on clusters when the features are weakly associated with the survival outcome since such features are discarded immediately after the initial screening. In order to overcome this issue, Gaynor and Bair (2013) proposes supervised sparse clustering as a modification to the sparse clustering algorithm of Witten and Tibshirani (2010b). The sparse clustering algorithm has a modified k-means score that uses distinct weights in the feature set; it initializes these feature weights using Cox scores (Cox, 1992) and optimizes the same objective.
Bair and Tibshirani (2004) and Gaynor and Bair (2013) assume the presence of termination signals. Additionally, there is a disconnect between the use of survival outcomes and the clustering step in these algorithms. In this paper, we provide a loss function that quantifies the divergence between lifetime distributions of the clusters, and we minimize said loss function using a neural network in order to obtain the optimal clusters.
Divergence measures. Since we need to optimize the divergence between two distributions using a finite sample of training data, this problem is related to two-sample tests in statistics. Recently, there is growing interest in loss functions for comparing two distributions, including the Wasserstein distance (Dudley, 2002, p. 420) and the Maximum Mean Discrepancy (MMD) (Gretton et al., 2009; Fukumizu et al., 2008), specifically in the context of training generative adversarial networks (GANs) (Arjovsky et al., 2017; Li et al., 2015b).
Traditional two-sample tests, such as the Kolmogorov-Smirnov (KS) test, are generally disregarded in learning tasks due to their alternating infinite series that makes gradient computations computationally expensive. The alternative state-of-the-art divergence approach that seems most suited to our problem is MMD Gretton et al. (2009), unfortunately, as seeing in our results, they can be less resilient to small sample sizes in empirical distributions.
Deep survival methods. Many deep learning approaches Luck et al. (2017); Katzman et al. (2018); Lee et al. (2018); Ren et al. (2018) have been proposed for predicting the lifetime distribution of a subject given her covariates, while effectively handling censored data that typically arise in survival tasks. DeepHit Lee et al. (2018) introduced a novel architecture and a ranking loss function in addition to the log-likelihood loss for lifetime prediction in the presence of multiple competing risks. Using a log-likelihood loss similar to DeepHit, Ren et al. (2018) propose a recurrent architecture to predict the survival distribution that captures sequential dependencies between neighbouring time points. To the best of our knowledge, there hasn’t been a work on using neural networks for a lifetime clustering task.
Frailty analysis. Extensive research has been done on what is known as frailty analysis, for predicting survival outcomes in the presence of clustered observations (Hougaard, 1995; Chuang et al., 2005; Huang and Wolfe, 2002). Although frailty models provide more flexibility in the presence of clustered observations, they do not provide a mechanism for obtaining the clusters themselves, which is our primary goal.

Review summaries:
This paper proposes a deep learning method for clustering subjects (patients, users of social networks) into clusters of different lifetime (survival time) distributions. Inputs for a subject include covariates (age, gender, etc) and an event sequence. A key difficulty in the problem is that the lifetime of a subject is often unobserved. The proposed method is to learn a termination function that is non-decreasing with time, which essentially treats prolonged inactivity as termination in a probabilistic way. Clustering is done in a discriminative manner where the objective function considers cluster sizes besides between-cluster differences in lifetime distribution. (The inclusion of cluster size information in the objective function avoids degenerate clustering results.)