# Clustering Interval-Censored Time-Series
for Disease Phenotyping

Irene Y. Chen Massachusetts Institute of Technology. Email:
<EMAIL_ADDRESS>Rahul G. Krishnan University of Toronto. Work done
partially while at Microsoft Research New England. Email:
<EMAIL_ADDRESS>David Sontag Massachusetts Institute of Technology.
Email<EMAIL_ADDRESS>
###### Abstract

Unsupervised learning is often used to uncover clusters in data. However,
different kinds of noise may impede the discovery of useful patterns from
real-world time-series data. In this work, we focus on mitigating the
interference of interval censoring in the task of clustering for disease
phenotyping. We develop a deep generative, continuous-time model of time-
series data that clusters time-series while correcting for censorship time. We
provide conditions under which clusters and the amount of delayed entry may be
identified from data under a noiseless model. On synthetic data, we
demonstrate accurate, stable, and interpretable results that outperform
several benchmarks. On real-world clinical datasets of heart failure and
Parkinson’s disease patients, we study how interval censoring can adversely
affect the task of disease phenotyping. Our model corrects for this source of
error and recovers known clinical subtypes.

##  1 Introduction

Cluster analysis of time-series data is a task of interest across a variety of
scientific disciplines including biology [Listgarten et al., 2007],
meteorology [Camargo et al., 2007], and astrophysics [Rebbapragada et al.,
2009]. Automating the discovery of latent patterns in real-world data can be
challenging due to noise. We focus on mitigating errata in pattern discovery
from interval censoring [Miller Jr, 2011].

Interval censoring arises when time-series data are only observed within a
known interval. The difference between an observed time value and the true
time value, known as left-censorship, and the lack of a common outcome against
which to align, known as right-censorship, can lead common techniques for
clustering to erroneous conclusions about the underlying patterns. To address
this, practitioners must often manually align data to a meaningful start time
in order to find non-trivial groups via unsupervised clustering—a process
whose difficulty can range from expensive and time-consuming in some problems
to infeasible in others. In this work, we develop a machine learning algorithm
that clusters time-series data while simultaneously correcting for interval
censoring. In doing so, we automate the time-consuming process of manual data
alignment and use our method to reveal structure that would otherwise not be
found by straightforward application of clustering analysis.

This paper is motivated by disease phenotyping in healthcare. Many diseases
are biologically heterogeneous despite a common diagnosis—for example autism
[Doshi-Velez et al., 2014], heart failure [Shah et al., 2015], and Parkinson’s
disease [Fereshtehnejad et al., 2017]. The variation in biomarkers (e.g.,
glucose or creatinine) across patients can stem from different patient
subtypes that manifest in distinct disease trajectories. Scientists seek to
understand this disease heterogeneity by identifying groups of people whose
biomarkers behave similarly. For example, cardiologists use a measurement
called ejection fraction as a heuristic to separate heart failure patients
into two categories [Owan et al., 2006], with at least one of the two subtypes
believed to be heterogeneous [Shah et al., 2015]. To better understand patient
heterogeneity, cardiologists may turn to longitudinal, observational, and
often irregularly-measured patient data for disease phenotype discovery.

(a) width=0.3 (b) width=0.4 (c) width=0.1

Figure 1:  (a) Patient data can be interval-censored, meaning longitudinal
data can be missing based on entry to the dataset, e.g., first diagnosis, and
the data may lack a common outcome against which to align. Patients may enter
the dataset at any stage of the disease. (b) Interval censoring can make
clustering patient time-series data challenging because data may be aligned
incorrectly, e.g., first diagnosis. We seek to understand disease
heterogeneity by inferring subtypes after correcting for misalignment. (c)
Graphical model of SubLign.

Interval censoring in healthcare data presents a significant challenge for
disease phenotyping. Left-censorship in clinical datasets can occur when
patients have delayed entry, meaning data are unavailable before a diagnosis
or first hospital visit. Factors including geographic proximity to a hospital
[Chan et al., 2006], financial access to care [Miller and Wherry, 2017] or
mistrust of the healthcare system [Brandon et al., 2005] can affect when a
patient seeks medical help and consequently the beginning of data availability
with respect to the underlying progression of their disease. Other datasets
align patients by death time, but right-censorship restricts sample size to
patients who have died. Since many factors affect mortality, other causes of
death may confound results. For heart failure, a chronic disease that
progresses over many years with a large range of onset ages and survival
outcomes, this interval censoring can confound attempts to analyse disease
heterogeneity using observational data.

For a simplified illustration of the problem, see Figure 1(a) which depicts
the common reality of observational health data whereas Figure 1(b) depicts
the idealized latent substructure we would like to identify. Existing
subtyping models applied to clinical data assume (potentially erroneously)
that patients are aligned at entry into the dataset or study. The problem with
such an assumption is that the disparity between true disease stage and
observed observation time can result in unsupervised learning algorithms
uncovering the wrong, or perhaps less interesting, structure. For example, a
naive clustering algorithm might simply return clusters corresponding to the
disease stage at entry into the study, which may simply recapitulate some of
the biases mentioned earlier.

To that end, our paper makes the following contributions:

  1. 1.

We introduce and formalize the problem of cluster recovery from interval-
censored time-series data.

  2. 2.

We introduce a practical algorithm, SubLign, based on variational learning of
a deep generative model, which:

     * •

Makes no assumption on the distribution of delayed entry alignment values
besides extrema.

     * •

Operates on multivariate time-series with varying lengths and missing values,
both characteristics often found in real-world datasets.

     * •

Is unsupervised, meaning that neither subtype labels nor alignment values are
provided during training.

  3. 3.

We prove an identifiability result showing that, in a noiseless setting, both
the degree of delayed entry from left-censorship and the subtype identity are
recoverable.

We show robust quantitative results on synthetic data where, over multiple
runs, our method outperforms baselines for subtyping and patient alignment. On
two real-world clinical datasets—heart failure and Parkinson’s disease—we
automatically recover known clinical subtypes without manual alignment or
expert knowledge.

##  2 Related Work

Learning alignment and clustering has been studied in fields across computer
vision, signal processing, and health. Approaches often make assumptions
including few discrete time steps [Young et al., 2018, Huopaniemi et al.,
2014]; a single piecewise linear function [Young et al., 2018] or Gaussian
mixture model [Huopaniemi et al., 2014]; significantly more samples per object
than number of objects [Mattar et al., 2012, Gaffney and Smyth, 2005]; very
small windows of potential misalignment [Liu et al., 2009, Listgarten et al.,
2007]; or known lag time [Li et al., 2011]. Methods that directly measure
similarity between time-series, e.g., dynamic time warping [Cuturi, 2011] or
methods that aggregate multiple imputation methods [Faucheux et al., 2021] can
also be used for clustering time-series data. Our method aims to cluster
interval-censored multivariate time series without these constraints.

Clinicians and scientists learn disease subtypes to better understand
heterogeneity in disease progression in a process known as disease
phenotyping. Existing approaches often rely on the assumption that the
observed measurements are aligned — and therefore not censored. Researchers
then apply clustering techniques like hierarchical clustering of time series
[Doshi-Velez et al., 2014], affinity clustering [Luo et al., 2020], or matrix
factorization [Udler et al., 2018, Perros et al., 2017]. Other models define
disease subtypes as stages of disease progression [Alaa and van der Schaar,
2019]. For this work, we define disease subtypes as distinct from disease
stage and jointly learn both.

##  3 SubLign: Subtype & Align

There are two stages to SubLign. First, we learn a generative model of the
observed data which disentangles variation in the observed data due to delayed
entry from variation related to subtype identity. Second, we infer subtype
representations and (optionally) cluster the representations to obtain the
explicit subtype identity for each time-series. Figure 1(c) describes the
graphical model, and Algorithm 1 depicts the pseudocode for this procedure.

Algorithm 1 SubLign

1: Input: Observation times $X\in\mathbb{R}^{N\times M}$, biomarkers
$Y\in\mathbb{R}^{N\times M\times D}$

2: Output: $\tau_{k}$ for each subtype and $\hat{\delta_{i}}$ for each patient

3: Step $1$: Learning

4: repeat

5: Encode time-series: $h_{i}=\text{RNN}([X_{i},Y_{i}])\;\forall
i\in\\{1,\ldots,N\\}$

6: Compute variational distribution
$q(Z_{i}|X_{i},Y_{i})=\mathcal{N}(\mu(h_{i};\phi_{2}),\Sigma(h_{i};\phi_{3}))\forall
i\in\\{1,\ldots,N\\}$

7: for patient $i=1$ to $N$ do

8: Run grid-search to find
$\hat{\delta_{i}}=\arg\max_{q(\delta_{i})}\mathcal{L}(Y_{i}|X_{i};\gamma,\phi,q(\delta_{i}))$
(Eq. 3)

9: end for

10: Update $\gamma,\phi$ via stochastic gradient ascent on
$\mathcal{L}(Y|X;\gamma,\phi,\hat{\delta})$

11: until convergence

12: Step $2$: Inference and Clustering

13: Infer $\mathcal{Z}=\\{z_{i}|z_{i}=\mu(h_{i};\phi_{2})\\}$ for
$X_{i},Y_{i}$

14: Find $K$ clusters using $k$-means on $\mathcal{Z}$ and compute cluster
centers $\mu_{k}$

15: Infer parameters of subtype trajectories $\tau_{k}=g(\mu_{k})$

#### SubLign Generative Model

Consider the following setup. We observe $N$ multivariate time-series (one for
each patient), each of length up to $M$:
$[(x_{1,1},y_{1,1}),\ldots,(x_{1,M},y_{1,M})],\ldots[(x_{N,1},y_{N,1}),\ldots,\\\
(x_{N,M},y_{N,M})]$. $y_{i,m}\in\mathbb{R}^{D}$ is a vector of observations
for time-series $i$ at time-stamp $x_{i,m}\in\mathbb{R}^{+}$. We denote
collections of observations as $Y_{i}=\\{y_{i,1},\ldots,y_{i,M}\\}$ and time-
stamps as $X_{i}=\\{x_{i,1},\ldots,x_{i,M}\\}$ for patient $i$.

Figure 1(c) depicts the graphical model corresponding to the latent-variable
generative model of continuous-time multivariate data:

$\displaystyle\forall i=\\{1,\ldots,N\\},\;\forall m\in\\{1,\ldots,M\\},$
$\displaystyle\forall d\in
D_{i,m},\;\delta_{i}\sim\text{Cat}(\mathcal{D}),z_{i}\sim\mathcal{N}(\mathbf{0},\mathbb{I}),\;\Theta=g(z_{i};\gamma),$
$\displaystyle\overline{y}_{i,m}[d]=f(\kappa(x_{i,m}+\delta_{i};\Theta[d])),\;y_{i,m}\sim\mathcal{N}(\overline{y}_{i,m}[d],\mathbf{I})$
(1)

We drop indices denoting patient and dimension where unnecessary. $D_{i,m}$
denotes the set (and $|D_{i,m}|$ denotes the number) of observed biomarkers
for patient $i$ at their $m$-th observation. To accommodate missing data, not
all biomarkers are required to be measured at every observation. Each delayed
entry value $\delta_{i}\in\mathbb{R^{+}}$ has a maximum alignment deviation
value $\delta^{+}$ over all time-series. We discretize the closed interval
$[0,\delta^{+}]$ as $\mathcal{D}=[0,\epsilon,2\epsilon,\ldots,\delta^{+}]$
with hyperparameter $\epsilon$ and use a categorical distribution over
$\mathcal{D}$ with uniform probabilities over each element as our prior over
$\delta_{i}$. Function $g:\mathbb{R}^{N_{z}}\to\mathbb{R}^{D\times(P+1)}$ has
parameters $\gamma$ and maps from the latent variable to
$\Theta\in\mathbb{R}^{D\times(P+1)}$, a matrix of parameters for $D$
polynomials, each of degree $P$. $f$ is a known link function that describes
how observed values relate to observed time-points.

Parameterization We discuss the specific parameterizations of Equation 3. In
the context of our motivating application of disease phenotyping, these
functions represent common characteristics in the progression of patient
biomarkers.

Link function and polynomials: For $f$ and $P$, we study the following
choices: Sigmoid: $P=1$ and $f(x)=\frac{1}{1+\exp(-x)}$ and Quadratic: $P=2$
and $f(x)=x$. The sigmoid function can represent bounded and monotonically
increasing clinical variables. The quadratic function represents cases where
disease severity, as measured by biomarkers, decreases (likely in response to
therapy) and then increases (once therapy fails), or vice versa. Other choices
for $P,f$ are permissible as long as they are differentiable with respect to
the model parameters. We allow for the possibility that individual biomarkers
have different parameterizations.

Modeling polynomial parameters: We parameterize $g(z_{i};\gamma)$ using a two
layer neural network with ReLU activation functions with parameters $\gamma$.
To be concrete, if $D=1,P=2$, and $f$ is the sigmoid function, then the
outputs of $g$ are $[\beta_{0}(z),\beta_{1}(z)]$ and
$y=\frac{1}{1+\exp^{-(\beta_{0}(z)x+\beta_{1}(z))}}$. Similarly, if $D=1,P=1$,
and $f$ is the quadratic function then the outputs of $g$ are
$[a(z),b(z),c(z)]$ and $y=a(z)x^{2}+b(z)x+c(z)$.

###  3.1 Step $1$: Learning

We learn the parameters $\gamma$ of the model in Equation 3 via maximum
likelihood estimation. Since the model is a non-linear latent variable model,
we maximize a variational lower bound on the conditional likelihood of data
given the time-stamps corresponding to observations.

$\displaystyle\log\prod_{m=1}^{M}\prod_{d\in
D_{i,m}}p(y_{i,m}[d]|x_{i,m};\gamma)=\log p(Y_{i}|X_{i};\gamma)$ (2)
$\displaystyle\geq\mathcal{L}(Y_{i}|X_{i};\gamma,\phi)$
$\displaystyle=\mathbb{E}_{q(Z_{i}|X_{i},Y_{i};\phi)}\Big{[}\log\sum_{\delta_{i}}p(Y_{i}|X_{i},\delta_{i},Z_{i};\gamma)p(\delta_{i})+$
$\displaystyle\quad\quad\quad+\log\frac{p(Z_{i})}{q(Z_{i}|X_{i},Y_{i};\phi)}\Big{]}\geq\mathcal{L}(Y_{i}|X_{i};\gamma,\phi,q(\delta_{i}))$
$\displaystyle=\mathbb{E}_{q(Z_{i}|X_{i},Y_{i};\phi)}\Big{[}\mathbb{E}_{q(\delta_{i})}\Big{[}\log
p(Y_{i}|X_{i},\delta_{i},Z_{i};\gamma)$
$\displaystyle\quad\quad\quad+\log\frac{p(\delta_{i})}{q(\delta_{i})}\Big{]}+\log\frac{p(Z_{i})}{q(Z_{i}|X_{i},Y_{i};\phi)}\Big{]}$
(3)

The first lower bound uses a variational distribution for $Z$ parameterized
via an inference network [Kingma and Welling, 2013, Rezende et al., 2014] with
parameters $\phi$. The second lower bound is a variational distribution over
$\delta$. The function $q(\delta_{i})$ parameterizes the space of one-hot
distributions in $\mathcal{D}$, i.e., a categorical distribution over discrete
choices of $\delta_{i}$.

Our learning algorithm alternates between two steps. We first maximize the
lower bound using subgradient ascent. To do so, we solve:
$\hat{\delta_{i}}=\arg\max_{q(\delta_{i})}\mathcal{L}(Y_{i}|X_{i};\gamma,\phi,q(\delta_{i}))$.
For our choice of variational distribution, this maximization can be performed
via a grid search. We then derive gradients
$\nabla_{\gamma,\phi}\mathcal{L}(Y_{i}|X_{i};\gamma,\phi,\hat{\delta}_{i})$ to
update the generative model and inference network via stochastic gradient
ascent.

###  3.2 Step $2$: Obtaining Learned Subtypes

After learning the model, we may re-use the inference network to predict the
latent variable $z_{i}=\mu(h_{i};\phi_{2})$ for each patient in the training
set. Combining $z_{i}$ across all time-series gives us the set $\mathcal{Z}$.
When reasonable, we refer to $\mu(h_{i};\phi_{2})$ as $\mu_{i}$.

Although latent variable $z_{i}$ encodes latent structure from each time-
series, we may be interested in explicit subtypes for a given value of $K$. To
obtain discrete subtypes, we can run clustering algorithms on $\mathcal{Z}$ to
obtain $K$ cluster centers $\\{\mu_{1},\ldots,\mu_{K}\\}$. Because we use a
Gaussian prior for our biomarker values, measuring distances in the space of
the latent variable can be done with the Euclidean norm, making the $k$-means
a reasonable choice of clustering algorithm.

We compute $\\{\tau_{1},\ldots,\tau_{K}\\}$ where $\tau_{k}=g(\mu_{k})$ as the
progression-patterns corresponding to each of the discrete subtypes of the
disease. For example, if $f\circ\kappa$ is linear, then we obtain $K$
different slopes and biases, each of which describes how the time-series
behaves in that subtype. In practice, $K$ may be chosen based on domain
knowledge; alternatively, qualitative results can be assessed for each version
of $K$, e.g., by plotting the corresponding $f$ functions.

###  3.3 Remarks on SubLign

Role of the latent variable: The latent variable $z$ plays an important role
in quantifying how each biomarker behaves. Each time-series’s latent variable
is used to predict the parameters of $D$ polynomial functions and
$f\circ\kappa$ maps from observation times onto the observed biomarkers. Time-
series whose representation space $z$ are _close_ hail from the same subtype,
and consequently manifest similar patterns in their biomarkers. This variation
in $z$ results in variation in the parameters $\Theta$ and therefore in
variation among the data as a function of the time-points. As an illustration,
in Figure 1(b) in the blue phenotype, $f$ is the sigmoid function. Depending
on the latent space, we could imagine a one-dimensional $z$ where $z<0$
represents the curve (and subtype) in blue and $z\geq 0$ to represent the
curve in red. The value of $\delta_{i}$ indicates the degree of delayed entry
associated with each time-series. The delayed entry from interval censoring is
corrected by applying the scalar $\delta_{i}$ element-wise to $X_{i}$ and then
transforming it by $f$.

Scalability: The runtime of SubLign is impacted by the grid search over model
parameter $\delta_{i}$. The corresponding lines 7-9 in Algorithm 1 have
complexity $\mathcal{O}(NSF)$ where $F$ is the complexity associated with a
single forward pass of the inference network and the generative model, $N$ is
the number of examples, and $S=\delta^{+}/\epsilon$ is the number of time
steps. The model practitioners may therefore balance computational resources
with $S$. In our experiments, we found comparable performance for $S$ as low
as 5.

Real-world clinical data: SubLign is motivated by, and designed to capture,
variation in clinical biomarkers while taking into account the challenges of
clinical data. Observational healthcare data are often irregularly spaced, and
contains missingness. The use of a continuous time model allows us to
naturally handle the former issue since we only maximize the likelihood of
data corresponding to time-points where they are observed. When a single
biomarker is missing while others are observed, it may be marginalized out (by
ignoring the corresponding loss term). r Accommodating different kinds of
censorship: Equation 3 naturally characterizes delayed entry arising from
left-censorship. SubLign can also accommodate right-censorship by reversing
the sequence of time-series and applying Algorithm 1 (resulting in our ability
to infer the degree of right censorship). When both left-censorship and right-
censorship are present, it corrects for left-censorship explicitly (using
$\delta$) while right-censorship is implicitly accounted for since we only
maximize the likelihood of data up to the point that we observe time-series.

##  4 Identifiability under a noiseless model

While SubLign presents a viable, practical model for clustering and aligning
censored time-series data, it is worth reflecting upon whether we can ever
identify subtype and alignment from observational data. In what follows, we
present theoretical conditions that show that there exists conditions under
which the problem we study is identifiable.

#### Identifiability

We assume distinct time stamps for the $M$ observations in $X_{i}$. The
generative process we assume for $Y_{i}$, conditional on $X_{i}$, is:

$\displaystyle\forall i=\\{1,\ldots,N\\},\;s_{i}\sim\text{Cat}(K),\;$
$\displaystyle\forall m\in\\{1,\ldots,M\\},\;d\in D_{i,m},$ $\displaystyle
y_{i,m}[d]=f(\kappa(x_{i,m}+\delta_{i};\theta^{P}[s_{i},d]))$ (4)

where $s_{i}\in\\{1,\ldots,K\\}$ is the subtype for time-series $i$, $D_{i,m}$
denotes the set of all observations at time-step $m$ for time-series $i$ where
$\forall i,m,|D_{i,m}|\leq D$, and $\delta_{i}$ is the delayed entry value.
The link function $f:\mathbb{R}\to\mathbb{R}$ has no parameters whereas
$\kappa:\mathbb{R}\to\mathbb{R}$ is an unknown polynomial function of degree
$P\in\mathbb{Z}^{+}$. We denote the parameters of $\kappa$, for each subtype
and dimension (e.g., biomarker), as $\theta^{P}\in\mathbb{R}^{K\times
D\times(P+1)}$. We denote $\theta^{P}[s_{i},d]$ as selecting the
$(s_{i},d)$-th vector of size $P+1$ from the tensor $\theta^{P}$. Similarly,
$\theta^{P}[s_{i}]$ selects the $s_{i}$-th matrix of size $X\times(P+1)$. We
define $\theta_{p}$ as the $p$-th coefficient of any polynomial function
parameter set $\theta$.

By construction, we have that
$s_{i}=s_{i^{\prime}}\iff\theta^{P}[s_{i}]=\theta^{P}[s_{i^{\prime}}]$, i.e.,
for two subtypes, the values of each time-series are described either by
$y=f(\kappa(x;\theta^{P}[s_{i}]))$ or
$y=f(\kappa(x;\theta^{P}[s_{i^{\prime}}]))$.

We begin with a set of assumptions for identifiability.

######  Assumption 1.

$f$ is invertible, and
$\kappa(x,\theta)=\theta_{0}+\sum_{p=1}^{P}\theta_{p}x^{p}$ describes a family
of polynomial functions in x with parameter $\theta$ and degree $P>0$. The
parameters of each subtype are unique.

######  Assumption 2.

$M\geq P+1$, i.e., for each patient time-series, there exists at least one of
the $D$ features where we observe at least $P+1$ values.

######  Assumption 3.

For each subtype $s_{k}$, there exists a time-series $i$ whose alignment
$\delta_{i}=0$, meaning no delayed entry.

######  Theorem 1.

Under assumptions 1, 2, 3 for the model in Equation 4, we can identify the
time-delays $\delta_{1},\ldots,\delta_{N}$. We can identify the polynomial
coefficients $\theta^{P}$ up-to a permutation of its rows and columns and the
identity of $s_{1},\ldots,s_{N}$ up-to a permutation over $K$ choices.

Proof sketch:  We defer the full proof for Theorem 1 to Appendix A but provide
the sketch here. Consider the case where we have a single biomarker for each
patient. The proof is constructive; first we transform the data using the
inverse of $f$ resulting in a set of data drawn from polynomial equations. The
polynomial coefficients may be estimated from the observed data; we can then
find the roots of these polynomials and pick the smallest root. These roots
exactly quantify the degree of delayed entry; i.e. they tell us how much each
polynomial has been shifted by. We can correct each time-series for this
shift, re-estimate polynomial coefficients from the shifted-polynomials and
cluster them to reveal the underlying subtype identity for each time-series
(and consequently each patient).

Remarks: Theorem 1 describes conditions under which delayed entry and the
polynomial parameters of cluster biomarker progression are identifiable. This
encouraging result demonstrates scenarios where the parameters of the model in
Equation 4 can provably be identified.

On the assumptions for identification: It is possible to relax Assumption 3 to
only require the existence of a single time-series from each subtype; this
modification only allows identifiability of $\delta_{1},\ldots,\delta_{N}$ up
to a translation within each subtype. The above result relies on the existence
of at least one biomarker for which there are sufficiently many observations –
this is a reasonable assumption in the context of clinical data since there is
often a _canonical_ biomarker tracked over time for each disease.

On the strategy for identification: We conjecture our analysis for
identification of subtype and alignment is of independent interest for
identifying causal effects in survival analysis where an important challenge
is how to handle confounding that jointly affects both survival time and
censorship. Related work in this field [Seaman et al., 2020, Choi and
O’Malley, 2017] has focused on restricting the class of models used to
characterize the survival function. Our work presents distinct parameteric
assumptions towards this goal.

##  5 Evaluation Setup

#### Synthetic data

We generate two classes of synthetic datasets from the _sigmoid_ and
_quadratic_ parameterizations in the ‘SubLign: Subtype & Align’ section. For
the sigmoid dataset, we generate data from $K=2$ subtypes and $\forall
i,m,|D_{i,m}|=3$ biomarker dimensions. For the quadratic dataset, we generate
data with $K=2$ and $\forall i,m,|D_{i,m}|=1$. See appendix for the data
generation process for the six quadratic datasets.

In both synthetic settings, we sample $N=1000$ patients with $M=4$
observations, variance $\sigma^{2}=0.25$, and max disease stage $T^{+}=10$.
For each patient, we sample subtype $s\sim\text{Bern}(0.5)$. The true disease
stage is drawn $t_{m}\sim\text{Unif}(0,T^{+})$ for observation
$m\in\\{1,\ldots,M\\}$. The biomarker values are drawn $y_{m}\sim
N(\lambda_{m},\sigma^{2})$ where
$\lambda_{m}=\sum_{k\in\\{1,\ldots,K\\}}\mathbbm{1}(s_{i}=k)f_{k}(t_{m})$. For
the sigmoid dataset, the first subtype generating function across three
dimensions is $f_{1}(t)=[\sigma(-4+t),\sigma(-1+t),\sigma(-8+8t)]$ and the
second subtype generating function is
$f_{2}(t)=[\sigma(-1+t),\sigma(-8+8t),\sigma(-25+3.5t)]$. The observed disease
time $x_{m}$ is shifted such that the first patient observation is at time 0.
Therefore $x_{m}=t_{m}-\zeta$ where $\zeta=\min_{j\in\\{1,\ldots,M\\}}t_{j}$
is the earliest true disease time for the patient.

#### Real-world clinical data

We study two real-world clinical datasets (Parkinson’s disease and heart
failure) to assess SubLign performance when data have realistic dynamics and
are missing over time. Our clinical datasets inherently contain delayed entry
due to external factors including patient access to healthcare. We include
preprocessing code and internal review board details in the appendix.

Heart failure (HF): We use electronic health records from a large health
system in the United States, currently redacted for double blind review. We
identify patients who enter the emergency department with a diagnosis of HF
and extract echocardiogram values—measurements from an ultrasound of the
heart—from the full patient history. We include echocardiogram features that
are present in more than 60% of echo studies. Our dataset includes $N=1534$
patients and $|D_{i,m}|\leq 12;\forall i,m$ features with $M=38$ maximum
observations per patient over a potential span of 10 years in the dataset. We
extract 27 baseline features including race, sex, and comorbidities (e.g.,
renal failure) to validate subtypes only and not for use in the model. The
dataset values are linearly scaled such that values are between 0 and 1 with
larger values denoting more abnormality.

Parkinson’s disease (PD):  We use publicly-available data from the Parkinson’s
Progression Markers Initiative (PPMI), an observational clinical study,
totalling $N_{t}=423$ PD patients and $N_{c}=196$ healthy controls where
$N=N_{t}+N_{c}$. We extract four biomarker measurements of autonomic, motor,
non-motor, and cognitive ability from $N=619$ total participants with $M=17$
maximum observations per patient. Our baseline data include $25$ features
including demographic information and patient history used to validate
subtypes. For PD patients, the first recorded visit is within 2 years of the
patient’s PD diagnosis. Measurements are scaled between 0 and 1 with larger
values corresponding to more abnormal values. We use the sigmoid
parameterization of the SubLign model for both datasets because HF and PD are
chronic and incurable diseases.

###  5.1 Hyperparameters and Baselines

We find optimal hyperparameters via grid search. For both synthetic and
clinical experiments, we search over hyperparameters including dimensions of
the latent space $z$ (2, 5, 10), the number of hidden units in the RNN (50,
100, 200), the number of hidden units in the multi-layer perceptron (50, 100,
200), the learning rate (0.001, 0.01, 0.1, 1.), regularization parameter (0.,
0.1, 1.), and regularization type (L1, L2). We select the hyperparameter
configuration with the best validation loss, as measured by Equation 3. Our
models are implemented in Python 3.7 using PyTorch [Paszke et al., 2019] and
are learned via Adam [Kingma and Ba, 2014] on a single NVIDIA k80 GPU for 1000
epochs. We set alignment extrema $\delta^{+}=10$ based on the maximum of the
synthetic dataset and the maxima of the HF and PD datasets. We search over 50
time steps with $\epsilon=0.1$.

We compare to seven different baselines. Our greedy baseline, denoted as
KMeans+Loss, first clusters the observed values using $k$-means clustering.
Then, using the inferred labels $s$, we simultaneously learn $\theta_{k}$ for
each subtype and $\delta_{i}$ for each patient by minimizing:

$\arg\min_{\theta,\delta}\sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{d=1}^{D}\sum_{k=1}^{K}\mathbbm{1}[s_{i}=k][y_{i,j,d}-f(x_{i,j,d}+\delta_{i};\theta_{k})]^{2}$
(5)

using Broyden–Fletcher–Goldfarb–Shanno. This naive clustering based approach
(in the space of the original data) attempts to correct for shifts in time. We
also compare to:

  1. 1.

SubNoLign: a modified SubLign with no alignment value. This model is
comparable to [Zhang et al., 2019], which learns a deep patient representation
while controlling for model architecture

  2. 2.

SuStaIn [Young et al., 2018]: a subtype and stage inference algorithm for
disease progression,

  3. 3.

BayLong: a Bayesian model of longitudinal clinical data [Huopaniemi et al.,
2014],

  4. 4.

PAGA [Wolf et al., 2019]: a state-of-the-art single-cell trajectory pseudo-
time method,

  5. 5.

Clustering with dynamic time warping (DTW) using kernel methods [Cuturi, 2011]
and soft-DTW [Cuturi and Blondel, 2017],

  6. 6.

SPARTan: A tensor factorization based approach for phenotyping from time-
series data [Perros et al., 2017]

We detail baseline implementations in the appendix.

###  5.2 Evaluation

We evaluate models on 5 trials, each with a different randomized data split
and random seed. For each trial, we learn on a training set (60%), find the
best performance across all hyperparameters on the validation set (20%), and
report the performance metrics on the held-out test set (20%). The same data
folds are used across all models for each trial.

We report the performance on the test set over three metrics. Adjusted Rand
index (ARI) measures whether pairs of samples are correctly assigned in the
same or different subtypes [Hubert and Arabie, 1985]. The Swaps metric reports
the number of swaps needed to sort the predicted disease times into the true
disease stages, expressed as percentage of total possible swaps, with full
equation in the appendix. The Pearson correlation coefficient expresses
correlation between the predicted and true disease stage. ARI measures the
clustering performance while the Swaps metric and the Pearson correlation
coefficient quantify how well the learning algorithm infers the alignment
values.

Model |  ARI $\uparrow$ |  Swaps $\downarrow$ |  Pearson $\uparrow$
---|---|---|---
SubLign | 0.94 $\pm$ 0.02 | 0.09 $\pm$ 0.00 | 0.85 $\pm$ 0.04
SubNoLign |  0.81 $\pm$ 0.21 | – | –
KMeans+Loss |  0.67 $\pm$ 0.04 |  0.21 $\pm$ 0.03 |  0.49 $\pm$ 0.01
SuStaIn |  0.66 $\pm$ 0.02 |  0.16 $\pm$ 0.00 |  0.30 $\pm$ 0.02
BayLong |  0.19 $\pm$ 0.18 |  0.48 $\pm$ 0.00 |  0.01 $\pm$ 0.02
PAGA |  0.32 $\pm$ 0.05 |  0.52 $\pm$ 0.07 |  0.04 $\pm$ 0.20
Soft-DTW |  0.06 $\pm$ 0.01 | – | –
Kernel-DTW |  0.06 $\pm$ 0.07 | – | –
SPARTan |  0.22 $\pm$ 0.18 | – | –
(a) width=0.8

Model | ARI
---|---
SubLign | 0.58 $\pm$ 0.12
SubNoLign |  0.42 $\pm$ 0.14
KMeans+Loss |  0.05 $\pm$ 0.04
SuStaIn |  0.12 $\pm$ 0.11
BayLong |  0.04 $\pm$ 0.17
PAGA |  0.02 $\pm$ 0.02
Soft-DTW |  0.46 $\pm$ 0.43
Kernel-DTW |  0.21 $\pm$ 0.36
SPARTan |  0.15 $\pm$ 0.10
(b)

Figure 2: Means and standard deviations over 5 trials for: (a) synthetic
sigmoid dataset with 1000 patients, 3 dimensions, and 4 observations per
patient, (b) 619 patients in the PPMI dataset including 423 Parkinson’s
disease patients and 196 healthy controls. Baseline methods include SuStaIn
[Young et al., 2018], BayLong [Huopaniemi et al., 2014], PAGA [Wolf et al.,
2019], Soft-DTW [Cuturi, 2011], Kernel-DTW [Dhillon et al., 2004], and SPARTan
[Perros et al., 2017].

#### Quantitative metrics on clinical datasets:

Because real world data often lack ground truth labels for subtype or
alignment, we create two semi-synthetic experiments with clinical datasets.
For HF, we evaluate SubLign’s ability to infer relative disease stage by
introducing additional censoring into the test sets. Specifically, we train
SubLign using 80% data (train and validation data) as usual. We then modify
the remaining data (20%) by removing the first year of patient observations,
creating distorted test set $(X^{\prime},Y^{\prime})$, and by removing the
last year of patient observations, creating
$(X^{\prime\prime},Y^{\prime\prime})$. The same amounts of observations are
removed from each set to control for length of observations. We infer
alignment values using the trained SubLign model: $\delta^{\prime}$ from
$(X^{\prime},Y^{\prime})$ and $\delta^{\prime\prime}$ from
$(X^{\prime\prime},Y^{\prime\prime})$. By construction,
$\delta^{\prime}>\delta^{\prime\prime}$. We report the percentage of patients
for which SubLign is able to recover this relationship. For PD, we report the
held-out clustering performance for healthy control patients and patients with
PD. We use disease status (PD patient or healthy control) as labels and $K=2$
subtypes.

##  6 Evaluation and Analysis

###  6.1 Recovering Subtypes with interval censoring

SubLign is able to recover subtypes despite interval censoring, outperforming
all baselines. For sigmoid synthetic data (Figure 2(a), ARI column), SubLign
can recover subtypes (mean ARI of 0.94) better than the KMeans+Loss baseline
(0.67) which assumes a greedy approach. Not correcting for alignment time
decreases the quality of inferred subtypes as can be seen in SubNoLign (0.81).

Similarly, SubLign recovers known subtypes in the PD dataset with
statistically significantly higher ARI over baselines (see Figure 2(b)). When
ARI performance intervals overlap, we use a t-test on pairwise differences
over trials to compute statistical significance. We specify the optimal
parameters for experiments and further describe statistical significance
methods in the appendix.

Some baselines appear to suffer because SubLign leverages the longitudinal
nature of patient data compared to the cross-sectional assumptions of PAGA and
SuStaIn. Other baselines have strong priors, i.e., [Huopaniemi et al., 2014],
which may explain its poor performance. Dynamic time warping methods appear to
perform poorly for datasets with few observations or high missingness rates.
Lastly, the tensor factorization method [Perros et al., 2017] fails, likely
because transforming data from continuous to discrete time results in a very
large and extremely sparse matrix factorization that is a tricky optimization
problem.

We include additional results in the appendix, including visualizations of the
SubLign subtypes compared to baselines, model misspecification analysis, and
experiments varying the level of missingness.

###  6.2 Recovering Known Alignment Values

For the synthetic sigmoid data, SubLign outperforms baselines in inferring
alignment values (Figure 2(a), Swaps and Pearson columns). SubLign recovers
alignment values better according to the Swaps metric (mean value of 0.09) and
the Pearson metric (mean value of 0.85) compared to the next best baselines of
KMeans+Loss and SuStaIn. Note that many baselines only recover subtypes and do
not learn patient alignment values. Although the real-world clinical datasets
do not contain true alignment values, we use the previously described HF setup
with an artificially censored test set. We find that SubLign predicts known
alignment relationships in an altered dataset. When evaluated on the
manipulated test data, SubLign recovers the constructed relationship of
$\delta^{\prime}>\delta^{\prime\prime}$ with a higher performance ($71\%\pm
2\%$) over five trials compared to K-Means+Loss ($57.8\pm 4\%$) and SuStaIn
($53.8\pm 3\%$). See appendix for full results.

###  6.3 Clinical Insights from Correcting for Misalignment

We validate SubLign subtypes learned from the HF and PD datasets using known
clinical findings. The appendix includes a table of baseline features, which
are not included as input to SubLign, with statistically significant
differences in subtypes: 7 features (out of 26) for PD and 11 features (out of
27) for HF.

#### Heart failure

Cardiologists classify patients into two groups based on ejection fraction: HF
with reduced ejection fraction (systolic HF) and HF with preserved ejection
fraction (diastolic HF). However, in our HF dataset, over 30% of patients
correspond to neither group based on clinical diagnosis. For SubLign, we set
$K=3$, one more than the number of known groups, to discover a potentially new
subtype. Without ground truth subtype labels, we observe that SubLign finds
systolic HF and diastolic HF as statistically significant baseline features.
Of the three subtypes, subtype C corresponds to systolic HF, and subtype A and
B correspond to diastolic HF, mirroring known clinical heterogeneity in
diastolic HF [Shah et al., 2015]. Of the two diastolic HF subtypes, subtype A
has a higher proportion of women while subtype B has a higher rate of obese
patients, both subgroups with documented heterogeneity in diastolic HF [Duca
et al., 2018, Tadic and Cuspidi, 2019]. In contrast, the subtypes found by the
KMeans+Loss baseline do not include known systolic HF and diastolic HF as
statistically significant features. See appendix for full results.

#### Parkinson’s disease

Biomarkers used to track PD are self-reported, which can be biased,
subjective, and noisy. Despite these obstacles, SubLign discovers subtypes
that match known clinical findings on two cohorts: 619 combined PD and healthy
control patients, and 423 PD patients only.

For PD and healthy control patients, we run SubLign with $K=2$ to uncover
characteristics of the two known groups. SubLign subtype A clearly corresponds
to healthy controls whereas subtype B designates PD patients. Statistically
significant baseline features include all components of the University of
Pennsylvania Smell Identification Test (UPSIT), which is a measure of smell
dysfunction and highly linked to PD [Haehner et al., 2009], and having a full
sibling or biological dad with PD, which aligns with research suggesting PD
may be hereditary [Klein and Westenberger, 2012]. For PD patients only, we set
$K=3$ in SubLign to explore and discover potential disease heterogeneity.
Statistically significant baseline features include race and gender, which
parallel recent clinical findings about heterogeneity in PD manifestation
[Taylor et al., 2007, Duca et al., 2018] and indicate new potential areas for
future work. See appendix for tables of statistically significant baseline
features stratified by the discovered subtypes for HF and PD.

##  7 Discussion

We study the task of clustering interval-censored time-series data. We present
our method, SubLign, to learn latent representations of disease progression
that correct for temporal misalignment in real-world observations and consider
conditions for identifiability of subtype and alignment values. Empirically,
our method outperforms seven baselines, and analysis of subtypes reveals
clinically plausible findings. Better modeling of disease heterogeneity
through alignment can help clinicians and scientists to better understand and
predict how chronic diseases with many subtypes may progress. We hope that our
model—in learning a continuous latent space to model heterogeneity—may be
applied to other domains where subtypes and temporal alignment are entangled,
for example gene expression analysis [Bar-Joseph et al., 2002] or cancer
pathways [Wu et al., 2015].

Our model introduces directions for future work. Practically, SubLign assumes
that $\delta$ and $z$ are marginally independent. Intuitively, this means
that, across all subtypes of a disease, the time at which the patient enters
the data is independent of any other factor. There are certainly cases where
this assumption may easily be violated, and therefore it remains an important
area for further improvement. Our results for theoretical and noiseless
identification do not naturally extend to the noisy setting since root
finding, without further assumptions, can be sensitive to noise in the
coefficients of the polynomials. Alternative strategies for identification
that generalize our result are fertile ground for future work.

## Acknowledgements

The authors thank Steven Horng for assistance with the heart failure dataset
and Christina Ji for assistance with the Parkinson’s disease dataset. The
authors thank Steven Horng and Rahul Deo for clinical guidance on heart
failure experiments. The authors thank Rebecca Boiarsky, Hunter Lang, and
Monica Agrawal for helpful comments. Data used in the preparation of this
article were obtained from the Parkinson’s Progression Markers Initiative
(PPMI) database. PPMI — a public-private partnership — is funded by The
Michael J. Fox Foundation for Parkinson’s Research and funding partners. This
work was supported by NSF CAREER award #1350965.

## References

  * [Alaa and van der Schaar, 2019] Alaa, A. M. and van der Schaar, M. (2019).  Attentive state-space modeling of disease progression.  In Advances in Neural Information Processing Systems, pages 11334–11344.
  * [Bar-Joseph et al., 2002] Bar-Joseph, Z., Gerber, G., Gifford, D. K., Jaakkola, T. S., and Simon, I. (2002).  A new approach to analyzing gene expression time series data.  In Proceedings of the sixth annual international conference on Computational biology, pages 39–48.
  * [Brandon et al., 2005] Brandon, D. T., Isaac, L. A., and LaVeist, T. A. (2005).  The legacy of tuskegee and trust in medical care: is tuskegee responsible for race differences in mistrust of medical care?  Journal of the National Medical Association, 97(7):951.
  * [Camargo et al., 2007] Camargo, S. J., Robertson, A. W., Gaffney, S. J., Smyth, P., and Ghil, M. (2007).  Cluster analysis of typhoon tracks. part i: General properties.  Journal of Climate, 20(14):3635–3653.
  * [Chan et al., 2006] Chan, L., Hart, L. G., and Goodman, D. C. (2006).  Geographic access to health care for rural medicare beneficiaries.  The Journal of Rural Health, 22(2):140–146.
  * [Choi and O’Malley, 2017] Choi, J. and O’Malley, A. J. (2017).  Estimating the causal effect of treatment in observational studies with survival time endpoints and unmeasured confounding.  Journal of the Royal Statistical Society. Series C, Applied statistics, 66(1):159.
  * [Cuturi, 2011] Cuturi, M. (2011).  Fast global alignment kernels.  In Proceedings of the 28th international conference on machine learning (ICML-11), pages 929–936.
  * [Cuturi and Blondel, 2017] Cuturi, M. and Blondel, M. (2017).  Soft-dtw: a differentiable loss function for time-series.  arXiv preprint arXiv:1703.01541.
  * [Dhillon et al., 2004] Dhillon, I. S., Guan, Y., and Kulis, B. (2004).  Kernel k-means: spectral clustering and normalized cuts.  In Proceedings of the tenth ACM SIGKDD international conference on Knowledge discovery and data mining, pages 551–556.
  * [Doshi-Velez et al., 2014] Doshi-Velez, F., Ge, Y., and Kohane, I. (2014).  Comorbidity clusters in autism spectrum disorders: An electronic health record time-series analysis.  Pediatrics, 133(1).
  * [Duca et al., 2018] Duca, F., Zotter-Tufaro, C., Kammerlander, A. A., Aschauer, S., Binder, C., Mascherbauer, J., and Bonderman, D. (2018).  Gender-related differences in heart failure with preserved ejection fraction.  Scientific reports, 8(1):1–9.
  * [Faucheux et al., 2021] Faucheux, L., Resche-Rigon, M., Curis, E., Soumelis, V., and Chevret, S. (2021).  Clustering with missing and left-censored data: A simulation study comparing multiple-imputation-based procedures.  Biometrical Journal, 63(2):372–393.
  * [Fereshtehnejad et al., 2017] Fereshtehnejad, S.-M., Zeighami, Y., Dagher, A., and Postuma, R. B. (2017).  Clinical criteria for subtyping parkinson’s disease: biomarkers and longitudinal progression.  Brain, 140(7):1959–1976.
  * [Gaffney and Smyth, 2005] Gaffney, S. J. and Smyth, P. (2005).  Joint probabilistic curve clustering and alignment.  In Advances in neural information processing systems, pages 473–480.
  * [Haehner et al., 2009] Haehner, A., Boesveldt, S., Berendse, H., Mackay-Sim, A., Fleischmann, J., Silburn, P., Johnston, A., Mellick, G., Herting, B., Reichmann, H., et al. (2009).  Prevalence of smell loss in parkinson’s disease–a multicenter study.  Parkinsonism & related disorders, 15(7):490–494.
  * [Hubert and Arabie, 1985] Hubert, L. and Arabie, P. (1985).  Comparing partitions.  Journal of classification, 2(1):193–218.
  * [Huopaniemi et al., 2014] Huopaniemi, I., Nadkarni, G., Nadukuru, R., Lotay, V., Ellis, S., Gottesman, O., and Bottinger, E. P. (2014).  Disease progression subtype discovery from longitudinal emr data with a majority of missing values and unknown initial time points.  In AMIA Annual Symposium Proceedings, volume 2014, page 709. American Medical Informatics Association.
  * [Kingma and Ba, 2014] Kingma, D. P. and Ba, J. (2014).  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980.
  * [Kingma and Welling, 2013] Kingma, D. P. and Welling, M. (2013).  Auto-encoding variational bayes.  The 2nd International Conference on Learning Representations (ICLR).
  * [Klein and Westenberger, 2012] Klein, C. and Westenberger, A. (2012).  Genetics of parkinson’s disease.  Cold Spring Harbor perspectives in medicine, 2(1):a008888.
  * [Li et al., 2011] Li, G., Liu, B., Qin, S. J., and Zhou, D. (2011).  Dynamic latent variable modeling for statistical process monitoring.  In Proc. IFAC World Congress, Milan, Italy, pages 12886–12891.
  * [Listgarten et al., 2007] Listgarten, J., Neal, R. M., Roweis, S. T., Puckrin, R., and Cutler, S. (2007).  Bayesian detection of infrequent differences in sets of time series with shared structure.  In Advances in neural information processing systems, pages 905–912.
  * [Liu et al., 2009] Liu, X., Tong, Y., and Wheeler, F. W. (2009).  Simultaneous alignment and clustering for an image ensemble.  In 2009 IEEE 12th International Conference on Computer Vision, pages 1327–1334. IEEE.
  * [Luo et al., 2020] Luo, Y., Eran, A., Palmer, N., Avillach, P., Levy-Moonshine, A., Szolovits, P., and Kohane, I. S. (2020).  A multidimensional precision medicine approach identifies an autism subtype characterized by dyslipidemia.  Nature Medicine, 26(9):1375–1379.
  * [Mattar et al., 2012] Mattar, M. A., Hanson, A. R., and Learned-Miller, E. G. (2012).  Unsupervised joint alignment and clustering using bayesian nonparametrics.  In Proceedings of the Twenty-Eighth Conference on Uncertainty in Artificial Intelligence, pages 584–593.
  * [Miller and Wherry, 2017] Miller, S. and Wherry, L. R. (2017).  Health and access to care during the first 2 years of the aca medicaid expansions.  New England Journal of Medicine, 376(10):947–956.
  * [Miller Jr, 2011] Miller Jr, R. G. (2011).  Survival analysis, volume 66.  John Wiley & Sons.
  * [Owan et al., 2006] Owan, T. E., Hodge, D. O., Herges, R. M., Jacobsen, S. J., Roger, V. L., and Redfield, M. M. (2006).  Trends in prevalence and outcome of heart failure with preserved ejection fraction.  New England Journal of Medicine, 355(3):251–259.
  * [Paszke et al., 2019] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al. (2019).  Pytorch: An imperative style, high-performance deep learning library.  In Advances in neural information processing systems, pages 8026–8037.
  * [Perros et al., 2017] Perros, I., Papalexakis, E. E., Wang, F., Vuduc, R., Searles, E., Thompson, M., and Sun, J. (2017).  Spartan: Scalable parafac2 for large & sparse data.  In Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 375–384.
  * [Rebbapragada et al., 2009] Rebbapragada, U., Protopapas, P., Brodley, C. E., and Alcock, C. (2009).  Finding anomalous periodic time series.  Machine learning, 74(3):281–313.
  * [Rezende et al., 2014] Rezende, D. J., Mohamed, S., and Wierstra, D. (2014).  Stochastic backpropagation and approximate inference in deep generative models.  In Proceedings of the 31st International Conference on International Conference on Machine Learning-Volume 32, pages II–1278.
  * [Seaman et al., 2020] Seaman, S., Dukes, O., Keogh, R., and Vansteelandt, S. (2020).  Adjusting for time-varying confounders in survival analysis using structural nested cumulative survival time models.  Biometrics, 76(2):472–483.
  * [Shah et al., 2015] Shah, S. J., Katz, D. H., Selvaraj, S., Burke, M. A., Yancy, C. W., Gheorghiade, M., Bonow, R. O., Huang, C.-C., and Deo, R. C. (2015).  Phenomapping for novel classification of heart failure with preserved ejection fraction.  Circulation, 131(3):269–279.
  * [Tadic and Cuspidi, 2019] Tadic, M. and Cuspidi, C. (2019).  Obesity and heart failure with preserved ejection fraction: a paradox or something else?  Heart Failure Reviews, 24(3):379–385.
  * [Taylor et al., 2007] Taylor, K. S. M., Cook, J. A., and Counsell, C. E. (2007).  Heterogeneity in male to female risk for parkinson’s disease.  Journal of Neurology, Neurosurgery & Psychiatry, 78(8):905–906.
  * [Udler et al., 2018] Udler, M. S., Kim, J., von Grotthuss, M., Bonás-Guarch, S., Cole, J. B., Chiou, J., on behalf of METASTROKE, C. D. A., the ISGC, Boehnke, M., Laakso, M., Atzmon, G., et al. (2018).  Type 2 diabetes genetic loci informed by multi-trait associations point to disease mechanisms and subtypes: a soft clustering analysis.  PLoS medicine, 15(9):e1002654.
  * [White et al., 2011] White, I. R., Royston, P., and Wood, A. M. (2011).  Multiple imputation using chained equations: issues and guidance for practice.  Statistics in medicine, 30(4):377–399.
  * [Wolf et al., 2019] Wolf, F. A., Hamey, F. K., Plass, M., Solana, J., Dahlin, J. S., Göttgens, B., Rajewsky, N., Simon, L., and Theis, F. J. (2019).  Paga: graph abstraction reconciles clustering with trajectory inference through a topology preserving map of single cells.  Genome biology, 20(1):59.
  * [Wu et al., 2015] Wu, H., Gao, L., and Kasabov, N. (2015).  Inference of cancer progression from somatic mutation data.  IFAC-PapersOnLine, 48(28):234–238.
  * [Yoon et al., 2018] Yoon, J., Zame, W. R., and van der Schaar, M. (2018).  Estimating missing data in temporal data streams using multi-directional recurrent neural networks.  IEEE Transactions on Biomedical Engineering, 66(5):1477–1490.
  * [Young et al., 2018] Young, A. L., Marinescu, R. V., Oxtoby, N. P., Bocchetta, M., Yong, K., Firth, N. C., Cash, D. M., Thomas, D. L., Dick, K. M., Cardoso, J., et al. (2018).  Uncovering the heterogeneity and temporal complexity of neurodegenerative diseases with subtype and stage inference.  Nature communications, 9(1):4273.
  * [Zhang et al., 2019] Zhang, X., Chou, J., Liang, J., Xiao, C., Zhao, Y., Sarva, H., Henchcliffe, C., and Wang, F. (2019).  Data-driven subtyping of parkinson’s disease using longitudinal clinical records: a cohort study.  Scientific reports, 9(1):1–12.

## Appendix

  * •

In Section A, we present the proof for conditions of identifiability and
variational lower bound derivation.

  * •

In Section B, we describe additional experiment setup details for the sigmoid,
quadratic, and clinical experiments.

  * •

In Section C, we present additional experimental results on discovering
subtypes including visualization of subtypes, model misspecification, and
missingness methods.

    * –

Visualizations of subtypes show that SubLign subtypes match the data
generating function closer than SubNoLign.

    * –

Model misspecification results show that SubLign is robust to misspecification
from cubic spline data generation.

    * –

Missingness method experiments show that SubLign outperforms baselines for
various missingness rates with comparable performance to some baselines for no
missingness.

  * •

In Section D, we present quadratic experiment results.

    * –

An additional six synthetic experiments show that SubLign outperforms
baselines.

    * –

In cases where alignment values are not identifiable (e.g., one subtype’s data
generating function is a flat line), SubLign cannot recover those alignment
values.

  * •

In Section E, we discuss internal review board and patient privacy details for
the clinical datasets used.

  * •

In Section F, we present the set of statistically significant features for the
heart failure and Parkinson’s disease clinical experiments.

    * –

Model misspecification results show that SubLign is robust to misspecification
from cubic spline data generation.

    * –

Missingness method experiments show that SubLign outperforms baselines for
various missingness rates with comparable performance to some baselines for no
missingness.

##  Appendix A Identifiability and Inference

###  A.1 Identifiability

Algorithm 2 Procedure for the identification of model parameters

1: Input: Observation times $X\in\mathbb{R}^{N\times M}$, biomarkers
$Y\in\mathbb{R}^{N\times M\times D}$, polynomial degree $P$, invertible
function $f$

2: Output: $\theta^{P}$, $\delta_{1},\ldots,\delta_{N}$, $s_{1},\ldots,s_{N}$
for each patient

3: Step 1: Transform the observed biomarkers; $Q=f^{-1}(Y)$

4: Step 2: Obtain time-shifts using a single biomarker;

5: a) For each patient $i$, estimate the parameters $\dot{\theta^{1}}_{i}$ of
$\kappa(x;\dot{\theta^{1}}_{i})$ using a single biomarker
$((x_{i,1},q_{i,1}),\ldots,(x_{i,M},q_{i,M})$ via polynomial regression,

6: b) Compute up to $P$ roots of polynomial $\kappa(x,\dot{\theta^{1}}_{i})$
for each patient $i$ as $R_{i}=\\{r_{1},\ldots,r_{P}\\}$ and set
$\xi_{i}=\min\text{Real}(R_{i})$ where Real denotes the real part of
(potentially complex) roots.

7: c) Estimate $\widetilde{\theta^{1}}_{i}$ for polynomials in a _canonical
position_ using $((x_{i,1}-\xi_{i},q_{i,1}),\ldots,(x_{i,M}-\xi_{i},q_{i,M})$
via polynomial regression,

8: d) Cluster $\widetilde{\theta^{j}}_{i}$ across patients via K-means
clustering to yield cluster identities $s_{1},\ldots,s_{N}$

9: e) $\forall k,\;\;\eta_{k}=\min\\{\xi_{i}\,|\,i\text{ s.t. }s_{i}=k\\}$ and
$\forall i,\;\;\delta_{i}=\xi_{i}-\eta_{s_{i}}$

10: Step 3: Estimate true polynomial coefficients using shifted observation
times;

11: for biomarker $j=1$ to $J$ do

12: For each patient, estimate the parameters $\theta^{j}_{i}$ of
$\kappa(x;\hat{\theta^{j}}_{i})$ using
$((x_{1,1}-\delta_{i},q_{1,1}[j]),\ldots,(x_{1,M}-\delta_{i},q_{1,M}[j])$ via
polynomial regression,

13: end for

14: Return
$\theta^{P}=[\theta^{1}|\ldots|\theta^{J}],\\{\delta_{1},\ldots,\delta_{N}\\},\\{s_{1},\ldots,s_{N}\\}$

We restate our assumptions.

######  Assumption 1.

$f$ is invertible, and
$\kappa(x,\theta)=\theta_{0}+\sum_{p=1}^{P}\theta_{p}x^{p}$ describes a family
of polynomial functions in x with parameter $\theta$ and degree $P>0$. The
parameters of each subtype are unique.

######  Assumption 2.

$M\geq P+1$, i.e., for each object, across all the $D$ features, we observe at
least $P+1$ values.

######  Assumption 3.

For each subtype $s_{k}$, there exists an object $i$ whose alignment
$\delta_{i}=0$.

We provide the proof for Theorem 1 below:

###### Proof.

The proof is constructive; i.e. we give an algorithm for the identification of
the parameters of the model in Equation 4. The algorithm for identification is
presented in Algorithm 2 and proceeds in three steps.

Step 1: The first step transforms the observed biomarkers by applying the
inverse of function $f$, which exists by Assumption 1. This leaves us with
data as:

$f^{-1}(y_{i,m})=\kappa(x_{i,m}+\delta_{i};\theta^{P}_{s_{i}})\;\;\forall i\in
N,m\in M$

i.e. for all bio-markers, across all patients, we have data arising from
different polynomial functions.

Step 2: Without loss of generality, the second step uses the first biomarker
to identify the values of $\delta_{i}$ for each patient.

  1. a)

First, we estimate the polynomial coefficients for each patient separately; we
are guaranteed exact recovery of the coefficients by Assumption 2.

  2. b)

Next we find the roots for each polynomial. If they are complex, consider
their real part, and define $\xi_{i}$ to be the smallest root of the
polynomial. At least one (real or complex) root is guaranteed to exist by the
Fundamental Theorem of Algebra for every non-constant polynomial (Assumption
1). Note that the choice of using the smallest root is arbitrary; what matters
is that a consistent choice of root is selected for each patient’s
polynomials.

  3. c)

The goal of this step to learn a new polynomial for each patient which is
translated to ensure that the root selected in step b) lies at $x=0$.

To do so, we first shift the observational time-steps by $\xi_{i}$, and we re-
estimate the coefficients of each _shifted_ polynomial.

We make use of the fact that if $\xi_{i}$ is the smallest complex root of a
polynomial $\kappa(x)$ then the polynomial $\kappa(x+\xi_{i})$ has its
smallest complex root at $0$. We can recover the parameters of this polynomial
exactly by shifting our observations and re-estimating the coefficients.

This operation recovers the coefficients of every patient’s polynomial in its
_canonical position_ i.e. a translated polynomial whose the smallest root (or
its real component) is at $x=0$.

This step can be viewed as a de-biasing step which allows us to re-estimate
$\widetilde{\theta}$ without while ignoring the effect that left-censorship
has on parameter estimates.

  4. d)

We cluster the coefficients estimated in step c). By construction, we know
that $s_{i}=s_{i^{\prime}}\iff\theta_{i}=\theta_{i^{\prime}}$ which guarantees
that clustering recovers the true-underlying subtype for each patient (up to a
permutation over $K$ choices).

  5. e)

Finally we stratify patients by their subtype, and we define $\delta_{i}$ as
the difference between their smallest root and the smallest value of $\xi_{i}$
among all other patients within that subtype.

By Assumption 3, we know that for each subtype, there exists a patient for
whom $\delta_{i}=0$, this reference patient will also be the one whose
polynomial has the smallest root. We note here that without Assumption 3, we
would still have identification of $\delta_{i}$ up to a constant.

Therefore, by shifting each patient’s smallest root by their reference
patient’s smallest root, we can recover the original time-shifts.

Step 3: Given the values of $\delta_{1},\ldots,\delta_{N}$ from Step 2, we can
now estimate the true values of the polynomial coefficients exactly in the
noiseless setting via polynomial regression. ∎

###  A.2 Variational Lower Bound

$\displaystyle\log$ $\displaystyle\;p(Y|X;\gamma)$
$\displaystyle=\log\int_{Z,\delta}p(Y,Z,\delta|X;\gamma)\,dZ\,d\delta$
$\displaystyle=\log\int_{Z,\delta}q(Z|X,Y;\phi)\frac{p(Y,Z,\delta|X;\gamma)}{q(Z|X,Y;\phi)}\,dZ\,d\delta$
$\displaystyle\geq\int_{Z,\delta}q(Z|X,Y;\phi)\log\frac{p(Y,Z,\delta|X;\gamma)}{q(Z|X,Y;\phi)}\,dZ\,d\delta$
(6)
$\displaystyle=\int_{\delta}q(\delta)\int_{Z}q(Z|X,Y;\phi)\log\frac{p(\delta)p(Y,Z|X,\delta;\gamma)}{q(\delta)q(Z|X,Y;\phi)}$
$\displaystyle=\int_{\delta}q(\delta)\int_{Z}q(Z|X,Y;\phi)\log\frac{p(Y,Z|X,\delta;\gamma)}{q(Z|X,Y;\phi)}+\int_{\delta}q(\delta)\int_{Z}q(Z|X,Y;\phi)\log\frac{p(\delta)}{q(\delta)}$
$\displaystyle=\mathbb{E}_{q(Z|X,Y;\phi)}\left[\log\frac{p(Y,Z,\delta|X,\delta;\gamma)}{q(Z|X,Y;\phi)}\right]$
(7)

##  Appendix B Experiment Setup Details

Here we describe additional details about experiment setup:

  * •

In Section B.1, we describe the swaps metric.

  * •

In Section B.2, we describe the specific parameters used for the baselines.

  * •

In Section B.3, we specify the optimal hyperparameters chosen for sigmoid
experiments and clinical experiments.

  * •

In Section B.4, we detail the biomarkers and baseline features used for both
clinical datasets: heart failure and Parkinson’s disease.

  * •

In Section B.5, we describe how SubLign handles missing values, including
irregularly lengths of time-series data.

  * •

In Section B.6, we describe how we compute statistical significance for
results.

###  B.1 Swaps Metric

For true sorted alignment values $a_{1},\ldots,a_{N}$ for $N$ patients, we
define the swaps metric $\mathcal{S}$ of proposed alignment values
$b_{1},\ldots,b_{N}$ as the number of swaps needed to sort the predicted
disease times into the true disease stages, expressed as percentage of total
possible swaps.

$\mathcal{S}=\frac{\sum_{i,j;i<j}\mathbbm{1}(a_{i}<b_{i},a_{j}>b_{j})}{N(N-1)/2}$

###  B.2 Baselines

We outline SubNoLign and KMeans+Loss in the main paper. Here we describe the
remaining baselines.

#### SuStaIn

SuStaIn [Young et al., 2018] is a disease progression algorithm that recovers
subtype and stage from cross-sectional data. We transform our longitudinal
data by dropping patient affiliation across visits. We transform the data by
subtracting the mean for each feature and dividing by the standard deviation
for each feature. We assume the Z-scored values have a max of 5. We run for
1,000,000 epochs for the Markov Chain Monte Carlo sampling and 1,000 epochs
for optimization. We use an open source implementation by the
authors.111https://github.com/ucl-pond/pySuStaIn

#### Bayesian approach

The Bayesian approach [Huopaniemi et al., 2014] assumes longitudinal data, but
there must be a small number of measured time points. We assume that there are
10 observed time points where observed data can begin as well as a window of
10 time points before the observed window where a patient’s values can be
aligned to. Because biomarker values are scaled between 0 and 10, we assume
that values change between time points based on a Gaussian with $\sigma=2$ and
that subtype means for each time point are drawn from a Gaussian with
$\sigma=5$. We draw 4000 samples and use the maximum a posteriori estimate to
determine stage and subtype for test patients. Because we could not find an
open-source option, we implemented the algorithm ourselves based on the
description in the paper.

#### PAGA

Partition-based graph abstraction, or PAGA, [Wolf et al., 2019] assumes cross-
sectional data, so we create separate visits for each patient visit. For
algorithm parameters, we set resolution to 0.05, number of neighbors to 15,
and connectivity cutoff of 0.05. We use an open source implementation by the
authors.222https://github.com/dynverse/ti_paga/blob/master/run.py

#### Dynamic time warping

Dynamic time warping (DTW) defines similarity between time series that can be
combined with clustering techniques. DTW methods include using soft-DTW
[Cuturi and Blondel, 2017] and kernel [Cuturi, 2011] before using K-means with
the chosen similarity metric. We use open source implementations of DTW
algorithms333https://pypi.org/project/dtw-python/ to generate our baseline
comparisons.

#### Tensor factorization

Sparse tensor factorization has been used for disease phenotyping. The
decomposition of large and sparse datasets using canonical polyadic
decomposition can create an interpretable output for phenotyping. We use the
Matlab open source implementation of SPARTan.
444https://github.com/kperros/SPARTan We found these baseline results to yield
poor clustering performance despite aggressive hyperparameter tuning. We
surmise this is because transforming our data from continuous to discrete time
resulted in a very large and extremely sparse matrix factorizing which is a
tricky optimization problem.

###  B.3 Optimal hyperparameters for sigmoid and clinical baselines

For all models, we run for 1000 epochs and use the model with the best
training loss over the 1000 epochs for evaluation. For the sigmoid dataset,
the optimal hyperparameters are latent space of dimension 5, 100 hidden units
in the RNN, 50 hidden units in the multi-layer perceptron, learning rate of
0.01, and no regularization.

For the Parkinson’s disease dataset, we searched on a slightly smaller set of
hyperparameters for SubLign and found optimal hyperparameters of $\beta=0.01$,
no regularization, 10 latent dimensions, 10 hidden units for the multi-layer
perceptron, 200 units for the recurrent neural network, and learning rate of
0.1.

For the heart failure dataset, we searched on a slightly smaller set of
hyperparameters for SubLign and found optimal hyperparameters of
$\beta=0.001$, no regularization, 10 latent dimensions, 20 hidden units for
the multi-layer perceptron, 50 units for the recurrent neural network, and
learning rate of 0.01.

###  B.4 Clinical dataset biomarkers and baseline features

The source of the heart failure dataset will be revealed after blind peer
review. For the heart failure dataset, we include the following biomarkers:
Aorta - Ascending, Aorta - Valve Level, Aortic Valve - Peak Velocity, Left
Atrium - Four Chamber Length, Left Atrium - Long Axis Dimension, Left
Ventricle - Diastolic Dimension, Left Ventricle - Ejection Fraction, Left
Ventricle - Inferolateral Thickness, Left Ventricle - Septal Wall Thickness,
Mitral Valve - E Wave, Mitral Valve - E Wave Deceleration Time, and Right
Atrium - Four Chamber Length.

From the Parkinson’s Progression Markers Initiative (PPMI) dataset, we include
four main biomarkers: 1) MOCA, a cognitive assessment, 2) SCOPA-AUT, an
autonomic assessment, 3) NUPDRS1, an assessment of non-motor symptoms, and 4)
a maximum taken over NUPDRS3 and NUPDRS2 as an assessment of motor symptoms.
We removed patients without extractable biomarker measurements and include
preprocessing code.

The baseline features considered for heart failure are: age, anemia,
atherosclerosis, atrial fibrillation, Black, body mass index, chronic kidney
disease, diastolic heart failure, esophageal reflux, female, hyperlipidemia,
hypertension, hypothyroidism, kidney disease, major depressive disorder,
obesity old myocardial infarction, other race, pulmonary heart disease,
pneumonia, renal failure, type 2 diabetes, urinary tract infection, and White.

The baseline features considered for Parkinson’s disease (PD) are: male,
Hispanic/Latino, White, Asian, Black, American Indian, Pacific Islander, not
specified race, biological mom with PD, biological dad with PD, full sibling
with PD, half sibling with PD, maternal grandparent with PD, paternal
grandparent with PD, maternal aunt/uncle with PD, paternal aunt/uncle with PD,
kids with PD, years of education, right handed, left handed, University of
Pennsylvania Smell Identification Test (UPSIT) part 1, UPSIT part 2, UPSIT
part 3, UPSIT part 4, and UPSIT total.

###  B.5 Missing values

SubLign allows for missing biomarker dimensions and missing patient visits to
accommodate the sparsity of clinical data. For missing visits, we adapt the
recognition network to handle variable sequence lengths. We mask out missing
observations so they have no contribution to the learning stage, except for
the recognition network input. For the recognition network input, we linearly
interpolate missing values for each patient. For baselines that cannot handle
missing data, we also linearly interpolate missing values for each patient.

###  B.6 Statistical significance

To estimate robustness of our models, we evaluate our held-out performance
over 5 trials. Each trial consists of randomized 60/20/20
training/validation/test data folds and a different random seed. In order to
compare models across 5 trials, we report the means and standard deviations
from the 5 trials. When the reported performance intervals overlap, we compute
the statistical significance of the pairwise differences using a t-test and a
Benjamini-Hochberg correction of 0.05.

##  Appendix C Additional Experiment Results

  * •

In Section C.1, we visualize the SubLign and SubNoLign subtypes for sigmoid
data.

  * •

In Section C.2, we present results on model misspecification.

  * •

In Section C.3, we present the empirical results with varying levels of
missingness.

###  C.1 SubLign and SubNoLign subtype visualization

In Figure 3, we show the visualization for SubLign subtypes compared to
SubNoLign for the first dimension of the sigmoid dataset. We find that the
visualization of the SubNoLign subtypes are not as close to the data
generating function as the SubLign subtypes. All other parameters, data
dimensions, and experimental conditions are held constant.

Figure 3: One of three dimensions of learned SubLign and SubNoLign subtypes
from sigmoid synthetic data plotted on top of original data generating
functions.

###  C.2 Model Misspecification

Because our model implicitly assumes a functional form (e.g., sigmoid or
quadratic), we investigate learning under model misspecification. With
synthetic datasets created using splines on 5 randomly generated control
points with the same noise rates, dimensions, and censoring as the original
experiments, we create two settings: one where the control points are
monotonically increasing (using SubLign with a sigmoid $f$ function) and one
with no restrictions (using SubLign with a quadratic $f$ function).

We run SubLign and use a sigmoid $f$ function the monotonically increasing
control points (labeled “Incr") and use quadratic $f$ function for splines
generated without restrictions (labeled “Any").

In Table 1, we present results on model misspecification over 2 settings and 5
trials of SubLign learning on data generated from splines with 5 control
points over 3 dimensions. We conclude that SubLign is robust against
reasonable model misspecification using piecewise polynomial data.

Model | ARI$\uparrow$ | Pearson $\uparrow$ | Swaps $\downarrow$
---|---|---|---
SubLign-Incr | 0.82 $\pm$ 0.17 | 0.83 $\pm$ 0.08 | 0.14 $\pm$ 0.04
SubNoLign-Incr | 0.77 $\pm$ 0.10 | – | –
KMeans+Loss-Incr | 0.58 $\pm$ 0.11 | 0.43 $\pm$ 0.09 | 0.21 $\pm$ 0.06
SubLign-Any | 0.46 $\pm$ 0.12 | 0.67 $\pm$ 0.39 | 0.22 $\pm$ 0.14
SubNoLign-Any | 0.29 $\pm$ 0.10 | – | –
KMeans+Loss-Any | 0.22 $\pm$ 0.07 | 0.23 $\pm$ 0.21 | 0.48 $\pm$ 0.11
Table 1: Model misspecification experiment means and standard deviations using
5 cubic splines datasets.

###  C.3 Missingness Experiments

Here we present experiments varying the amount of missingness in synthetic
datasets. Following the Parkinson’s disease dataset (PPMI) where 47-60% of the
biomarkers are missing with a maximum of 17 observations for each patient, we
modify our synthetic sigmoid dataset and remove biomarkers uniformly randomly
with different missingness rates and $M=17$.

Additionally we experiment with two missingness imputation methods for the
baseline models. In addition to the linear interpolation used in the results
presented in the main paper, we experiment with two other imputation methods:
chained equations (MICE) [White et al., 2011] and a multi-directional
recurrent neural network (MRNN) designed for multivariate time-series [Yoon et
al., 2018].

We find that SubLign outperforms baselines across higher missingness rates
although with no missing observations, some baselines are comparable. In Table
2, we show results with 50% of the data missing uniformly random. In Table 3,
we show results with 25% of the data missing uniformly random. In Table 4, we
show results with none of the data missing uniformly random.

Table 2: Experiments on synthetic data with 50% of the data missing. Baselines include SuStaIn [Young et al., 2018], BayLong [Huopaniemi et al., 2014], PAGA [Wolf et al., 2019], SPARTan [Perros et al., 2017], clustering using Soft-DTW [Cuturi, 2011], and clustering using Kernel-DTW [Dhillon et al., 2004]. Imputation methods include MICE [White et al., 2011] and MRNN [Yoon et al., 2018]. Imputation Method | Model | ARI $\uparrow$ | Swaps $\downarrow$ | Pearson $\uparrow$
---|---|---|---|---
– | SubLign | 0.813 $\pm$ 0.024 | 0.299 $\pm$ 0.019 | 0.613 $\pm$ 0.049
– | SubNoLign | 0.789 $\pm$ 0.058 | – | –
MICE | KMeans+Loss | 0.780 $\pm$ 0.046 | 0.327 $\pm$ 0.048 | 0.503 $\pm$ 0.018
| SuStaIn | 0.459 $\pm$ 0.010 | 0.243 $\pm$ 0.004 | 0.160 $\pm$ 0.030
| BayLong | 0.028 $\pm$ 0.003 | 0.480 $\pm$ 0.002 | 0.009 $\pm$ 0.003
| PAGA | 0.003 $\pm$ 0.002 | 0.494$\pm$ 0.028 | 0.034 $\pm$ 0.001
– | Soft-DTW | 0.081 $\pm$ 0.004 | – | –
– | Kernel-DTW | 0.013 $\pm$ 0.002 | – | –
| SPARTan | 0.081 $\pm$ 0.013 | – | –
MRNN | KMeans+Loss | 0.783 $\pm$ 0.071 | 0.321 $\pm$ 0.120 | 0.562 $\pm$ 0.042
| SuStaIn | 0.450 $\pm$ 0.120 | 0.304 $\pm$ 0.120 | 0.434 $\pm$ 0.120
| BayLong | 0.028 $\pm$ 0.003 | 0.480 $\pm$ 0.002 | 0.009 $\pm$ 0.003
| PAGA | 0.004 $\pm$ 0.001 | 0.492$\pm$ 0.031 | 0.032 $\pm$ 0.003
| Soft-DTW | 0.094 $\pm$ 0.005 | – | –
| Kernel-DTW | 0.000 $\pm$ 0.003 | – | –
| SPARTan | 0.091 $\pm$ 0.005 | – | –
Table 3: Experiments on synthetic data with 25% of the data missing. Baselines include SuStaIn [Young et al., 2018], BayLong [Huopaniemi et al., 2014], PAGA [Wolf et al., 2019], SPARTan [Perros et al., 2017], clustering using Soft-DTW [Cuturi, 2011], and clustering using Kernel-DTW [Dhillon et al., 2004]. Imputation methods include MICE [White et al., 2011] and MRNN [Yoon et al., 2018]. Imputation Method | Model | ARI $\uparrow$ | Swaps $\downarrow$ | Pearson $\uparrow$
---|---|---|---|---
– | SubLign | 0.811 $\pm$ 0.406 | 0.309 $\pm$ 0.028 | 0.554 $\pm$ 0.070
– | SubNoLign | 0.743 $\pm$ 0.040 | – | –
MICE | KMeans+Loss | 0.741 $\pm$ 0.041 | 0.332 $\pm$ 0.048 | 0.508 $\pm$ 0.014
| SuStaIn | 0.611 $\pm$ 0.406 | 0.309 $\pm$ 0.028 | 0.554 $\pm$ 0.070
| BayLong | 0.111 $\pm$ 0.010 | 0.464 $\pm$ 0.071 | 0.011 $\pm$ 0.058
| PAGA | 0.010 $\pm$ 0.003 | 0.433 $\pm$ 0.007 | 0.048 $\pm$ 0.030
| Soft-DTW | 0.151 $\pm$ 0.031 | – | –
| Kernel-DTW | 0.002 $\pm$ 0.004 | – | –
| SPARTan | 0.168 $\pm$ 0.008 | – | –
MRNN | KMeans+Loss | 0.653 $\pm$ 0.029 | 0.308 $\pm$ 0.358 | 0.497 $\pm$ 0.025
| SuStaIn | 0.615 $\pm$ 0.112 | 0.244 $\pm$ 0.004 | 0.577 $\pm$ 0.020
| BayLong | 0.108 $\pm$ 0.016 | 0.461 $\pm$ 0.005 | 0.010 $\pm$ 0.041
| PAGA | 0.011 $\pm$ 0.002 | 0.413 $\pm$ 0.005 | 0.031 $\pm$ 0.002
| Soft-DTW | 0.103 $\pm$ 0.011 | – | –
| Kernel-DTW | 0.006 $\pm$ 0.004 | – | –
| SPARTan | 0.171 $\pm$ 0.041 | – | –
Table 4: Experiments on synthetic data with 0% of the data missing. Model | ARI $\uparrow$ | Swaps $\downarrow$ | Pearson $\uparrow$
---|---|---|---
SubLign | 0.980 $\pm$ 0.000 | 0.273 $\pm$ 0.012 | 0.714 $\pm$ 0.022
SubNoLign | 0.809 $\pm$ 0.382 | – | –
KMeans+Loss | 0.980 $\pm$ 0.016 | 0.057 $\pm$ 0.112 | 0.480 $\pm$ 0.039
SuStaIn [Young et al., 2018] | 0.765 $\pm$ 0.012 | 0.144 $\pm$ 0.004 | 0.477 $\pm$ 0.020
BayLong [Huopaniemi et al., 2014] | 0.201 $\pm$ 0.166 | 0.451 $\pm$ 0.050 | 0.011 $\pm$ 0.021
PAGA [Wolf et al., 2019] | 0.251 $\pm$ 0.031 | 0.481 $\pm$ 0.003 | 0.015 $\pm$ 0.002
Soft-DTW [Cuturi, 2011] | 0.974 $\pm$ 0.015 | – | –
Kernel-DTW [Dhillon et al., 2004] | 0.949 $\pm$ 0.040 | – | –
SPARTan [Perros et al., 2017] | 0.251 $\pm$ 0.031 | – | –

##  Appendix D Quadratic Data Results

We describe an additional set of experiments using the quadratic functional
family. These experiments were designed to better understand where SubLign is
able to learn clustering and alignment metrics well.

  * •

In Section D.1, we detail the dataset creation.

  * •

In Section D.2, we outline the optimal hyperparameters for the quadratic
experiments.

  * •

In Section D.3, we describe the empirical results for the quadratic
experiments.

###  D.1 Setup

For the quadratic dataset, we generate data from 2 subtypes and 1 dimension
with generating functions. See Table 5 for subtype generating functions.
Similar to the sigmoid synthetic dataset, for each patient in the datasets, we
draw subtype $\kappa\sim Bern(0.5)$. The true disease stage is drawn
$t_{m}\sim U(0,T^{+})$ for observation $m\in[M]$. The biomarker values are
drawn $y_{m}\sim N(\mu_{m},S)$ where
$\mu_{m}=\sum_{k\in[K]}\mathbbm{1}(\kappa=k)f_{k}(t_{m})$. The observed
disease time $x_{m}$ is shifted such that the first patient observation is at
time 0. Therefore $x_{m}=t_{m}-\tau$ where $\tau=\min_{j\in[M]}t_{j}$ and is
the earliest true disease time for the patient. We sample $N=1000$ patients
with $M=4$ patient observations times with noise $S=0.25$ and upper time bound
$T^{+}=10$.

Figure | Description | Subtype Generating Functions
---|---|---
6 | Quadratic curve and flat line, separable | $f_{1}(t)=0.25t^{2}-2.2t+5,$
|  | $f_{2}(t)=2$
7 | Quadratic curve and flat line, overlapping | $f_{1}(t)=0.25t^{2}-2.2t+5,$
|  | $f_{2}(t)=-2$
8 | Quadratic curve and sloped line, separable | $f_{1}(t)=0.25t^{2}-2.2t+5,$
|  | $f_{2}(t)=0.4t$
9 | Quadratic curve and sloped line, overlapping | $f_{1}(t)=0.25t^{2}-2.2t+5,$
|  | $f_{2}(t)=0.4t-5$
10 | Quadratic curves in opposite directions, separable | $f_{1}(t)=0.25t^{2}-2.2t+3,$
|  | $f_{2}(t)=-0.25t^{2}+2.2-5$
11 | Quadratic curves in opposite directions, overlapping | $f_{1}(t)=0.25t^{2}-2.2t+7,$
|  | $f_{2}(t)=-0.25t^{2}+2.2-5$
Table 5: Quadratic dataset subtype generating functions and corresponding
figure numbers

We construct our quadratic experiments such that we examine different model
classes (i.e. flat, linear, quadratic) as well as examine subtypes that are
overlapping or separable.

We include baseline results for the quadratic datasets. Note that SuStaIn
[Young et al., 2018] assumes monotonically increasing functions and is
therefore omitted. We denote degenerate solutions with dashes.

###  D.2 Optimal hyperparameters for quadratic datasets

For the synthetic quadratic dataset corresponding to Figure 4, we found the
optimal hyperparameters for SubLign as no regularization, 5 hidden dimensions
for the multi-layer perceptron, 200 latent dimensions, 200 units for the
recurrent neural network, and learning rate of 0.001.

For the synthetic quadratic dataset corresponding to Figure 5, we found the
optimal hyperparameters for SubLign as no regularization, 5 hidden dimensions
for the multi-layer perceptron, 200 latent dimensions, 200 units for the
recurrent neural network, and learning rate of 0.001.

For the synthetic quadratic dataset corresponding to Figure 6, we found the
optimal hyperparameters for SubLign as no regularization, 5 hidden dimensions
for the multi-layer perceptron, 200 latent dimensions, 200 units for the
recurrent neural network, and learning rate of 0.001.

For the synthetic quadratic dataset corresponding to Figure 7, we found the
optimal hyperparameters for SubLign as no regularization, 5 hidden dimensions
for the multi-layer perceptron, 200 latent dimensions, 200 units for the
recurrent neural network, and learning rate of 0.001.

For the synthetic quadratic dataset corresponding to Figure 8, we found the
optimal hyperparameters for SubLign as no regularization, 5 hidden dimensions
for the multi-layer perceptron, 200 latent dimensions, 200 units for the
recurrent neural network, and learning rate of 0.001.

For the synthetic quadratic dataset corresponding to Figure 8, the optimal
hyperparameters are latent space of dimension 10, 20 hidden units in the RNN,
50 hidden units in the multi-layer perceptron, learning rate of 0.01, and no
regularization.

For the synthetic quadratic dataset corresponding to Figure 9, the optimal
hyperparameters are latent space of dimension 5, 100 hidden units in the RNN,
50 hidden units in the multi-layer perceptron, learning rate of 0.01, and no
regularization.

###  D.3 Results

In Figures 4 to 9, we present the quantitative results of 6 different
quadratic cases as well as a plot of example data and the data generating
subtypes. In each, we see that SubLign or SubNoLign outperforms the baselines.

When the subtypes are separable (i.e. Fig 5, 7, and 9), SubLign handily
recovers the subtypes. When the subtypes are not separable (i.e. Fig 4, 6, and
8), SubLign still outperforms baselines.

We note that alignment metrics are especially challenging to recover when one
subtype is a flat or sloped line as in with Figure 4 to 9. Because the
alignment metric is entirely unidentifiable, the swaps and Pearson metrics
suffer. Note that for the swaps metric, 0.5 corresponds to random guessing, so
the lack of identifiability of one of the subtypes would cause a swaps metric
of 0.25. When the second subtype has a changing slope, as in Figure 8, the
alignment metrics are more recoverable.

When the model is degenerate and does not return the alignment values, we
denote this with an empty cell.

Model | ARI $\uparrow$ | Swaps $\downarrow$ | Pearson $\uparrow$
---|---|---|---
SubLign | 0.277 $\pm$ 0.040 | 0.277 $\pm$ 0.014 | 0.516 $\pm$ 0.007
SubNoLign | 0.103 $\pm$ 0.002 | – | –
KMeans+Loss | 0.213 $\pm$ 0.009 | 0.498 $\pm$ 0.022 | 0.016 $\pm$ 0.051
SuStaIn [Young et al., 2018] | 0.151 | 0.203 | 0.000
Bayesian [Huopaniemi et al., 2014] | 0.000 $\pm$ 0.000 | 0.501 $\pm$ 0.017 | 0.018 $\pm$ 0.125
PAGA [Wolf et al., 2019] | 0.027 $\pm$ 0.001 | – | –

Figure 4: Synthetic results over 5 trials. Top: Data generating functions for
two subtypes (thick lines) and example aligned patients (dots and thin lines).
Bottom: SubLign outperforms baselines while KMeans+Loss recovers subtypes (ARI
metric) better than SubNoLign, but alignment metrics are difficult to recover
because of the horizontal subtype

MODEL | ARI $\uparrow$ | SWAPS $\downarrow$ | PEARSON $\uparrow$
---|---|---|---
SubLign | 0.980 $\pm$ 0.000 | 0.253 $\pm$ 0.001 | 0.527 $\pm$ 0.011
SubNoLign | 0.980 $\pm$ 0.000 | – | –
KMeans+Loss | 0.883 $\pm$ 0.000 | 0.471 $\pm$ 0.011 | 0.064 $\pm$ 0.067
SuStaIn [Young et al., 2018] | 0.228 $\pm$ 0.039 | 0.182 $\pm$ 0.010 | 0.000 $\pm$ 0.000
Bayesian [Huopaniemi et al., 2014] | 0.198 $\pm$ 0.189 | 0.446 $\pm$ 0.052 | 0.157 $\pm$ 0.286
PAGA [Wolf et al., 2019] | 0.227 $\pm$ 0.035 | – | –

Figure 5: Synthetic results over 5 trials. Top: Data generating functions for
two subtypes (thick lines) and example aligned patients (dots and thin lines).
Bottom: SubLign and SubNoLign have near-perfect clustering accuracy (ARI)
while alignment metrics (swaps, Pearson) are difficult recover because of the
horizontal subtype.

Model | ARI $\uparrow$ | Swaps $\downarrow$ | Pearson $\uparrow$
---|---|---|---
SubLign | 0.122 $\pm$ 0.006 | 0.272 $\pm$ 0.001 | 0.621 $\pm$ 0.020
SubNoLign | 0.145 $\pm$ 0.006 | – | –
KMeans+Loss | 0.031 $\pm$ 0.030 | 0.302 $\pm$ 0.026 | 0.498 $\pm$ 0.010
SuStaIn [Young et al., 2018] | 0.138 $\pm$ 0.019 | 0.119 $\pm$ 0.006 | 0.000 $\pm$ 0.000
Bayesian [Huopaniemi et al., 2014] | 0.001 $\pm$ 0.002 | – | –
PAGA [Wolf et al., 2019] | 0.009 $\pm$ 0.000 | – | –

Figure 6: Synthetic results over 5 trials. Top: Data generating functions for
two subtypes (thick lines) and example aligned patients (dots and thin lines).
Bottom: SubLign outperforms baselines in clustering and alignment metrics
although the task is challenging.

Model | ARI $\uparrow$ | Swaps $\downarrow$ | Pearson $\uparrow$
---|---|---|---
SubLign | 0.968 $\pm$ 0.016 | 0.253 $\pm$ 0.015 | 0.604 $\pm$ 0.041
SubNoLign | 0.968 $\pm$ 0.016 | – | –
KMeans+Loss | 0.964 $\pm$ 0.008 | 0.486 $\pm$ 0.032 | 0.044 $\pm$ 0.101
SuStaIn [Young et al., 2018] | 0.220 $\pm$ 0.011 | 0.196 $\pm$ 0.008 | 0.000 $\pm$ 0.000
Bayesian [Huopaniemi et al., 2014] | 0.221 $\pm$ 0.214 | 0.448 $\pm$ 0.070 | 0.165 $\pm$ 0.193
PAGA [Wolf et al., 2019] | 0.205 $\pm$ 0.012 | – | –

Figure 7: Synthetic results over 5 trials. Top: Data generating functions for
two subtypes (thick lines) and example aligned patients (dots and thin lines).
Bottom: SubLign, SubNoLign, and KMeans+Loss perform well on clustering.

Model | ARI $\uparrow$ | Swaps $\downarrow$ | Pearson $\uparrow$
---|---|---|---
SubLign | 0.729 $\pm$ 0.049 | 0.153 $\pm$ 0.006 | 0.843 $\pm$ 0.016
SubNoLign | 0.721 $\pm$ 0.035 | – | –
KMeans+Loss | 0.540 $\pm$ 0.034 | 0.490 $\pm$ 0.020 | 0.043 $\pm$ 0.057
SuStaIn [Young et al., 2018] | 0.198 $\pm$ 0.024 | 0.200 $\pm$ 0.008 | 0.000 $\pm$ 0.000
Bayesian [Huopaniemi et al., 2014] | 0.003 $\pm$ 0.007 | – | –
PAGA [Wolf et al., 2019] | 0.059 $\pm$ 0.006 | – | –

Figure 8: Synthetic results over 5 trials. Top: Data generating functions for
two subtypes (thick lines) and example aligned patients (dots and thin lines).
Bottom: SubLign learns subtypes and recovers alignment better than baselines.

Model | ARI $\uparrow$ | Swaps $\downarrow$ | Pearson $\uparrow$
---|---|---|---
SubLign | 1.000 $\pm$ 0.000 | 0.134 $\pm$ 0.012 | 0.907 $\pm$ 0.009
SubNoLign | 0.984 $\pm$ 0.008 | – | –
KMeans+Loss | 1.000 $\pm$ 0.000 | 0.508 $\pm$ 0.034 | 0.014 $\pm$ 0.095
SuStaIn [Young et al., 2018] | 0.148 $\pm$ 0.032 | 0.247 $\pm$ 0.016 | 0.000 $\pm$ 0.000
Bayesian [Huopaniemi et al., 2014] | 0.261 $\pm$ 0.349 | 0.409 $\pm$ 0.052 | 0.303 $\pm$ 0.119
PAGA [Wolf et al., 2019] | 0.233 $\pm$ 0.017 | – | –

Figure 9: Synthetic results over 5 trials. Top: Data generating functions for
two subtypes (thick lines) and example aligned patients (dots and thin lines).
Bottom: SubLign learns subtypes and alignment values.

##  Appendix E Data Privacy and Ethics

The heart failure dataset was collected from a large health system in the
United States and shared with us under a research data use agreement. The
hospital obtained the relevant consent from individuals to use their data, and
the dataset is covered by the hospital’s institutional review board (IRB).
Since the dataset was de-identified and provided with a data use agreement,
our institution’s IRB ruled it as exempt. The dataset use agreement was
approved by the legal teams from both our academic institution and the
hospital.

The Parkinson’s Progression Markers Initiative (PPMI) provides open and full
access to the study data, which is intended for researchers to study the
disease. The PPMI dataset only includes patients who consent to including
their data in the study, and the patients are de-identified in the dataset.

##  Appendix F SubLign Clinical Insights

In Figure 11, we show the statistically significant baseline features from our
SubLign results on heart failure patients. In Figure 10, we show the
statistically significant baseline features for our SubLign results on both
Parkinson’s disease patients and healthy controls as well as only Parkinson’s
disease patients.

Feature | A (321) | B (298)
---|---|---
Biological Dad With PD | 0.02 | 0.06
Sibling With PD | 0.01 | 0.05
UPSIT Part 1 | 7.55 | 5.49
UPSIT Part 2 | 7.64 | 5.69
UPSIT Part 3 | 6.98 | 5.23
UPSIT Part 4 | 7.53 | 5.62
UPSIT Total | 29.73 | 22.05
(a) width=0.45

Feature | A (156) | B (112) | B (155)
---|---|---|---
Male | $0.51$ | $0.68$ | $0.79$
White | $0.98$ | $0.95$ | $0.92$
UPSIT Part 1 | $6.01$ | $5.07$ | $5.43$
UPSIT Part 2 | $6.65$ | $5.30$ | $5.52$
UPSIT Part 3 | $5.92$ | $4.79$ | $4.90$
UPSIT Part 4 | $6.28$ | $5.20$ | $5.42$
UPSIT Total | $24.87$ | $20.36$ | $21.26$
(b) width=0.45

Figure 10: Subtypes found by SubLign for (a) 619 Parkinson’s disease patients and healthy controls and (b) 423 Parkinson’s disease patients only. Only statistically significant means between subtypes according to an ANOVA test with $p<0.05$ with a Benjamini-Hochberg correction are listed. Feature | A (321) | B (298)
---|---|---
Biological Dad With PD | 0.028 | 0.068
Sibling With PD | 0.010 | 0.058
UPSIT Part 1 | 7.558 | 5.493
UPSIT Part 2 | 7.648 | 5.695
UPSIT Part 3 | 6.988 | 5.238
UPSIT Part 4 | 7.539 | 5.624
UPSIT Total | 29.73 | 22.05
Table 6: Subtypes found by SubLign from Parkinson’s disease patients and healthy controls using sparsely collected biomarkers. Only statistically significant means between subtypes according to an ANOVA test with $p<0.05$ are listed. Feature | A (674) | B (444) | C (416)
---|---|---|---
Age | 75.98 | 74.73 | 69.43
Female | 0.71 | 0.23 | 0.43
Anemia | 0.23 | 0.16 | 0.14
Atherosclerosis | 0.28 | 0.34 | 0.40
Atrial Fibrillation | 0.44 | 0.55 | 0.43
Chronic Kidney Disease | 0.27 | 0.34 | 0.34
Diastolic Heart Failure | 0.50 | 0.36 | 0.06
Obese | 0.56 | 0.65 | 0.46
Old Myocardial Infarction | 0.12 | 0.14 | 0.24
Pulmonary Heart Disease | 0.29 | 0.22 | 0.19
Systolic HF | 0.09 | 0.27 | 0.53
Figure 11: Subtypes found by SubLign from heart failure patients using
echocardiogram biomarkers. Only statistically significant means between
subtypes according to an ANOVA test with $p<0.05$ with a Benjamini-Hochberg
correction are listed.

In comparison, we show the KMeans+Loss subtypes on the heart failure dataset
(Table 7).

Feature | A (240) | B (802) | C (492)
---|---|---|---
Age | 71.567 | 74.565 | 73.793
Hyperlipidemia | 0.529 | 0.448 | 0.541
Chronic Kidney Disease | 0.346 | 0.273 | 0.370
Esophageal Reflux | 0.375 | 0.259 | 0.289
Pulmonary Heart Disease | 0.367 | 0.204 | 0.256
Kidney Disease | 0.254 | 0.200 | 0.278
Atherosclerosis | 0.196 | 0.131 | 0.201
Anemia | 0.217 | 0.163 | 0.213
Obese | 0.688 | 0.500 | 0.608
Table 7: Heart Failure KMeans+Loss subtypes (patient counts in parentheses),
described by mean baseline features. Only statistically significant features
are listed and do not include systolic and diastolic HF, two known phenotypes
of HF.