A Similarity-Driven Sequential Model for
Predicting Medical Treatment Initiation via
Graph LaplacianTong Wu1,2 , Yue Wang1 , Yunlong Wang1 , Cao Xiao1 , Emily Zhao1 , Lucas
Glass1 , and Jimeng Sun3
1

2

IQVIA Inc., Plymouth Meeting, PA, USA
Biomedical Engineering, University of Minnesota, Minneapolis, MN, USA
3
Georgia Institute of Technology, Atlanta, GA, USA

Abstract. A diverse range of computational models have been proposed
to extract temporal patterns from clinical time series data for each patient and among patient group for predictive healthcare, including multivariate regression, neural network, and decision tree, etc. These models
typically learn an analytical mapping from the numerical or categoricalpatients’ medical data to the targeted labels. However, the common relations among patients (e.g., share the same doctor, or receive the same
medical treatments) were rarely considered. In this paper, we represent
the relation between patients and clinicians by bipartite graphs to explicitly address relations such as from which clinician a patient is diagnosed
or offered a medical service. We then solve for the top eigenvectors of
the Laplacian matrix of the bipartite graph, and include the eigenvectors as latent representations of the similarity between patient-clinician
pairs into a time-sensitive prediction model. We conducted experiments
using real-world data to predict the initiation of first-line treatment for
Chronic Lymphocytic Leukemia (CLL) patients. Results show that therelational similarity uncovered from the patient-clinician bipartite graphcan enhance prediction accuracies over multiple baselines, for example,
a 5% increment over the long short-term memory sequential model interms of area under precision-recall curve.
Keywords: Healthcare data mining · Time Series Analysis · Graph
Laplacian.

1

Introduction

Recent years have witnessed an explosion in the growth of the amount of digital
information in electronic health records, which provides great opportunities forapplications such as healthcare analytics and clinical informatics [25, 24]. Amongthe electronic health data, clinical time series is one major type. More specifically,
the medical records of a patient visiting one or more clinical sites, including
diagnoses, laboratory tests, prescriptions, etc., are time-stamped and can bearranged in sequential orders.

2T. Wu et al.

Over the years, many computational models have been proposed for disease detection [10, 11], disease progression [1], patient subtyping [2, 4], and many
more, based on clinical time series data. A core challenge regarding studies oftime series data is to determine the similarity of one time series with respect
to another, which is central to tasks such as forecasting, mining, and clustering.
Over the years several similarity measures for time series data have been proposed, such as cross-correlation, dynamic time warping [20], and edit distance[12, 5]. A good candidate of similarity measure is expected to be able to handle
high-dimensional data, while being fast, scalable, and efficient [21].
A disease progression model is designed to predict the development of potential treatments for many slowly progressing diseases, e.g. Alzheimer’s disease,by detecting more granular stages as compared to those defined in clinical diagnosis [22]. A primary goal of oncology pharmaceutical market research is todevelop a quantitative prediction model that leverages machine learning methods and patients medical records to perform automated and accurate predictionof disease progression, thus can help raise awareness of treatment options to therelevant clinicians, and is also beneficial to early disease intervention. Longitudinal clinical time series along with related patient or physician data is important
in informing disease progression patterns, and motivating many deep learning
based disease progression models including recurrent neural network (RNN),
attention model, and graph embedding [7, 9, 1, 4, 23].
Despite these initial success, the relational structure between patients andclinicians has long been overlooked and is of great potentials to enhance prediction accuracy. The assumption is that patients with the same disease who visit
the same clinicians tend to receive similar treatments, which can be arguably attributed to that clinicians follow a set of common medical knowledge, and more
likely make similar decisions for patients with the same disease.
We propose to model the patient-clinician relational structure as a bipartite
graph, in which the two disjoint sets of vertices represent the patients and theclinicians, respectively; the set of edges records the number of visits made byeach patient with each clinician. We note that many patients visit different
clinicians for diagnosis and follow-up treatments. Thus we create two graphs,
one for patients and their diagnosis clinicians, and another for patients and theirfollow-up clinicians. An example of the graph using patients’ clinical data fromIQVIA’s prescription and claim database is given in Figure 1. Here we show six
patients diagnosed with Chronic Lymphocytic Leukemia (CLL), and their visitedclinicians for both diagnosis and follow-up treatments. CLL is a type of cancer
in which the bone marrow makes too many lymphocytes. For CLL patients,
there is normally an observational period before starting first-line treatment.
We restrict the follow-up clinicians to oncologists and hematologists for better
relevance. We then apply spectral graph analysis to solve for the top eigenvectorsof the graph Laplacian, and take the eigenvectors as latent representations of
the similarity between patient-clinician pairs. The extracted features capture thelatent proximity of patients who visit the same clinicians, and can enhance theprediction accuracy of a variety of disease progression models.

Title Suppressed Due to Excessive Length3

CLL Patient
Diagnosis
ClinicianFollow-up
ClinicianFig. 1. Illustration of the bipartite graph between CLL patients and diagnosis/followup clinicians. Follow-up clinicians are restricted to oncologists and hematologists. Diagnosis and follow-up clinicians can overlap.

As a key contribution, to our best knowledge, this work is the first diseaseprogression/detection model that leverages graph theory to exploit the relational
similarity of clinician-visiting from patients’ clinical time series data. Using thepatients’ clinical medical data from IQVIA’s database, we show that the proposed similarity features can improve the performance of a wide choice of machine learning models, from XGBoost [6] to more recent deep learning basedmodels, e.g. convolutional neural network (CNN) and long short-term memory(LSTM).

2

Related Works

Medical treatment prediction is a core research task of disease progression modeling. Recently, many deep learning models have made rapid advancements on
this topic. In [9], a two-level attention model was designed to detect influential
past visits and significant clinical variables for better prediction accuracy andinterpretability. In [8], a graph-based attention model was proposed to extract
hierarchical information from medical oncologies and improve RNN-based raredisease prediction. In [15], a bi-directional RNN was designed to remember information of both the past and future visits based on three attention mechanism
to measure the relationship of different visits for prediction. In [4], a RNN architecture through dynamically matching temporal patterns was proposed to learn
the similarity between two longitudinal patient record sequences for personalized prediction of Parkinson’s Disease. Other similar approaches have also been
proposed [17, 16].

4

T. Wu et al.
Table 1. List of notations defined in this paper.
Notation

Description

D
S
I
U, V
M = |U|, N = |V|N
E = {{wij }M
i=1 }j=1

Clinical medical records
Patients’ sequences of visits
Patients’ demographics
Sets of clinicians and patients.
Sizes of U and V.
Set of edges connecting U and V; wij is thecount of patient j visiting clinician i.
Bipartite graph for patients and their visited clinicians
Adjacency matrix of G
Degree matrix of A
Laplacian matrix of G
Relational similarity features

G = {U, V, E}
A
DL{Xi }K
i=1

3

Methods

3.1

Problem Formulation

We let D denote patients’ clinical medical records. For a patient p, D(p) =(p) (p)(p){S(p) , I (p) }, where S(p) = [s1 , s2 , . . . , sT ] is the sequence of visits and I (p)(p)
includes demographics and related medical features. Each visit st consists of
information such as diagnoses, procedures, prescriptions, visited clinicians, etc.
The goal of prediction is to learn f : D(p) 7→ y (p) , where y (p) is the label indicating if the patient will start treatment in the next time window.
We model the patient-clinician relation as a bipartite graph G = {U, V, E}
over D, where U, V are the sets of clinicians and patients, E is the set of edges
connecting U and V, with each weight wij denoting the count of patient j visiting
clinician i. Table 1 lists the notations we used in the paper.
3.2

Extracting Relational Similarities from Patients’ Clinical Time
Series via Graph Laplacian

The algorithm for extracting relational similarity is motivated by spectral clustering [19], in which data points are considered as nodes of a similarity graphand mapped to a low-dimensional space where they can be segregated to form
clusters.
Firstly, we construct K bipartite graphs {Gi }K
i=1 from D, where each Gi =|Ui | |V|
{Ui , V, Ei } and Ei = {{wij }i=1
}j=1 . We let M = |Ui | and N = |V| for simplicity.
V represents patients thus remains the same across all graphs.
Secondly, for each graph Gi we construct the adjacency matrix A. Since in
Gi we ignore interaction within patients (or clinicians), A is block sparse, andTitle Suppressed Due to Excessive Length5

Algorithm 1 Graph-based Similarity Feature Extraction
Input: Bipartite graphs {Gi }K
i=1
Output: {Xi }K
i=1
for i = 1 to K do
Gi = {Ui , V, Ei }, M = |Ui |, N = |V|
0M,M BA=(where B[i, j] = wij )BT 0N,NMP+N+N
D = diag({
A[i, j]}M
i=1 )
j=1

L = D−1/2 AD−1/2
Xi = [e1 , e2 , . . . , ek ] = eigendecomp(L)for α = 1 to M + N doP
Xi [α, :] = Xi [α, :]/( β Xi [α, β]2 )1/2
end forend for
0M,M Btakes the form as, where B ∈ RM ×N is the matrix representation
BT 0N,Nof Ei .
Thirdly, we compute the Laplacian matrix of Gi as L = D−1/2 AD−1/2 , where
D is a diagonal matrix whose (i, i)-th element is the sum of the i-th row of A.
Lastly, we compute the top k eigenvectors of L, i.e. Xi = [e1 , e2 , . . . , ek ] ∈R(M +N )×k , using algorithms such as Implicitly Restarted Arnoldi Method [14]that is suitable to process large sparse matrices. Xi is normalized such that eachrow of Xi has an L2 -norm of 1.
For a patient p, we extract the corresponding row from {Xi }K
i=1 as the relational similarity features. The complete algorithm is given in Algorithm 1.

4

Experiments and Results

4.1

CohortWe extract data from IQVIA longitudinal prescription (Rx) and medical claims(Dx) database, including hundreds of millions of patient clinical records. In thisstudy, we selected all the patients diagnosed with CLL from 01-2017 to 12-2018from the IQVIA database and kept only the patients with complete Rx/Dx
information. We split the time period from 07-2017 to 12-2018 into 3 equalintervals. Within each interval, we defined the positive cohort as patients who
were diagnosed with CLL before the interval and started treatment within theinterval, and the negative cohort as patients who were diagnosed with CLL before
the interval, but did not start treatment during the interval. The final positive
and negative cohorts have 11,259 and 109,563 patient profiles, each of whichcontains a medical sequence and a feature vector including demographics andrelational similarity features, respectively.
In time series forecasting, it is common practice to reserve the last part ofeach time series for testing, and use the rest of the series for training [3]. This is6

T. Wu et al.

to avoid using information from the future to predict past events. As a result, wecannot apply conventional train/test split of data such as k-fold cross-validationthat would shuffle data randomly. Therefore, the 07-2018 to 12-2018 interval ishold out for model testing, and the rest intervals are used for model training.
4.2

Feature preparation

BOW features: For each interval, we pulled one-year patient clinical records(diagnosis, procedure, etc.) before the interval. We divided the one-year lookback into four quarters, and used Bag-of-Words (BOW) [18] to extract features(counts of occurrences of services) from the records within each quarter. Medical services with less than 1% occurrences in the entire dataset were ignored.
Eventually, we obtained 329 BOW features.
Demographics and medical features: Demographic features include age
and gender. Medical features are 11 clinical services from patients’ medical history, which are chosen by clinical experts and have been shown relevance totreatment initiation. Both demographics and medical features are converted tocategorical variables.
Vectorized clinical time series: The medical service codes in the patients’clinical time series were converted to 300-dimensional dense vectors through a
pre-trained embedding.
Relational similarity: We constructed two bipartite graphs, one for patients and their diagnosis clinicians, and another for patients and their follow-up
clinicians. We took the top 5 eigenvectors from each of the two Laplacian matrices derived from the two graphs using Algorithm 1, resulting in a 10-dimensional
vector as the relational similarity feature for one patient.
4.3

Evaluation

Because of the imbalance of positive and negative samples in the cohort, we use
precision-recall area-under-curve (PR-AUC) and precision@k as the metrics toevaluate model performance.
All models are built using the training data set and evaluated using thetesting data set. We use Adam optimizer [13] with a default learning rate of
1e-4. The number of training epochs for each model is 50 and an early stopping
criterion is invoked if the performance does not improve in 10 consecutive epochs.
All models are implemented in Keras with Tensorflow backend and tested ona system equipped with 128GB RAM, 16 Intel(R) Core Xeon(R) E5-2683 v4
2.10GHz CPUs, and Nvidia Tesla P100-PCIE-16GB.
4.4

Models and Performance

To demonstrate the robustness of the proposed method, we have applied it on
three baseline models for comparison: XGB, CNN, and LSTM. We evaluate theperformance of these models with and without the relational similarity feature.

Title Suppressed Due to Excessive Length

7

Precision-Recall Curve
1.0

Area Under Curve
XGB: 0.2235
XGB-Proposed: 0.2427
CNN: 0.2322
CNN-Proposed: 0.2437
LSTM: 0.2484
LSTM-Proposed: 0.2652

Precision

0.8

0.6

0.4

0.2

0.0

0.0

0.2

0.4

Recall

0.6

0.8

1.0

Fig. 2. Comparison of baseline models with and without relational similarity feature
measured in PR-AUC.

For fairly comparison, we also introduce the diagnosis clinicians and follow-up
clinicians IDs to the baseline models, to guarantee that the proposed method
and baseline model leverage the same amount of information.
XGB: A tree boosting regression model implemented in XGBoost with 500
estimators. Input features include BOW features, demographics, medical features, and relational similarity features.
CNN: A 4-layer CNN with a structure of 1×1(128)−2×2(128)−3×3(128)−
5 × 5(128) (kernel size/number of kernels), followed by two dense layers. Clinicaltime series are fed into the first CNN layer. Features of relational similarity,
demographics, and medical are concatenated with flattened CNN features at thefirst dense layer.
LSTM: A bi-directional LSTM handles the clinical time series. The hidden dimension of the LSTM is 256. The maximum hidden states at each time
step of the top-level LSTM are concatenated with relational similarity features,
demographics, and medical features, and processed by two dense layers.
Performance of the models with and without the relational similarity feature
is evaluated in PR-AUC and precision@k, and shown in Figure 2 and 3. Among
all the models, the feature number of proposed and baseline are on comparable
level, e.g., XGB-baseline has 1,331, and XGB-proposed has 1,339 features.
A sampling of the precision@k curve is provided in Table 2. The results further confirm that by capturing the relational similarity feature, we can improve
the prediction accuracy of all baseline models, suggesting the universal effectiveness of relational similarity feature for a wide range of model structures.
When K goes large, the proposed method consistently outperforms the baseline

8

T. Wu et al.

Precision@k

1.0

XGB
XGB-Proposed
CNN
CNN-Proposed
LSTM
LSTM-Proposed

True treated patients

0.8

0.6

0.4

0.2

0.0

1000

2000

3000

Potential patients predicted

4000

5000

Fig. 3. Comparison of baseline models with and without relational similarity feature
measured in precision@k.

model, although the incremental gradually become marginal, e.g., in LSTM the
improvement is 2.6%, 2.2%, 0.6%, when K equals 5,000, 10,000, 50,000, respectively.

5

Conclusions

In this paper, we proposed a graph-based algorithm to extract latent relational
similarity from patients’ clinical time series. Experimental results using realworld data show that the proposed feature can improve the prediction accuracyof a wide range of model structures. We envision that the relational similarity canalso enhance model performance on other tasks, such as rare disease detection
or patient subtyping. In its current form, the algorithm is operating on a twodimensional feature space, i.e. patients and clinicians. In the future, we will add
medical services as another dimension into the feature space, and thus enable
the application of more advanced signal processing techniques such as tensor
decomposition to uncover more useful information from the multidimensionalfeature inputs.

References
1. Bai, T., Zhang, S., Egleston, B.L., Vucetic, S.: Interpretable representation learningfor healthcare via capturing disease progression through time. In: Proceedings ofTitle Suppressed Due to Excessive Length

9

Table 2. The number of successfully predicted treatments from 07-2018 to 12-2018
among the K selected patients.
K

600

1,200

1,800

3,600

XGB
XGB-Proposed
Improvement
CNN
CNN-Proposed
Improvement
LSTM
LSTM-Proposed
Improvement

346
360
4.0%
325
332
2.2%
346
382
10.4%

534
559
4.7%
526
562
6.8%
557
598
7.4%

669
725
8.4%
685
724
5.7%
710
762
7.3%

975
1,057
8.4%
1,060
1,079
1.8%
1,059
1,110
4.8%

the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data
Mining. pp. 43–51. ACM (2018)
2. Baytas, I.M., Xiao, C., Zhang, X., Wang, F., Jain, A.K., Zhou, J.: Patient subtyping via time-aware lstm networks. In: Proceedings of the 23rd ACM SIGKDD
international conference on knowledge discovery and data mining. pp. 65–74. ACM(2017)
3. Bergmeir, C., Benı́tez, J.M.: On the use of cross-validation for time series predictor
evaluation. Information Sciences 191, 192–213 (2012)
4. Che, C., Xiao, C., Liang, J., Jin, B., Zho, J., Wang, F.: An rnn architecture withdynamic temporal matching for personalized predictions of parkinson’s disease.
In: Proceedings of the 2017 SIAM International Conference on Data Mining. pp.
198–206. SIAM (2017)
5. Chen, L., Özsu, M.T., Oria, V.: Robust and fast similarity search for moving object
trajectories. In: Proceedings of the 2005 ACM SIGMOD international conference
on Management of data. pp. 491–502. ACM (2005)
6. Chen, T., Guestrin, C.: Xgboost: A scalable tree boosting system. In: Proceedingsof the 22nd acm sigkdd international conference on knowledge discovery and data
mining. pp. 785–794. ACM (2016)
7. Choi, E., Bahadori, M.T., Schuetz, A., Stewart, W.F., Sun, J.: Doctor ai: Predicting
clinical events via recurrent neural networks. In: Machine Learning for Healthcare
Conference. pp. 301–318 (2016)
8. Choi, E., Bahadori, M.T., Song, L., Stewart, W.F., Sun, J.: Gram: graph-based
attention model for healthcare representation learning. In: Proceedings of the 23rd
ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. pp. 787–795. ACM (2017)
9. Choi, E., Bahadori, M.T., Sun, J., Kulas, J., Schuetz, A., Stewart, W.: Retain: An
interpretable predictive model for healthcare using reverse time attention mechanism. In: Advances in Neural Information Processing Systems. pp. 3504–3512(2016)
10. Choi, E., Schuetz, A., Stewart, W.F., Sun, J.: Using recurrent neural networkmodels for early detection of heart failure onset. Journal of the American Medical
Informatics Association 24(2), 361–370 (2016)
11. Choi, E., Xiao, C., Stewart, W., Sun, J.: MiME: Multilevel medical embeddingof electronic health records for predictive healthcare. In: Bengio, S., Wallach, H.,
Larochelle, H., Grauman, K., Cesa-Bianchi, N., Garnett, R. (eds.) Advances in10

T. Wu et al.

Neural Information Processing Systems 31, pp. 4547–4557. Curran Associates, Inc.
(2018)
12. Fischer, A., Suen, C.Y., Frinken, V., Riesen, K., Bunke, H.: A fast matching algorithm for graph-based handwriting recognition. In: International Workshop on
Graph-Based Representations in Pattern Recognition. pp. 194–203. Springer (2013)
13. Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. arXiv preprint
arXiv:1412.6980 (2014)
14. Lehoucq, R.B., Sorensen, D.C.: Deflation techniques for an implicitly restarted
arnoldi iteration. SIAM Journal on Matrix Analysis and Applications 17(4), 789–
821 (1996)
15. Ma, F., Chitta, R., Zhou, J., You, Q., Sun, T., Gao, J.: Dipole: Diagnosis prediction in healthcare via attention-based bidirectional recurrent neural networks.
In: Proceedings of the 23rd ACM SIGKDD international conference on knowledge
discovery and data mining. pp. 1903–1911. ACM (2017)
16. Ma, F., Wang, Y., Xiao, H., Yuan, Y., Chitta, R., Zhou, J., Gao, J.: A generalframework for diagnosis prediction via incorporating medical code descriptions. In:
2018 IEEE International Conference on Bioinformatics and Biomedicine (BIBM).
pp. 1070–1075. IEEE (2018)17. Ma, F., You, Q., Xiao, H., Chitta, R., Zhou, J., Gao, J.: Kame: Knowledge-based
attention model for diagnosis prediction in healthcare. In: Proceedings of the 27th
ACM International Conference on Information and Knowledge Management. pp.
743–752. ACM (2018)18. Mikolov, T., Chen, K., Corrado, G., Dean, J.: Efficient estimation of word representations in vector space. arXiv preprint arXiv:1301.3781 (2013)
19. Ng, A.Y., Jordan, M.I., Weiss, Y.: On spectral clustering: Analysis and an algorithm. In: Advances in neural information processing systems. pp. 849–856 (2002)
20. Sakoe, H., Chiba, S., Waibel, A., Lee, K.: Dynamic programming algorithm optimization for spoken word recognition. Readings in speech recognition 159, 224(1990)
21. Serra, J., Arcos, J.L.: An empirical evaluation of similarity measures for time seriesclassification. Knowledge-Based Systems 67, 305–314 (2014)22. Sukkar, R., Katz, E., Zhang, Y., Raunig, D., Wyman, B.T.: Disease progression
modeling using hidden markov models. In: 2012 Annual International Conferenceof the IEEE Engineering in Medicine and Biology Society. pp. 2845–2848. IEEE
(2012)
23. Suresh, H., Hunt, N., Johnson, A., Celi, L.A., Szolovits, P., Ghassemi, M.: Clinicalintervention prediction and understanding with deep neural networks. In: Machine
Learning for Healthcare Conference. pp. 322–337 (2017)
24. Topol, E.J.: High-performance medicine: the convergence of human and artificial
intelligence. Nature Medicine (2019)
25. Xiao, C., Choi, E., Sun, J.: Opportunities and challenges in developing deep learning models using electronic health records data: a systematic review. Journal of
the American Medical Informatics Association 25(10), 1419–1428 (2018)