Title: DEEP LIFETIME CLUSTERING

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.

1. INTRODUCTION:

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.

2. FORMAL FRAMEWORK:

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 = {X k , {(A k,i , M k,i , Y k,i )} i∈N , Θ k } over discrete times t = 0, 1, . . .,

where X k is the random variable representing the covariates of a random subject in cluster k, Θ k is the time to the zeroth activity event (joining), M k,i represent covariates of event i, Y k,i is the inter-event time between the i-th and the (i -1)-st activity events (e.g., logins), and A k,i = 1 indicates an event with a termination signal (death), otherwise A k,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, t m ]. 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 A k,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

T k := max i i ≤i Y k,i i <i (1 -A k,i ) .

The true lifetime T k is unobserved if: (a) the terminal event A k,i = 1 does not occur within the observation time period [0, t m ], 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 S k (t

) := P [T k > t] = 1 -F k (t), t ∈ N ∪ {0}

, where F k (t) is the underlying cumulative distribution function (CDF) of T k . Note again that all the subjects of a cluster share the same lifetime distribution. A i = 1 denotes whether i-th activity event is terminal

Ŝ[t] Ŝ indexed by t ∈ {0, 1, 2 . . .}

Training data. Our training data D consists of subjects from K underlying (hidden) clusters with distinct lifetime distributions S k 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, t m ]. 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) := t m -Θ (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 t m , 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 t m (eg., u 3 ).

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 t m , 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.,

EQUATION

)

where K is a set of all mappings, Ŝk (κ) is the empirical lifetime distribution of subjects in D mapped to cluster k through κ, and ∆ is an empirical distribution divergence measure. u) is the true lifetime of subject u (may be unobserved due to right censoring or due to unobservability of termination signals). H (u) is her observed lifetime, the period between the first and the last observed event. χ (u) is the time between the last observed event and tm, and Q (u) (tm) is the number of events of u after her joining and before tm.

3. THE DEEPCLIFE MODEL:

In this section, we propose a practical lifetime clustering approach using neural networks that optimizes Equation (1). Let D 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, . . . , t max }, where t max = max u∈D H (u) is the maximum observed lifetime of any subject u ∈ D. Note that it is sufficient to define the empirical distribution till t max , since we have not observed any subject with lifetime greater than t max .

3.1 CLUSTER ASSIGNMENTS : α

(u) k (W 1 )

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, α

k (W 1 ) for all k ∈ {1, . . . K},

EQUATION

where W 1 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

EQUATION

} i in order to make it compatible with the feedforward architecture.

3.2 PROBABILITY OF TERMINATION:

β (u) (W 2 )

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 W 2 , we define β (u) (W 2 ) 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) (W 2 ) := A (u) Q (u) (tm) (i.e., W 2 is ignored).

When such termination signals are unobservable (e.g., social network), existing survival methods commonly use a timeout window of predefined size W fixed over the inactive period χ (u) to specify the probability of termination as β (u) (W fixed ) := 1[χ (u) >W fixed ]. However, such specification is arbitrary and precludes any learning of the window size parameter W fixed . Instead, we model the latent termination probabilities β (u) (W 2 ) 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) where W 2 ≥ 0.

β (u) (W 2 ) := 1 -e -W2•χ (

3.3. EMPIRICAL LIFETIME DISTRIBUTION OF CLUSTER:

k : Ŝk (W 1 , W 2 ; D)

Given the training data D and model parameters W 1 and W 2 , we can obtain the soft cluster assignments α (u) k (W 1 ) and the probability of termination β (u) (W 2 ) 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 (W 1 )} u∈D , and the probabilities of termination {β (u) (W 2 )} u∈D , the maximum likelihood estimate of the empirical lifetime distribution of cluster k is given by,

EQUATION

for all t ∈ {0, 1, . . . , t max }, where,

s k (W 1 ; D)[j] = u∈D 1[H (u) ≥ j] • α (u)

k (W 1 ), is the expected number of subjects in cluster k who are at risk (of termination) at time j, and,

d k (W 1 , W 2 ; D)[j] = u∈D 1[H (u) = j] • β (u) (W 2 ) • α (u)

k (W 1 ), 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 : ∆( Ŝa , Ŝb ):

We rewrite the objective function of lifetime clustering from Definition 3 with respect to the model parameters W 1 and W 2 as follows,

EQUATION

where Ŝi is the empirical distribution, a shorthand for the vector Ŝi (W 1 , W 2 ; 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 ∆( Ŝa , Ŝb ):=log(KD( Ŝa , Ŝb )), where KD is the p-value from the Kuiper test between Ŝa and Ŝ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. 

EQUATION

λ a,b = ( M a,b + 0.155 + 0.24 √ M a,b

)V and M a,b = nan b na+n b is the effective sample size. Then, the upper bound3 for the Kuiper p-value is,

EQUATION

)

where v(r, λ)=(4r 2 λ 2 -1)e -2r 2 λ 2 , w(r, λ)=-re -2r 2 λ 2 , r (lo) a,b = 1 √ 2λ a,b

, and, r

(up) a,b = 1 √ 2λ a,b

.

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 10 4 terms (Press et al., 1996) . We optimize Equation (4) with ∆( Ŝa , Ŝb ) := log(KD (up) ( Ŝa , Ŝ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.

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 W 1 and W 2 . We show in the supplementary material that the time and space complexity per iteration of the proposed approach are O(KB + K 2 t max ) and O(Kt max ) 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.

4. RELATED WORK:

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.

5. RESULTS:

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) 

(W fixed ) = 1[χ (u) > W fixed ].

During evaluation, we use the same W fixed 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 W fixed , 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 Ŝ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., Ŝ(u ) := Ŝκ(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 Ŝ(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, t max ]. 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.

5.1. EXPERIMENTS:

Synthetic experiment. We test our method with a synthetic dataset 6 for which we have the true clusters as ground truth. We generate 3 clusters C 1 , C 2 and C 3 with different lifetime distributions such that C 2 and C 3 have proportional hazards, but lifetime distribution of C 1 crosses the other two curves (shown in Figure 2b ). We choose an arbitrary time of measurement, t m =150, to imitate right censoring. We sample 10 4 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: Friendster experiment. Friendster dataset 5 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 W fixed = 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.

EQUATION

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., Ŝ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 dataset 5 (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 .

6. CONCLUSION:

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.

A APPENDIX:

A.1 PROOF OF PROPOSITION 1

Proposition (Proposition 1). Given the training data D, a cluster k, the cluster assignment probabilities {α (u) k (W 1 )} u∈D , and the probabilities of termination {β (u) (W 2 )} u∈D , the maximum likelihood estimate of the empirical lifetime distribution of cluster k is given by,

Ŝk (W 1 , W 2 ; D)[t] = t j=0 s k (W 1 ; D)[j] -d k (W 1 , W 2 ; D)[j] s k (W 1 ; D)[j] ,

for all t ∈ {0, 1, . . . , t max }, where,

s k (W 1 ; D)[j] = u∈D 1[H (u) ≥ j] • α (u)

k (W 1 ), is the expected number of subjects in cluster k who are at risk (of termination) at time j, and,

d k (W 1 , W 2 ; D)[j] = u∈D 1[H (u) = j] • β (u) (W 2 ) • α (u) k (W 1 )

, 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) (W 2 )) 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 (W 1 ), . . . , α (u) K (W 1

)) 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 (T k = t) := λ k,t for t ∈ {0, 1, . . . t max }, where T k 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,

L k = u∈D P (T k = H (u) ) e (u) P (T k > H (u) ) 1-e (u) κ (u) k

.

Next, we take expectation over the cluster assignments and termination signals to get,

EQUATION

)

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,

EQUATION

where, we have used

P (T k = j) = λ k,j , s k [j] is a shorthand for s k (W 1 , W 2 ; D)[j] = u∈D 1[H (u) ≥ j] • α (u)

k (W 1 ), the expected number of subjects at risk at time j, and,

d k [j] is a shorthand for d k (W 1 , W 2 ; D)[j] = u∈D 1[H (u) = j] • β (u) (W 2 ) • α (u)

k (W 1 ), 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, . . . t max }, and computing the lifetime distributions (CCDF) Ŝk from the corresponding lifetime PMFs finishes the proof.

A.2 PROOF OF PROPOSITION 2:

Proposition (Proposition 2). Given two empirical lifetime distributions Ŝa and Ŝb with discrete support and sample sizes n a and n b respectively, define the maximum positive and negative separations between them as, D+ a,b = sup t∈{0,...tmax}

( Ŝa [t] -Ŝb [t]) , D- a,b = sup t∈{0,...tmax} ( Ŝb [t] -Ŝa [t]) .

The Kuiper test p-value Kuiper (1960) gives the probability that Λ, the empirical deviation for n a and n b observations under the null hypothesis S a = S b , exceeds the observed value Then, the lower and upper bounds for the Kuiper p-value are, KD(

EQUATION

EQUATION

)

where v(r, λ) = (4r 2 λ 2 -1)e -2r 2 λ 2 , with w(r, λ) = -re -2r 2 λ 2 , and r

EQUATION

Proof. To prove these bounds, we first regard that the indefinite integral w(r, λ) = v(r, λ)dr = -re -2r 2 λ 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) = 1 √ 2λ , r (up) = 1 √ 2λ and J [l,m] (λ) = m r=l v(r, λ).

Here, we assume that r (lo) = 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,

EQUATION

)

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,

EQUATION

where we assume that r (lo) ≥ 1, otherwise, J [1,r (lo) ] (λ) is trivially zero.

Combining equations equation 11 and equation 12, we get,

EQUATION

)

We use the fact that KD( ) ,∞] (λ), and w(r, λ) = v(r, λ)dr, and that KD( Ŝa , Ŝb ) ≤ 1 to complete the proof.

Ŝa , Ŝb ) = 2 • J [1,∞] (λ a,b ), J [1,∞] (λ) = 1[r (lo) ≥ 1] • J [1,r (lo) ] (λ) + J [r (up

A.3 ADDITIONAL RESULTS:



A.3.1 FRIENDSTER EXPERIMENT:

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. 

A.3.2 MIMIC III EXPERIMENT:

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

A.4 OPTIMAL NUMBER OF CLUSTERS : K:

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 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. where K is a set of all mappings, Ŝ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. divergence between Ŝ1 and Ŝ3 while keeping Ŝ2 = Ŝ1 and Ŝ3 = Ŝ4 . On the other hand, since min i =j ∆( Ŝi , Ŝj ) = 0, the model should maximize the divergence between Ŝ1 and Ŝ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).

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

A.7 TIME/SPACE COMPLEXITY OF DEEPCLIFE-KUIPERUB:

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, . 

D + ← max( Ŝa -Ŝb )

Ŝa , Ŝb are vectors of length t max +1

3: lo) , λ) + V(r (up) , λ) -W(r (up) , λ)

D -← max( Ŝb -Ŝa ) 4: V ← D + + D - Kuiper statistic 5: M ← nan b na+n b Effective sample size 6: λ ← ( √ M + 0.155 + 0.24 √ M )V 7: r (lo) ← 1 √ 2λ 8: r (up) ← 1 √ 2λ 9: if r (lo) ≥ 1 then 10: KD-UB ← W(r (lo) , λ) -W(1, λ) + V(r (

11:

else 12:

KD-UB ← V(r (up) , λ) -W(r (up) , λ)

13:

return log KD-UB 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 C 1 , C 2 and C 3 with different lifetime distributions, as depicted in Figure 13 . Note that while the lifetime distributions C 2 and C 3 follow proportional hazards assumption, the lifetime distribution of C 1 violates this assumption by crossing the other survival curves.

We simulate 10 4 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, t m = 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 µ

EQUATION

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 D C 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.

A.9 METRICS Here we describe the metrics used for evaluating the clusters obtained from the methods. 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. 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.

:

See Appendix A.5 for why we wish to maximize the minimum divergence across all pairs instead of the sum of divergences across all pairs.

The test is typically defined using CDFs but it is equivalent to use lifetime distributions (CCDFs) instead.

Lower bound is presented in the supplementary material.

↑ indicates higher is better, ↓ indicates lower is better. Detailed description in the supplementary material.4 The values in bold are statistically significant with p-value < 0.01 using paired t-test over validation folds.

More details about all the datasets and preprocessing steps can be found in the supplementary material.

i =j ∆( Ŝi , Ŝj ) = 4 • ∆( Ŝ1 , Ŝ3 ). The sum could be further maximized by simply increasing

Figures and Tables:

Figure 1: Depicting two clusters (low-risk and high-risk; as shown by the true lifetime distributions) following different RMPPs, each with two subjects. T (u) is the true lifetime of subject u (may be unobserved due to right censoring or due to unobservability of termination signals). H(u) is her observed lifetime, the period between the first and the last observed event. χ(u) is the time between the last observed event and tm, and Q (u) (tm) is the number of events of u after her joining and before tm.

Figure 2: (a) Feedforward neural network g(; W1) outputs the cluster assignments for a batch of users. The cluster assignments along with the probability of termination is used to obtain lifetime distributions of each cluster using Kaplan-Meier estimator. Finally, logarithm of Kuiper p-value upper bound is used as the divergence loss ∆. (b) Lifetime distributions can have different shapes and can cross each other, violating proportional hazards assumptions. (c) Divergence metric must account for the uncertainty in the distributions, otherwise divergence maximization leads to imbalanced clusters. (d-e) Upper and lower bounds of the logarithm of Kuiper p-value when varying the Kuiper statistic D + +D -.

(Bounds for Kuiper p-value.) Given two empirical lifetime distributions Ŝa and Ŝb with discrete support and sample sizes n a and n b respectively, define the maximum positive and negative separations between them, t] -Ŝa [t]) . The Kuiper test p-value Kuiper (1960) gives the probability that Λ, the empirical deviation for n a and n b observations under the null hypothesis S a = S b 2 , exceeds the observed value V = D+ a,b + Da,b :

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) DeepCLife-KuiperUB, 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)

Figure 3: (Friendster) Empirical lifetime distributions of clusters obtained from different methods for K=3 (legend shows cluster sizes n1, n2, n3). Baseline methods (a-c) employ a two-stage clustering process and do not guarantee clusters with maximally different lifetime distributions. (d) DeepCLife-MMD suffers from sample anomalies (n2 = 0). (e) DeepCLife-KuiperUB obtains clusters with significantly different lifetime distributions (best Logrank scores). 0 5 10 15 20 25 30 Time (in days) 0.0 0.2 0.4 0.6 0.8 1.0 Probability (CCDF) S 1 (n 1 = 2957) S 2 (n 2 = 3955)

where λ a,b = M a,b + 0.155 + 0.24 √ M a,b V and M a,b = nan b na+n b is the effective sample size.

Figure 5: (Friendster) Empirical lifetime distributions of clusters obtained from different methods for K=2 (legend shows cluster sizes n1, n2).

Figure 6: (Friendster) Empirical lifetime distributions of clusters obtained from different methods for K=3 (legend shows cluster sizes n1, n2, n3).

Figure 8: (Friendster) Empirical lifetime distributions of clusters obtained from different methods for K=5 (legend shows cluster sizes).

CHOICE OF MINIMUM VS SUM IN EQUATIONS (1) AND (4) We define the clustering problem in Definition 3 as finding a mapping κ such that, κ = arg max κ∈K min i,j∈{1...K}, i =j ∆( Ŝi (κ), Ŝj (κ)) ,

Figure 9: (MIMIC III) Empirical lifetime distributions of clusters obtained from DeepCLife-KuiperUB with K=2 . . . 6 (legend shows cluster sizes).

Figure 10: (Friendster) Logrank scores (error bars denote standard errors) on a held-out validation set for the clusters obtained from DeepCLife-KuiperUB with different values of K. K = 4 seems to give the best mean performance, while K = 5 is close behind; increasing K further drastically reduces the Logrank score.

Figure 13: (Toy) True lifetime distributions of simulated clusters (with t m = 150).

Deep survival methods. Many deep learning approaches Luck et al. (2017); Katzman et al. (2018); Lee et al. (2018); Ren et al. (

Table of notations.

Table content: <table><tr><td colspan="2">Notation Description</td><td colspan="2">Notation Description</td></tr><tr><td>(u)</td><td>Variables indexed by a subject u</td><td>T</td><td>True lifetime of the subject (unobserved)</td></tr><tr><td>k</td><td>Random variables indexed by a random</td><td>H</td><td>Observed lifetime of the subject</td></tr><tr><td/><td>subject of cluster k</td><td/><td/></tr><tr><td>Θ</td><td>Joining time of the subject</td><td>χ</td><td>Time elapsed since the last observed ac-</td></tr><tr><td/><td/><td/><td>tivity event before t m</td></tr><tr><td>X</td><td>Covariates of the subject</td><td>F</td><td>CDF of the true lifetime T</td></tr><tr><td>Y i</td><td>Inter-event time between i-th and (i-1)-</td><td>S</td><td>True lifetime distribution (≡ CCDF of</td></tr><tr><td/><td>st activity events</td><td/><td>true lifetime T )</td></tr><tr><td>M i</td><td>Covariates of event i</td><td>Ŝ</td><td>Empirical lifetime distribution</td></tr><tr><td>A i</td><td/><td/><td/></tr></table>

(Synthetic) C-index (%) and Adjusted Rand index (%) for clusters with standard errors in parentheses for different methods 5 .

Table content: <table><tr><td/><td colspan="2">D {C1,C2}</td><td colspan="2">D {C1,C3}</td><td colspan="2">D {C1,C2,C3}</td></tr><tr><td>Method</td><td>C-index ↑ (%)</td><td>Adj. Rand Index ↑ (%)</td><td>C-index ↑ (%)</td><td>Adj. Rand Index ↑ (%)</td><td>C-index ↑ (%)</td><td>Adj. Rand Index ↑ (%)</td></tr><tr><td>SSC-Bair</td><td>62.75 (0.35)</td><td>74.66 (0.48)</td><td>62.99 (0.26)</td><td>56.86 (0.63)</td><td>63.77 (0.24)</td><td>47.67 (0.24)</td></tr><tr><td>SSC-Gaynor</td><td>56.34 (0.50)</td><td>19.88 (0.51)</td><td>57.21 (0.43)</td><td>16.60 (0.56)</td><td>56.75 (0.32)</td><td>5.84 (0.11)</td></tr><tr><td>DeepHit+GMM</td><td>63.31 (0.32)</td><td>85.05 (1.01)</td><td>65.23 (0.21)</td><td>78.47 (0.94)</td><td>59.59 (1.98)</td><td>38.77 (7.16)</td></tr><tr><td>DeepCLife-MMD</td><td>64.35 (0.33)</td><td>98.47 (0.28)</td><td>67.25 (0.26)</td><td>99.68 (0.16)</td><td>62.17 (0.71)</td><td>36.75 (1.10)</td></tr><tr><td>DeepCLife-KuiperUB</td><td>64.37 (0.32)</td><td>99.02 (0.14)</td><td>67.24 (0.26)</td><td>99.94 (0.06)</td><td>68.96 (0.38)</td><td>73.61 (0.62)</td></tr></table>

(Friendster) C-index (%), Integrated Brier Score (%) and Logrank score with standard errors in parentheses for different methods 4 and K=2, 4 clusters with number of training examples N (tr) =10 5 .

Table content: <table><tr><td/><td/><td>K = 2</td><td/><td>K = 4</td><td/></tr><tr><td>Method</td><td>C-index ↑ (%)</td><td>I.B.S ↓ (%)</td><td>Logrank Score ↑</td><td>C-index ↑ (%) I.B.S. ↓ (%)</td><td>Logrank Score ↑</td></tr><tr><td>SSC-Bair</td><td colspan="2">64.42 (0.15) 22.18 (0.02)</td><td>5479.27 ( 38.06)</td><td colspan="2">67.18 (0.13) 21.55 (0.03) 13013.86 ( 182.07)</td></tr><tr><td>SSC-Gaynor</td><td colspan="2">64.42 (0.18) 22.17 (0.02)</td><td>5557.41 ( 38.43)</td><td colspan="2">69.99 (0.28) 21.62 (0.01) 15204.81 ( 41.29)</td></tr><tr><td>DeepHit+GMM</td><td colspan="2">64.33 (1.80) 22.04 (0.23)</td><td>9207.46 (5031.48)</td><td colspan="2">76.55 (0.12) 20.64 (0.02) 40703.01 ( 477.25)</td></tr><tr><td>DeepCLife-MMD</td><td colspan="3">67.49 (0.11) 22.07 (0.02) 27642.60 (1301.85)</td><td colspan="2">70.93 (1.80) 22.40 (0.07) 33030.93 (4277.24)</td></tr><tr><td>DeepCLife-KuiperUB</td><td colspan="3">75.58 (0.15) 20.13 (0.02) 47837.25 ( 297.63)</td><td colspan="2">77.04 (0.88) 18.99 (0.20) 59236.36 (2126.55)</td></tr></table>

(Friendster) C-index (%), Integrated Brier Score (%) and Logrank score with standard errors in parentheses for different methods 3 and K=3, 5 clusters with number of training examples N (tr) =10 5 .

Table content: <table><tr><td/><td/><td>K = 3</td><td/><td>K = 5</td></tr><tr><td>Method</td><td>C-index ↑ (%)</td><td>I.B.S ↓ (%)</td><td colspan="2">Logrank Score ↑ C-index ↑ (%) I.B.S. ↓ (%)</td><td>Logrank Score ↑</td></tr><tr><td>SSC-Bair</td><td colspan="2">64.23 (0.24) 22.19 (0.02)</td><td>5277.82 ( 43.63)</td><td colspan="2">75.80 (0.16) 20.26 (0.02) 43605.16 ( 91.39)</td></tr><tr><td>SSC-Gaynor</td><td colspan="2">64.13 (0.21) 22.18 (0.02)</td><td>5400.75 ( 39.59)</td><td colspan="2">75.66 (0.20) 20.48 (0.01) 41824.28 ( 151.85)</td></tr><tr><td>DeepHit+GMM</td><td colspan="3">74.81 (0.11) 20.97 (0.06) 33880.55 ( 426.74)</td><td colspan="2">76.17 (0.27) 20.54 (0.11) 39254.67 (1721.56)</td></tr><tr><td>DeepCLife-MMD</td><td colspan="3">72.75 (1.24) 18.94 (0.65) 50173.94 (4809.80)</td><td colspan="2">63.08 (2.69) 21.55 (0.54) 14308.90 (9149.41)</td></tr><tr><td>DeepCLife-KuiperUB</td><td colspan="3">78.48 (0.32) 19.70 (0.55) 54191.28 ( 436.72)</td><td colspan="2">76.51 (1.25) 18.99 (0.16) 56130.21 (4036.02)</td></tr></table>

. . t max . First, we can compute the empirical lifetime distributions of all the K clusters in O(KB) time using O(Kt max ) space. Computing upper bound of the Kuiper p-value for two lifetime distributions using Algorithm 1 requires O(t max ) time and O(t max ) space. Since we compute the upper bound for all K 2 pairs, the time complexity then becomes O(K 2 t max ) whereas the space complexity remains O(t max ). The total time and space complexity of the proposed algorithm are O(KB + K 2 t max ) and O(Kt max ) respectively. Require: Two empirical lifetime distributions Ŝa and Ŝb with sample sizes n a and n b respectively 1: function KUIPERUB( Ŝa , Ŝb , n a , n b )

Table content: <table><tr><td>Algorithm 1 Kuiper Divergence Upper Bound (Log)</td></tr></table>

Peer review (iclr2020):
The authors propose an approach to cluster subjects into K clusters according to their underlying, often unknown, life distribution. The loss function to the neural-network-based approach is based on a tight upper bound to the two-sample Kuiper test. Experiments are conducted on artificial and two real life datasets. I enjoyed reading the paper, the authors propose an approach to address an important yet relatively under-explored problem in survival analysis. It is not entirely clear how to handle the case when after the model is trained, H^(u) for a new subject u is larger than t_max when the model is trained. In such case, Chi^(u) will be negative, thus what happens with _x0008_eta^(u)(W_2)? When the termination signals are not available (friendster data) termination signals are set artificially when training and evaluation, thus, is the model learning artificially set termination signals and thus artificial survival functions? what is the point of doing this? Moreover, if the termination signals are artificial, what is the point of calculating C-Index and Integrated Brier score? I understand the motivation to compare to existing methods, it is really not clear how meaningful they are at evaluating performance. In general, using C-Index and Brier scores in the context of the presented application may be misleading (or non-applicable) because these metrics are usually applied to survival analysis scenarios where one seeks to estimate likelihood of survival or time to event (over a usually infinite time horizon). Here, _x0008_eta^(u)(W_2) is a function of Chi^(u) which is known (also not comparable to survival analysis), artificially tied to a pre-specified time-horizon t_max, and not dependent on covariates as in survival analysis. Further, SSC-Bair and SSC-Gaynor are clustering models, how are survival estimates obtained for these? In practice, how is K selected? based on the friendster results for K=2,3,4,5 they all seem to produce distinct clusters, so which one should be used? This raises a question: how are the clusters, their members or the number of clusters informing the use case? In the MIMIC III experiment, where events are observed, the model finds two clusters without strikingly different survival functions (relative to friendster). Can one really consider S_2 as the high-risk group? Can one get higher risk clusters with larger K? Minor: - A_i^(u) is missing from the definition in the training data. - In Section 3.2, if xi^(u) is shared, why the superscript indicating that is subject specific?