Year: 2019

Abstract:
Due to the wider availability of modern electronic health records (EHR), patient care data is often being stored in the form of time-series. Clustering such timeseries data is crucial for patient phenotyping, anticipating patients’ prognoses by identifying “similar” patients, and designing treatment guidelines that are tailored to homogeneous patient subgroups. In this paper, we develop a deep learning approach for clustering time-series data, where each cluster comprises patients who share similar future outcomes of interest (e.g., adverse events, the onset of comorbidities, etc.). The clustering is carried out by using our novel loss functions that encourage each cluster to have homogeneous future outcomes. We adopt actorcritic models to allow “back-propagation” through the sampling process that is required for assigning clusters to time-series inputs. Experiments on two realworld datasets show that our model achieves superior clustering performance over state-of-the-art benchmarks and identifies meaningful clusters that can be translated into actionable information for clinical decision-making.

text: Due to the wider availability of modern electronic health records (EHR), patient care data is often being stored in the form of time-series. Clustering such timeseries data is crucial for patient phenotyping, anticipating patients’ prognoses by identifying “similar” patients, and designing treatment guidelines that are tailored to homogeneous patient subgroups. In this paper, we develop a deep learning approach for clustering time-series data, where each cluster comprises patients who share similar future outcomes of interest (e.g., adverse events, the onset of comorbidities, etc.). The clustering is carried out by using our novel loss functions that encourage each cluster to have homogeneous future outcomes. We adopt actorcritic models to allow “back-propagation” through the sampling process that is required for assigning clusters to time-series inputs. Experiments on two realworld datasets show that our model achieves superior clustering performance over state-of-the-art benchmarks and identifies meaningful clusters that can be translated into actionable information for clinical decision-making.

heading: 1 INTRODUCTION
text: Chronic diseases – such as cystic fibrosis, dementia, and diabetes – are heterogeneous in nature, with widely differing outcomes even in narrow patient subgroups. Disease progression manifests through a broad spectrum of clinical factors, collected as a sequence of measurements over time in electronic health records (EHR), which gives a rise to complex progression patterns among patients (Samal et al., 2011). For example, cystic fibrosis evolves slowly, allowing for the development of related comorbidities and bacterial infections, and creating distinct behaviors/responses to therapeutic interventions, which in turn makes the survival and quality of life substantially different (Ramos et al., 2017). Identifying patient subgroups with similar progression patterns can be advantageous for understanding such heterogeneous underlying diseases. This allows clinicians to anticipate patients’ prognoses by comparing “similar” patients in order to design treatment guidelines that are tailored to homogeneous patient subgroups (Zhang et al., 2019).
Temporal clustering has been recently used as a data-driven framework to partition patients with time-series observations into a set of clusters (i.e., into subgroups of patients). Recent research has typically focused on either finding fixed-length and low-dimensional representations (Zhang et al., 2019; Rusanov et al., 2016) or on modifying the similarity measure (Giannoula et al., 2018; Luong and Chandola, 2017) both in an attempt to apply the conventional clustering algorithms (e.g., K-means (Lloyd, 1982)) to time-series observations. However, clusters identified from these approaches these approaches are purely unsupervised – they do not account for each patient’s observed outcome (e.g., adverse events, the onset of comorbidities, etc.) – which leads to heterogeneous clusters if the clinical presentation of the disease differs even for similar patients. Thus, a common prognosis in each cluster remains unknown which can mystify the understanding of the underlying disease progression (Boudier et al., 2019). For instance, patients who appear to have similar timeseries observations may develop different sets of comorbidities in the future which, in turn, require different treatment guidelines to reduce such risks (Wami et al., 2013). To overcome this limitation, we focus on predictive clustering (Blockeel et al., 2017) which combines prediction with clustering. Therefore, the cluster assignments are optimized such that patients in a cluster share similar future outcomes to provide a prognostic value.
In this paper, we propose an actor-critic approach for temporal predictive clustering, which we call AC-TPC.1 Our model consists of three neural networks – an encoder, a selector, and a predictor – and a set of centroid candidates. More specifically, the encoder maps an input time-series into a latent encoding; the selector utilizes the encoding and assigns a cluster to which the time-series belongs to via a sampling process; and the predictor estimates the future outcome distribution conditioned on either the encoding or the centroid of the selected cluster. The following three contributions render our model to identify the predictive clusters. First, to encourage each cluster to have homogeneous future outcomes, we define a clustering objective based on the Kullback-Leibler (KL) divergence between the predictor’s output given the input time series, and the predictor’s output given estimated cluster assignments. Second, we transform solving a combinatorial problem of identifying cluster into iteratively solving two sub-problems: optimization of the cluster assignments and optimization of the cluster centroids. Finally, we allow “back-propagation” through the sampling process of the selector by adopting the training of actor-critic models (Konda and Tsitsiklis, 2000).
Throughout the experiments, we show significant performance improvements over the state-of-theart clustering methods on two real-world medical datasets. Then, to demonstrate the practical significance of our model, we consider a more realistic scenario where the future outcomes of interest are high-dimensional – such as, development of multiple comorbidities in the next year – and interpreting all possible combinations is intractable. Our experiments show that our model can identify meaningful clusters that can be translated into actionable information for clinical decision-making.

heading: 2 PROBLEM FORMULATION
text: Let X ∈ X and Y ∈ Y be random variables for an input feature and an output label (i.e., one or a combination of future outcome(s) of interest) with a joint distribution pXY (and marginal distributions are pX and pY , respectively) where X is the feature space and Y is the label space. Here, we focus our description on C-class classification tasks, i.e., Y = {1, · · · , C}.2 We are given a time-series dataset D = {(xnt , ynt )T n
t=1}Nn=1 comprising sequences of realizations (i.e., observations) of the pair (X, Y ) for N patients. Here, (xnt , y n t ) Tn t=1 is a sequence of T n observation pairs that correspond to patient n and t ∈ T n , {1, · · · , Tn} denotes the time stamp at which the observations are made.From this point forward, we omit the dependency on n when it is clear in the context and denote x1:t = (x1, · · · ,xt) for ease of notation. Our aim is to identify a set of K predictive clusters, C = {C(1), · · · , C(K)}, for time-series data. Each cluster consists of homogeneous data samples, that can be represented by its centroid, based on a certain similarity measure. There are two main distinctions from the conventional notion of clustering. First, we treat subsequences of each times-series as data samples and focus on partitioning {{xn1:t}T n
t=1}Nn=1 into C. Hence, we define a cluster as C(k) = {xn1:t|t ∈ T n, snt = k} for k ∈ K , {1, · · · ,K} where snt ∈ K is the cluster assignment for a given xn1:t. This is to flexibly update the cluster assignment (in real-time) to which a patient belongs as new observations are being accrued over time. Second, we define the similarity measure with respect to the label distribution and associate it with clusters to provide a prognostic value. More specifically, we want the distribution of output label for subsequences in each cluster to be homogeneous and, thus, can be well-represented by the centroid of that cluster. Let S be a random variable for the cluster assignment – that depends on a given subsequence x1:t – and Y |S = k be a random variable for the output given cluster k. Then, such property of predictive clustering can be achieved by minimizing the following Kullback-Leibler (KL) divergence:
KL ( Yt|X1:t = x1:t ∥∥Yt|St = k) for x1:t ∈ C(k) (1) where KL ( Yt|X1:t = x1:t
∥∥Yt|St = k) = ∫y∈Y p(y|x1:t)( log p(y|x1:t) − log p(y|st))dy. Here, p(y|x1:t) and p(y|st) are the label distributions conditioned on a subsequence x1:t and a cluster assignment st, respectively. Note that (1) achieves its minimum when the two distributions are equivalent.
1Source code available at https://github.com/ICLR2020-ACTPC/ACTPC_submission.git 2In this paper, we focus our description on C-class classification task, i.e., Y = {1, · · · , C}; in Appendix A, we discuss simple modifications of our model for regression and M -dimensional binary classification tasks, i.e., Y = R and Y = {0, 1}M , respectively.
Finally, we establish our goal as identifying a set of predictive clusters C that optimizes the following objective:
minimize C ∑ k∈K ∑ x1:t∈C(k) KL ( Yt|X1:t = x1:t ∥∥Yt|St = k). (2) Unfortunately, the optimization problem in (2) is highly non-trivial. We need to estimate the objective function in (2) while solving a non-convex combinatorial problem of finding the optimal cluster assignments and cluster centroids.

heading: 3 ACTOR-CRITIC APPROACH FOR TEMPORAL PREDICTIVE CLUSTERING
text: To effectively estimate the objective function in (2), we introduce three networks – an encoder, a selector, and a predictor – and an embedding dictionary as illustrated in Figure 1. These components together provide the cluster assignment and the corresponding centroid based on a given sequence of observations and enable us to estimate the probability density p(y|st). More specifically, we define each component as follows:
• The encoder, fθ : ∏t i=1 X → Z , is a RNN (parameterized by θ) that maps a (sub)sequence of
a time-series x1:t to a latent representation (i.e., encoding) zt ∈ Z where Z is the latent space. • The selector, hψ : Z → ∆K−1, is a fully-connected network (parameterized by ψ) that pro-
vides a probabilistic mapping to a categorical distribution from which the cluster assignment st ∈ K is being sampled. • The predictor, gφ : Z → ∆C−1, is a fully-connected network (parameterized by φ) that esti-
mates the label distribution given the encoding of a time-series or the centroid of a cluster.
• The embedding dictionary, E = {e(1), · · · , e(K)} where e(k) ∈ Z for k ∈ K, is a set of cluster centroids lying in the latent space which represents the corresponding cluster.
Here, ∆D−1 = {q ∈ [0, 1]K : q1 + · · ·+ qD = 1} is a (D− 1)-simplex that denotes the probability distribution for a D-dimensional categorical (class) variable.
At each time stamp t, the encoder maps a input (sub)sequence x1:t into a latent encoding zt , fθ(x1:t). Then, based on the encoding zt, the cluster assignment st is drawn from a categorical distribution that is defined by the selector output, i.e., st ∼ Cat(πt) where πt = [πt(1), · · · , πt(K)] , hψ(zt). Once the assignment st is chosen, we allocate the latent encoding zt to an embedding e(st) in the embedding dictionary E . Since the allocated embedding e(st) corresponds to the centroid of the cluster to which x1:t belongs, we can, finally, estimate the density p(y|st) in (2) as the output of the predictor given the embedding e(st), i.e., ȳt , gφ(e(st)).

heading: 3.1 LOSS FUNCTIONS
text: In this subsection, we define loss functions to achieve our objective in (2); the details of how we train our model will be discussed in the following subsection.
Predictive Clustering Loss: Since finding the cluster assignment of a given sequence is a probabilistic problem due to the sampling process, the objective function in (2) must be defined as an expectation over the cluster assignment. Thus, we can estimate solving the objective problem in (2) as minimizing the following loss function:
L1(θ, ψ, φ, E) = Ex,y∼pXY [ T∑ t=1 Est∼Cat(πt) [ `1(yt, ȳt) ]] (3)
where `1(yt, ȳt) = − ∑C c=1 y c t log ȳ c t . Here, we slightly abuse the notation and denote y = [y1 · · · yC ] as the one-hot encoding of y, and yc and ȳc indicates the c-th component of y and ȳ, respectively. It is worth to highlight that minimizing `1 is equivalent to minimizing the KL divergence in (2) since the former term of the KL divergence is independent of our optimization procedure.
One critical question that may arise is how to avoid trivial solutions in this unsupervised setting of identifying the cluster assignments and the centroids (Yang et al., 2017). For example, all the embeddings in E may collapse into a single point or the selector simply assigns equal probability to all the clusters regardless of the input sequence. In both cases, our model will fail to correctly estimate p(y|st) and, thus, end up finding a trivial solution. To address this issue, we introduce two auxiliary loss functions that are tailored to address this concern. It is worth to highlight that these loss functions are not subject to the sampling process and their gradients can be simply back-propagated.
Sample-Wise Entropy of Cluster Assignment: To motivate sparse cluster assignment such that the selector ultimately selects one dominant cluster for each sequence, we introduce sample-wise entropy of cluster assignment which is given as
L2(θ, ψ) = Ex∼pX
[ −
T∑ t=1 ∑ k∈K πt(k) log πt(k)
] (4)
where πt = [πt(1) · · ·πt(K)] = hψ(fθ(x1:t)). The sample-wise entropy achieves its minimum when πt becomes an one-hot vector.
Embedding Separation Loss: To prevent the embeddings in E from collapsing into a single point, we define a loss function that encourages the embeddings to represent different label distributions, i.e., gφ(e(k)) for k ∈ K, from each other:
L3(E) = − ∑ k 6=k′ `1(gφ(e(k)), gφ(e(k ′))) (5)
where `1 is reused to quantify the distance between label distributions conditioned on each cluster. We minimize (5) when updating the embedding vectors e(1), · · · , e(K).

heading: 3.2 OPTIMIZATION
text: The optimization problem in (2) is a non-convex combinatorial problem because it comprises not only minimizing the KL divergence but also finding the optimal cluster assignments and centroids. Hence, we propose an optimization procedure that iteratively solves two subproblems: i) optimizing the three networks – the encoder, selector, and predictor – while fixing the embedding dictionary and ii) optimizing the embedding dictionary while fixing the three networks. Pseudo-code of ACTPC can be found in Appendix F.
3.2.1 OPTIMIZING THE THREE NETWORKS – fθ , hψ , AND gφ
Finding predictive clusters incorporates the sampling process which is non-differentiable. Thus, to render “back-propagation”, we utilize the training of actor-critic models (Konda and Tsitsiklis, 2000). More specifically, we view the combination of the encoder (fθ) and the selector (hψ) as the “actor” parameterized by ωA = [θ, ψ], and the predictor (gφ) as the “critic”. The critic takes as input the the output of the actor (i.e., the cluster assignment) and estimates its value based on the samplewise predictive clustering loss (i.e., `1(yt, ȳt)) given the chosen cluster. This, in turn, renders the actor to change the distribution of selecting a cluster to minimize such loss. Thus, it is important for the critic to perform well on the updated output of the actor while it is important for the actor
to perform well on the updated loss estimation. As such, the parameters for the actor and the critic need to be updated iteratively.
Given the embedding dictionary E fixed (thus, we will omit the dependency on E), we train the actor, i.e., the encoder and the selector, by minimizing a combination of the predictive clustering loss L1 and the entropy of cluster assignments L2, which is given by
LA(θ, ψ, φ) = L1(θ, ψ, φ) + αL2(θ, ψ) (6)
where α ≥ 0 is a coefficient chosen to balance between the two losses. To derive the gradient of this loss with respect ωA = [θ, ψ], we utilize the ideas from actor-critic models (Konda and Tsitsiklis, 2000) as follows; please refer to Appendix B for the detailed derivation:
∇ωALA(θ, ψ, φ) = Ex,y∼pXY [ ∇ωA ( T∑ t=1 Est∼Cat(πt) [ `1(yt, ȳt) ])] + α∇ωAL2(θ, ψ)
= Ex,y∼pXY [ T∑ t=1 Est∼Cat(πt) [ `1(yt, ȳt)∇ωA log πt(st) ]] + α∇ωAL2(θ, ψ). (7)
Note that since no sampling process is considered in L2(θ, ψ), we can simply derive∇ωAL2(θ, ψ). Iteratively with training the actor, we train the critic, i.e., the predictor, by minimizing the predictive clustering loss L1 as the following:
LC(φ) = L1(θ, ψ, φ) (8)
whose gradient with respect to φ can be givens as ∇φLC(φ) = ∇φL1(θ, ψ, φ). Note that since the critic is independent of the sampling process, the gradient can be simply back-propagated.

heading: 3.2.2 OPTIMIZING THE CLUSTER CENTROIDS
text: Now, once the parameters for the three networks (θ, ψ, φ) are fixed (thus, we omit the dependency on θ, ψ, and φ), we updated the embeddings in E by minimizing a combination of the predictive clustering loss L1 and the embedding separation loss L3, which is given by
LE(E) = L1(E) + βL3(E) (9)
where β ≥ 0 is a coefficient chosen to balance between the two losses.

heading: 3.2.3 INITIALIZING AC-TPC VIA PRE-TRAINING
text: Since we transform the non-trivial combinatorial optimization problem in (2) into iteratively solving two sub-problems, initialization is crucial to achieve better optimization as a similar concern has been addressed in (Yang et al., 2017).
Therefore, we initialize our model based on the following procedure. First, we pre-train the encoder and the predictor by minimizing the following loss function based on the predicted label distribution given the latent encodings of input sequences, i.e., ŷt , gφ(zt) = gφ(fθ(x1:t)), as the following:
LI(θ, φ) = Ex,y∼pXY
[ −
T∑ t=1 `1(yt, ŷt)
] . (10)
Minimizing (10) encourages the latent encoding to be enriched with information for accurately predicting the label distribution. Then, we perform K-means (other clustering method can be also applied) based on the learned representations to initialize the embeddings E and the cluster assignments {{snt }T n
t=1}Nn=1. Finally, we pre-train the selector hψ by minimizing the cross entropy treating the initialized cluster assignments as the true clusters.

heading: 4 RELATED WORK
text: Temporal clustering, also known as time-series clustering, is a process of unsupervised partitioning of the time-series data into clusters in such a way that homogeneous time-series are grouped together
based on a certain similarity measure. Temporal clustering is challenging because i) the data is often high-dimensional – it consists of sequences not only with high-dimensional features but also with many time points – and ii) defining a proper similarity measure for time-series is not straightforward since it is often highly sensitive to distortions (Ratanamahatana et al., 2005). To address these challenges, there have been various attempts to find a good representation with reduced dimensionality or to define a proper similarity measure for times-series (Aghabozorgi et al., 2015).
Recently, Baytas et al. (2017) and Madiraju et al. (2018) proposed temporal clustering methods that utilize low-dimensional representations learned by RNNs. These works are motivated by the success of applying deep neural networks to find “clustering friendly” latent representations for clustering static data (Xie et al., 2017; Yang et al., 2017). In particular, Baytas et al. (2017) utilized a modified LSTM auto-encoder to find the latent representations that are effective to summarize the input time-series and conductedK-means on top of the learned representations as an ad-hoc process. Similarly, Madiraju et al. (2018) proposed a bidirectional-LSTM auto-encoder that jointly optimizes the reconstruction loss for dimensionality reduction and the clustering objective. However, these methods do not associate a target property with clusters and, thus, provide little prognostic value in understanding the underlying disease progression.
Our work is most closely related to SOM-VAE (Fortuin et al., 2019). This method jointly optimizes a static variational auto-encoder (VAE), that finds latent representations of input features, and a self-organizing map (SOM), that allows to map the latent representations into a more interpretable discrete representations, i.e., the embeddings. However, there are three key differences between our work and SOM-VAE. First, SOM-VAE aims at minimizing the reconstruction loss that is specified as the mean squared error between the original input and the reconstructed input based on the corresponding embedding. Thus, similar to the aforementioned methods, SOM-VAE neither associates future outcomes of interest with clusters. In contrast, we focus on minimizing the KL divergence between the outcome distribution given the original input sequence and that given the corresponding embedding to build association between future outcomes of interest and clusters. Second, to overcome non-differentiability caused by the sampling process (that is, mapping the latent representation to the embeddings), Fortuin et al. (2019) applies the gradient copying technique proposed by (van den Oord et al., 2017), while we utilize the training of actor-critic model (Konda and Tsitsiklis, 2000). Finally, while we flexibly model time-series using LSTM, SOM-VAE handles time-series by integrating a Markov model in the latent representations. This can be a strict assumption especially in clinical settings where a patient’s medical history is informative for predicting his/her future clinical outcomes (Ranganath et al., 2016).

heading: 5 EXPERIMENTS
text: In this section, we provide a set of experiments using two real-world time-series datasets. We iteratively update the three networks – the encoder, selector, and predictor – and the embedding dictionary as described in Section 3.2. For the network architecture, we constructed the encoder utilizing a single-layer LSTM (Hochreiter and Schmidhuber, 1997) with 50 nodes and constructed the selector and predictor utilizing two-layer fully-connected network with 50 nodes in each layer, respectively. The parameters (θ, ψ, φ) are initialized by Xavier initialization (Glorot and Bengio, 2010) and optimized via Adam optimizer (Kingma and Ba, 2014) with learning rate of 0.001 and keep probability 0.7. We chose the balancing coefficients α, β ∈ {0.1, 1.0, 3.0} utilizing grid search that achieves the minimum validation loss in (3); the effect of different loss functions are further investigated in the experiments. Here, all the results are reported using 5 random 64/16/20 train/validation/test splits.

heading: 5.1 REAL-WORLD DATASETS
text: We conducted experiments to investigate the performance of AC-TPC on two real-world medical datasets; detailed statistics of each dataset can be found in Appendix C:
• UK Cystic Fibrosis registry (UKCF)3: This dataset records annual follow-ups for 5,171 adult patients (aged 18 years or older) enrolled in the UK CF registry over the period from 2008 and 2015, with a total of 25,012 hospital visits. Each patient is associated with 89 variables (i.e., 11 static and 78 time-varying features), including information on demographics and genetic
3https://www.cysticfibrosis.org.uk/the-work-we-do/uk-cf-registry
mutations, bacterial infections, lung function scores, therapeutic managements, and diagnosis on comorbidities. We set the development of different comorbidities in the next year as the label of interest at each time stamp.
• Alzheimer’s Disease Neuroimaging Initiative (ADNI)4: This dataset consists of 1,346 patients in the Alzheimer’s disease study with a total of 11,651 hospital visits, which tracks the disease progression via follow-up observations at 6 months interval. Each patient is associated with 21 variables (i.e., 5 static and 16 time-varying features), including information on demographics, biomarkers on brain functions, and cognitive test results. We set predictions on the three diagnostic groups – normal brain functioning, mild cognitive impairment, and Alzheimer’s disease – as the label of interest at each time stamp.

heading: 5.2 BENCHMARKS
text: We compare AC-TPC with clustering methods ranging from conventional approaches based on Kmeans to the state-of-the-art approaches based on deep neural networks. All the benchmarks compared in the experiments are tailored to incorporate time-series data as described below:
• Dynamic time warping followed by K-means: Dynamic time warping (DTW) is utilized to quantify pairwise distance between two variable-length sequences and, then,K-means is applied (denoted as KM-DTW).
• K-means with deep neural networks: To handle variable-length time-series data, we utilize our encoder and predictor that are trained based on (10) for dimensionality reduction; this is to provide fixed-length and low-dimensional representations for time-series. Then, we apply K-means on the latent encodings z (denoted as KM-E2P (Z)) and on the predicted label distributions ŷ (denoted as KM-E2P (Y)), respectively.
• Extensions of DCN (Yang et al., 2017): Since the DCN is designed for static data, we replace their static auto-encoder with a sequence-to-sequence network to incorporate time-series data (denoted as DCN-S2S).5 In addition, to associated with the label distribution, we compare a DCN whose static auto-encoder is replaced with our encoder and predictor (denoted as DCNE2P) to focus dimensionality reduction while preserving information for predicting the label.
• SOM-VAE (Fortuin et al., 2019): We compare with SOM-VAE – though, this method is oriented towards visualizing input data via SOM – since it naturally clusters time-series data (denoted as SOM-VAE). In addition, we compare with a variation of SOM-VAE by replacing the decoder with our predictor in order to find embeddings that capture information for predicting the label (denoted as SOM-VAE-P). For both cases, we set the dimension of SOM to K.
It is worth highlighting that the label information is provided for training DCN-E2P, KM-E2P, and SOM-VAE-P while the label information is not provided for training KM-DTW, DCN-S2S, and SOM-VAE. Please refer to Appendix D for the summary of major components of the tested benchmarks and the implementation details.

heading: 5.3 PERFORMANCE METRICS
text: To assess the clustering performance, we applied the following three standard metrics for evaluating clustering performances when the ground-truth cluster label is available: purity score, normalized mutual information (NMI) (Vinh et al., 2010), and adjusted Rand index (ARI) (Hubert and Arabie, 1985). More specifically, the purity score assesses how homogeneous each cluster is (ranges from 0 to 1 where 1 being a cluster consists of a single class), the NMI is an information theoretic measure of how much information is shared between the clusters and the labels that is adjusted for the number of clusters (ranges from 0 to 1 where 1 being a perfect clustering), and ARI is a corrected-for-chance version of the Rand index which is a measure of the percentage of correct cluster assignments (ranges from -1 to 1 where 1 being a perfect clustering and 0 being a random clustering). In a nutshell, all the three performance metrics are commonly used but all have its pros and cons; for instance, the
4https://adni.loni.usc.edu 5This extension is a representative of recently proposed deep learning approaches for clustering of both static data (Xie et al., 2017; Yang et al., 2017) and time-series data (Baytas et al., 2017; Madiraju et al., 2018) since these methods are built upon the same concept – that is, applying deep networks for dimensionality reduction to conduct conventional clustering methods, e.g., K-means.
purity score easily converges to 1 when there are as many clusters as data samples. Thus, using them together suffices to demonstrate the effectiveness of the clustering methods.
To assess the prediction performance of the identified predictive clusters, we utilized both area under receiver operator characteristic curve (AUROC) and area under precision-recall curve (AUPRC) based on the label predictions of each cluster and the ground-truth binary labels on the future outcomes of interest. Note that the prediction performance is available only for the benchmarks that incorporate the label information during training.

heading: 5.4 CLUSTERING PERFORMANCE
text: We start with a simple scenario where the true class (i.e., the ground-truth cluster label) is available and the number of classes is tractable. In particular, we setC = 23 = 8 based on the binary labels for the development of three common comorbidities of cystic fibrosis – diabetes, ABPA, and intestinal obstruction – in the next year for the UKCF dataet and C = 3 based on the mutually exclusive three diagnostic groups for the ADNI dataset. We compare AC-TPC against the aforementioned benchmarks with respect to the clustering and prediction performance in Table 1.
As shown in Table 1, AC-TPC achieved performance gain over all the tested benchmarks in terms of both clustering and prediction performance – where most of the improvements were statistically significant with p-value < 0.01 or p-value < 0.05 – for both datasets. Importantly, clustering methods – i.e., KM-DTW, DCN-S2S, and SOM-VAE – that do not associate with the future outcomes of interest identified clusters that provide little prognostic value on the future outcomes (note that the true class is derived from the future outcome of interest). This is clearly shown by the ARI value near 0 which indicates that the identified clusters have no difference with random assignments. Therefore, similar sequences with respect to the latent representations tailored for reconstruction or with respect to the shape-based measurement using DTW can have very different future outcomes.
In Figure 2, we further investigate the purity score, NMI, and ARI by varying the number of clusters K from 4 to 16 on the UKCF dataset in the same setting with that stated above (i.e., C = 8).
Here, the three methods – i.e., KM-DTW, DCN-S2S, and SOM-VAE – are excluded for better visualization. As we can see in Figure 2, our model rarely incur performance loss in both NMI and ARI while the benchmarks (except for SOM-VAE-P) showed significant decrease in the performance as K increased (higher than C). This is because the number of clusters identified by AC-TPC (i.e., the number of activated clusters where we define cluster k is activated if |C(k)| > 0) was the same with C most of the times, while the DCN-based methods identified exactly K clusters (due to the K-means). Since the NMI and ARI are adjusted for the number of clusters, a smaller number of identified clusters yields, if everything being equal, a higher performance. In contrast, while our model achieved the same purity score for K ≥ 8, the benchmark showed improved performance as K increased since the purity score does not penalize having many clusters. This is an important property of AC-TPC that we do not need to know a priori what the number of cluster is which is a common practical challenge of applying the conventional clustering methods (e.g., K-means).
The performance gain of our model over SOM-VAE-P (and, our analysis is the same for SOMVAE) comes from two possible sources: i) SOM-VAE-P mainly focuses on visualizing the input with SOM which makes both the encoder and embeddings less flexible – this is why it performed better with higherK – and ii) the Markov property can be a too strict assumption for time-series data especially in clinical settings where a patient’s medical history is informative for predicting his/her future clinical outcomes (Ranganath et al., 2016).

heading: 5.5 CONTRIBUTIONS OF THE AUXILIARY LOSS FUNCTIONS
text: As described in Section 3.1, we introduced two auxiliary loss functions – the sample-wise entropy of cluster assignment (4) and the embedding separation loss (5) – to avoid trivial solution that may arise in identifying the predictive clusters. To analyze the contribution of each auxiliary loss function, we report the average number of activated clusters, clustering performance, and prediction performance on the UKCF dataset with 3 comorbidities as described in Section 5.4. Throughout the experiment, we set K = 16 – which is larger than C – to find the contribution of these loss functions to the number of activated clusters.
As we can see in Table 2, both auxiliary loss functions make important contributions in improving the quality of predictive clustering. More specifically, the sample-wise entropy encourages the selector to choose one dominant cluster. Thus, as we can see results with α = 0, without the sample-wise entropy, our selector assigns an equal probability to all 16 clusters which results in a trivial solution. We observed that, by augmenting the embedding separation loss (5), AC-TPC identifies a smaller number of clusters owing to the well-separated embeddings; in Appendix E, we further investigate the usefulness of (5) in identifying the number of clusters in a data-driven fashion.

heading: 5.6 TARGETING MULTIPLE FUTURE OUTCOMES – A PRACTICAL SCENARIO
text: In this experiment, we focus on a more practical scenario where the future outcome of interest is high-dimensional and the number of classes based on all the possible combinations of future outcomes becomes intractable. For example, suppose that we are interested in the development of M comorbidities in the next year whose possible combinations grow exponentially C = 2M . Interpreting such a large number of patient subgroups will be a daunting task which hinders the understanding of underlying disease progression. Since different comorbidities may share common driving factors (Ronan et al., 2017), we hope our model to identify much smaller underlying (latent) clusters that govern the development of comorbidities. Here, to incorporate with M comorbidities (i.e., M binary labels), we redefine the output space as Y = {0, 1}M and modify the predictor and loss functions, accordingly, as described in Appendix A.
Throughout this experiment, we aim at identifying subgroups of patients that are associated with the next-year development of 22 different comorbidities in the UKCF dataset. In Table 3, we reported 12 identified clusters – on average, the number of activated clusters were 13.6 – and the top three frequent comorbidities developed in the next year since the latest observation and the corresponding frequency; please refer to Appendix E for a full list. Here, the frequency is calculated in a clusterspecific fashion based on the true label. As we can see in Table 3, the identified clusters displayed very different label distributions; that is, the combination of comorbidities as well as their frequency were very different across the clusters. For example, patients in Cluster 1 experienced low-risk of developing any comorbities in the next year while 85% of patients in Cluster 0 experienced diabetes in the next year.
In Figure 3, we further investigated subgroups of patients – Cluster 0, 5, 7, 8, and 10 – who had high risk of developing diabetes in the next year. Although all these clusters displayed high risk of diabetes, the frequency of other co-occurred comorbidities was significantly different across the clusters. In particular, around 89% of the patients in Cluster 5 experienced asthma in the next year while it was less than 3% of the patients in the other cluster. Interestingly, “leukotriene” – a medicine commonly used to manage asthma – and “FEV1% predicted” – a measure of lung function – were the two most different input features between patients in Cluster 5 and those in the other clusters. We observed similar findings in Cluster 7 with ABPA, Cluster 8 with liver disease, and Cluster 10 with osteopenia. Therefore, by grouping patients who are likely to develop a similar set of comorbidities, our method identified clusters that can be translated into actionable information for clinical decision-making.

heading: 6 CONCLUSION
text: In this paper, we introduced AC-TPC, a novel deep learning approach for predictive clustering of time-series data. We carefully defined novel loss functions to encourage each cluster to have homogeneous future outcomes (e.g., adverse events, the onset of comorbidities, etc.) and designed optimization procedures to address challenges and to avoid trivial solution in identifying such cluster assignments and the centroids. Throughout the experiments on two real-world datasets, we showed that our model achieves superior clustering performance over state-of-the-art and identifies meaningful clusters that can be translated into actionable information for clinical decision-making.

heading: B DETAILED DERIVATION OF (7)
text: To derive the gradient of the predictive clustering loss in (7) with respect ωA = [θ, ψ], we utilized the ideas from actor-critic models (Konda and Tsitsiklis, 2000). The detailed derivation of the former term in (7) is described below (for notational simplicity, we omit the expectation on Ex,y∼pXY ):
∇ω ( T∑ t=1 Est∼Cat(πt) [ `1(yt, ȳt) ]) = ∇ω ( T∑ t=1 ∑ st∈K πt(st)`1(yt, ȳt) )
= T∑ t=1 ∑ st∈K ∇ωπt(st)`1(yt, ȳt)
= T∑ t=1 ∑ st∈K ∇ωπt(st) πt(st) πt(st)`1(yt, ȳt)
= T∑ t=1 ∑ st∈K πt(st)`1(yt, ȳt)∇ω log πt(st)
= T∑ t=1 Est∼Cat(πt) [ `1(yt, ȳt)∇ω log πt(st) ]
(11)

heading: C DETAILS ON THE DATASETS
text: C.1 UKCF DATASET
UK Cystic Fibrosis registry (UKCF)6 records annual follow-ups for 5,171 adult patients (aged 18 years or older) over the period from 2008 and 2015, with a total of 25,012 hospital visits. Each patient is associated with 89 variables (i.e., 11 static and 78 time-varying features), including information on demographics and genetic mutations, bacterial infections, lung function scores, therapeutic managements, and diagnosis on comorbidities. The detailed statistics are given in Table 4.
C.2 ADNI DATASET
Alzheimer’s Disease Neuroimaging Initiative (ADNI)7 study consists of 1,346 patients with a total of 11,651 hospital visits, which tracks the disease progression via follow-up observations at 6 months interval. Each patient is associated with 21 variables (i.e., 5 static and 16 time-varying features), including information on demographics, biomarkers on brain functions, and cognitive test results.
6https://www.cysticfibrosis.org.uk/the-work-we-do/uk-cf-registry 7https://adni.loni.usc.edu
The three diagnostic groups were normal brain functioning (0.55), mild cognitive impairment (0.43), and Alzheimer’s disease (0.02). The detailed statistics are given in Table 5.

heading: D DETAILS ON THE BENCHMARKS
text: We compared AC-TPC in the experiments with clustering methods ranging from conventional approaches based on K-means to the state-of-the-art approaches based on deep neural networks. The details of how we implemented the benchmarks are described as the following: • Dynamic time warping followed by K-means8: Dynamic time warping (DTW) is utilized to
quantify pairwise distance between two variable-length sequences and, then,K-means is applied (denoted as KM-DTW).
• K-means with deep neural networks: To handle variable-length time-series data, we utilized an encoder-predictor network as depicted in Figure 4b and trained the network based on (10) for dimensionality reduction; this is to provide fixed-length and low-dimensional representations for time-series. Then, we appliedK-means on the latent encodings z (denoted as KM-E2P (Z)) and on the predicted label distributions ŷ (denoted as KM-E2P (Y)), respectively. We implemented the encoder and predictor of KM-E2P with the same network architectures with those of our model: the encoder is a single-layer LSTM with 50 nodes and the decoder is a two-layered fully-connected network with 50 nodes in each layer.
• Extensions of DCN9 (Yang et al., 2017): Since the DCN is designed for static data, we utilized a sequence-to-sequence model in Figure 4a for the encoder-decoder network as an extension to
8https://github.com/rtavenar/tslearn 9https://github.com/boyangumn/DCN
incorporate time-series data (denoted as DCN-S2S) and trained the network based on the reconstruction loss (using the reconstructed input sequence x̂1:t). For implementing DCN-S2S, we used a single-layer LSTM with 50 nodes for both the encoder and the decoder. And, we augmented a fully-connected layer with 50 nodes is used to reconstruct the original input sequence from the latent representation of the decoder. In addition, since predictive clustering is associated with the label distribution, we compared a DCN whose encoder-decoder structure is replaced with our encoder-predictor network in Figure 4b (denoted as DCN-E2P) to focus the dimensionality reduction – and, thus, finding latent encodings where clustering is performed – on the information for predicting the label distribution. We implemented the encoder and predictor of DCN-E2P with the same network architectures with those of our model as described in Section 5.
• SOM-VAE10 (Fortuin et al., 2019): We compare with SOM-VAE – though, this method is oriented towards visualization of input data via SOM – since it naturally clusters time-series data assuming Markov property (denoted as SOM-VAE). We replace the original CNN architecture of the encoder and the decoder with three-layered fully-connected network with 50 nodes in each layer, respectively. The network architecture is depicted in Figure 4c where x̂t and x̄t indicate the reconstructed inputs based on the encoding zt and the embedding et at time t, respectively. In addition, we compare with a variation of SOM-VAE by replacing the decoder with the predictor to encourage the latent encoding to capture information for predicting the label distribution (denoted as SOM-VAE-P). For the implementation, we replaced the decoder of SOM-VAE with our predictor which is a two-layered fully-connected layer with 50 nodes in each layer to predict the label distribution as illustrated in Figure 4d. Here, ŷt and ȳt indicate the predicted labels based on the encoding zt and the embedding et at time t, respectively. For both cases, we used the default values for balancing coefficients of SOM-VAE and the dimension of SOM to be equal to K.
We compared and summarized major components of the benchmarks in Table 6.

heading: E ADDITIONAL EXPERIMENTS
text: E.1 ADDITIONAL RESULTS ON TARGETING MULTIPLE FUTURE OUTCOMES
Throughout the experiment in Section 5.6, we identified 12 subgroups of patients that are associated with the next-year development of 22 different comorbidities in the UKCF dataset. In Table 7, we reported 12 identified clusters and the full list of comorbidities developed in the next year since the latest observation and the corresponding frequency. Here, the frequency is calculated in a clusterspecific fashion based on the true label.
E.2 TRADE-OFF BETWEEN CLUSTERING AND PREDICTION PERFORMANCE
In predictive clustering, the trade-off between the clustering performance (for better interpretability) – which quantifies how the data samples are homogeneous within each cluster and heterogeneous across clusters with respect to the future outcomes of interest – and the prediction performance is a common issue. The most important parameter that governs this trade-off is the number of clusters. More specifically, increasing the number of clusters will make the predictive clusters have higher diversity to represent the output distribution and, thus, will increase the prediction performance while decreasing the clustering performance. One extreme example is that there are as many clusters as data samples which will make the identified clusters fully individualized; as a consequence, each cluster will lose the interpretability as it no longer groups similar data samples.
To highlight this trade-off, we conducted experiments under the same experimental setup with that of Section 5.6 where our aim is to identify underlying (unknown) clusters when the future outcome of interest is high-dimensional. For the performance measures, we utilized the AUROC and AUPRC to assess the prediction performance, and utilized the average Silhouette index (SI) (Rousseeuw, 1987) – a widely used measure of how similar a member is to its own cluster (homogeneity within a cluster) compared to other clusters (heterogeneity across clusters ) when the ground-truth cluster labels are
10https://github.com/ratschlab/SOM-VAE
not available – to assess the identified clusters. Formally, the SI for a subsequence xn1:t ∈ Ck can be given as follows:
SI(n) = b(n)− a(n) max ( a(n), b(n) ) (12) where a(n) = 1|Ck|−1 ∑ m 6=n ‖ynt − ymt ‖1 and b(n) = mink′ 6=k 1 |Ck′ | ∑ m∈Ck′ ‖ynt − ymt ‖1. Here, we used the L1-distance between the ground-truth labels of the future outcomes of interest since our goal is to group input subsequences with similar future outcomes. To control the number of identified clusters (i.e., the activated clusters) of our method, we set β = 0 (since the embedding separation loss in (5) controls the activation of clusters) and reported the performance by increasing the number of possible clusters K.
As can be seen in Figure 5, the prediction performance increased with an increased number of identified clusters due to the higher diversity to represent the label distribution while making the identified clusters less interpretable (i.e., the cohesion and separation among clusters become ambiguous as seen in the low average SI). On the other hand, when we set β = 1.0 (which is selected based on the validation loss in 3), our method consistently identified a similar number of clusters for K > 20, i.e., 13.8 on average, in a data-driven fashion and provided slightly reduced prediction performance with significantly better interpretability, i.e., the average SI 0.120 on average. This highlights the usefulness of (5) which helps to identify clusters with different label distributions.

heading: F PSEUDO-CODE OF AC-TPC
text: As illustrated in Section 3.2, AC-TPC is trained in an iterative fashion. We provide the pseudo-code for optimizing our model in Algorithm 1 and that for initializing the parameters in Algorithm 2.
Algorithm 1 Pseudo-code for Optimizing AC-TPC
Input: Dataset D = {(xnt , ynt )T n t=1}Nn=1, number of clusters K, coefficients (α, β), learning rate (ηA, ηC , ηE), mini-batch size nmb, and update step M Output: AC-TPC parameters (θ, ψ, φ) and the embedding dictionary E Initialize parameters (θ, ψ, φ) and the embedding dictionary E via Algorithm 2
repeat Optimize the Encoder, Selector, and Predictor for m = 1, · · · ,M do
Sample a mini-batch of nmb data samples: {(xnt , ynt )T n t=1}nmbn=1 ∼ D for n = 1, · · · , nmb do
Calculate the assignment probability: πnt = [πnt (1) · · ·πnt (K)]← hψ(fθ(xn1:t)) Draw the cluster assignment: snt ∼ Cat(πnt ) Calculate the label distributions: ȳnt ← gφ(e(snt )) and ŷnt ← gφ(fθ(xn1:t))
end for Update the encoder fθ and selector hψ:
θ ← θ − ηA
( 1
nmb nmb∑ n=1 Tn∑ t=1 `1(y n t , ȳ n t )∇θ log πnt (snt )− α∇θ K∑ k=1 πnt (k) log π n t (k)
)
ψ ← ψ − ηA
( 1
nmb nmb∑ n=1 Tn∑ t=1 `1(y n t , ȳ n t )∇ψ log πnt (snt )− α∇ψ K∑ k=1 πnt (k) log π n t (k) ) Update the predictor gφ:
φ← φ− ηC 1
nmb nmb∑ n=1 Tn∑ t=1 ∇φ`1(ynt , ȳnt )
end for
Optimize the Cluster Centroids for m = 1, · · · ,M do
Sample a mini-batch of nmb data samples: {(xnt , ynt )T n t=1}nmbn=1 ∼ D for n = 1, · · · , nmb do
Calculate the assignment probability: πnt = [πnt (1) · · ·πnt (K)]← hψ(fθ(xn1:t)) Draw the cluster assignment: snt ∼ Cat(πnt ) Calculate the label distributions: ȳnt ← gφ(e(snt ))
end for for k = 1, · · · ,K do
Update the embeddings e(k):
e(k)← e(k)− ηE
( 1
nmb nmb∑ n=1 Tn∑ t=1 ∇e(k)`1(ynt , ȳnt )− γ K∑ k′=1 k′ 6=k ∇e(k)`1 ( gφ(e(k)), gφ(e(k
′)) ))
end for Update the embedding dictionary: E ← {e(1), . . . e(K)}
end for until convergence
Algorithm 2 Pseudo-code for pre-training AC-TPC
Input: Dataset D = {(xnt , ynt )T n t=1}Nn=1, number of clusters K, learning rate η, mini-batch size nmb Output: AC-TPC parameters (θ, ψ, φ) and the embedding dictionary E Initialize parameters (θ, ψ, φ) via Xavier Initializer
Pre-train the Encoder and Predictor repeat
Sample a mini-batch of nmb data samples: {(xnt , ynt )T n t=1}nmbn=1 ∼ D for n = 1, · · · , nmb do
Calculate the label distributions: ŷnt ← gφ(fθ(xn1:t)) end for
θ ← θ − η 1 nmb nmb∑ n=1 Tn∑ t=1 ∇θ`1(ynt , ŷnt ) φ← φ− η 1 nmb nmb∑ n=1 Tn∑ t=1 ∇φ`1(ynt , ŷnt )
until convergence
Initialize the Cluster Centroids Calculate the embedding dictionary E and initial cluster assignments cnt
E , {{cnt }T n t=1}Nn=1 ← K-means({{znt }T n t=1}Nn=1,K)
Pre-train the Selector repeat
Sample a mini-batch of nmb data samples: {(xnt , ynt )T n t=1}nmbn=1 ∼ D for n = 1, · · · , nmb do
Calculate the assignment probability: πnt = [πnt (1) · · ·πnt (K)]← hψ(fθ(xn1:t)) end for Update the selector hψ:
ψ ← ψ + η 1 nmb nmb∑ n=1 Tn∑ t=1 K∑ k=1 cnt (k) log π n t (k)
until convergence

Review summaries:
This work proposes a method to perform clustering on a time-series data for prediction purposes, unlike the classical clustering where it is done in an unsupervised manner. The authors use an encoder (RNN) to process the time-series medical records, a selector to sample the cluster label for each encoding, and a predictor to predict the labels based on the selected cluster centroids. Since the sampling process prevents the authors from using back-prop, they employ an actor-critic method.