Journal of Biomedical Informatics 135 (2022) 104218 Contents lists available at ScienceDirect Journal of Biomedical Informatics journal homepage: www.elsevier.com/locate/yjbin Original Research Longitudinal deep learning clustering of Type 2 Diabetes Mellitus trajectories using routinely collected health records Enrico Manzini a,b,c ,∗,1 , Bogdan Vlacho d ,1 , Josep Franch-Nadal d,e,f , Joan Escudero g , Ana Génova g , Elisenda Reixach h , Erik Andrés h , Israel Pizarro i , José-Luis Portero i , Dídac Mauricio d,e,j ,∗∗,2 , Alexandre Perera-Lluna a,b,c ,2 a Universitat Politecnica de Catalunya, Barcelona, Spain b Networking Biomedical Research Centre in the subject area of Bioengineering, Biomaterials and Nanomedicine, Madrid, Spain c Institut de Recerca Sant Joan de Deu, Barcelona, Spain d DAP-Cat Group, Unitat de Suport a la Recerca, Fundaciò Institut Universitari per a la recerca a l’Atenciò Primària de Salut Jordi Gol i Gurina (IDIAPJGol), Barelona, Spain e Center for Biomedical Research on Diabetes and Associated Metabolic Diseases (CIBERDEM CB15/00071), Instituto de Salud Carlos III, 28029, Madrid, Spain f Primary Health Care Center Raval Sud, Institut Català de la Salut, Barcelona, Spain g Grupo Pulso, Spain h Fundació TIC Salut Social, Departament de Salut, Generalitat de Catalunya, Barcelona, Spain i Novo Nordisk, Spain j Department of Medicine, University of Vic - Central University of Catalonia, Vic, Spain ARTICLE INFO ABSTRACT Keywords: Type 2 diabetes mellitus (T2DM) is a highly heterogeneous chronic disease with different pathophysiological Type 2 diabetes and genetic characteristics affecting its progression, associated complications and response to therapies. The Deep learning advances in deep learning (DL) techniques and the availability of a large amount of healthcare data allow us Longitudinal cluster to investigate T2DM characteristics and evolution with a completely new approach, studying common disease AutoEncoder trajectories rather than cross sectional values. We used an Kernelized-AutoEncoder algorithm to map 5 years Diabetic complications Electronic health records of data of 11,028 subjects diagnosed with T2DM in a latent space that embedded similarities and differences between patients in terms of the evolution of the disease. Once we obtained the latent space, we used classical clustering algorithms to create longitudinal clusters representing different evolutions of the diabetic disease. Our unsupervised DL clustering algorithm suggested seven different longitudinal clusters. Different mean ages were observed among the clusters (ranging from 65.3±11.6 to 72.8±9.4). Subjects in clusters B (Hypercholesteraemic) and E (Hypertensive) had shorter diabetes duration (9.2±3.9 and 9.5±3.9 years respectively). Subjects in Cluster G (Metabolic) had the poorest glycaemic control (mean glycated hemoglobin 7.99±1.42%), while cluster E had the best one (mean glycated hemoglobin 7.04±1.11%). Obesity was observed mainly in clusters A (Neuropathic), C (Multiple Complications), F (Retinopathy) and G. A dashboard is available at dm2.b2slab.upc.edu to visualize the different trajectories corresponding to the 7 clusters. 1. Introduction 463 million [2]. Moreover, projections for the future indicate that the global impact of diabetes is likely to continue increasing considerably. Type 2 diabetes mellitus(T2DM) is a chronic and highly prevalent Nowadays, the evidence suggests that T2DM is a complex highly disease ranking as eighth in terms of the overall burden of diseases heterogeneous metabolic disease that encompasses different pathophys- measure developed by WHO [1]. According to the International Dia- iological and genetic pathways [3]. The presentation and progression betes Federation, estimates of the disease’s prevalence showed alarming of the disease can vary between subjects leading to poor glycemic and increases since 2000, tripling the number of affected subjects in 2019 to metabolic control [4]. Indeed in Europe alone, the European Medicines ∗ Corresponding author at: Universitat Politecnica de Catalunya, Barcelona, Spain. ∗∗ Corresponding author. E-mail addresses: enrico.manzini@upc.edu (E. Manzini), didacmauricio@gmail.com (D. Mauricio). 1 These authors have contributed equally to this work and share first authorship. 2 These authors have contributed equally to this work and share senior authorship. https://doi.org/10.1016/j.jbi.2022.104218 Received 20 June 2022; Received in revised form 8 September 2022; Accepted 3 October 2022 Available online 8 October 2022 1532-0464/© 2022 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by- nc-nd/4.0/). E. Manzini et al. Journal of Biomedical Informatics 135 (2022) 104218 Agency (EMA), for the period between 2005 and 2017, approved 40 lifestyles information, medications prescriptions/dispensations and re- new drugs to treat diabetes [5]. In Catalonia (Spain), according to ferrals to specialists, and laboratory test results, introduced by health previously published data, 43.9% of the T2DM subjects did not reach professionals during routine health surveillance and health care. an adequate glycemic control [6]. This issue might be since T2DM We selected only subjects diagnosed with type 2 diabetes (ICD-10 diagnosis is oversimplified via the assessment of blood levels of glucose diagnostic codes: E11 and E14 and their sub codes) [18]. In addition, only. subjects were not eligible if they were under 18 years old or they had The availability of routinely collected data from Electronic Health other types of diabetes such as type 1 diabetes, secondary or gestational Records (EHR) in recent years had a huge impact on clinical, pharma- diabetes. For each eligible patient, clinical and analytical variables coepidemiologic and health services research, but the accuracy, quality were available, such as glycated hemoglobin (HbA1c), body mass index and heterogeneity of these data remain challenging to handle [7]. In (BMI), diastolic and systolic blood pressure (DBP and SBP), lipid profile the case of the diabetic disease, for example, they have been recently and renal function. Furthermore, data on antidiabetic treatment (drugs used to point out the complexity of T2DM and its pathophysiological in the A10 group of ATC classification) [19] and diagnosis of the most phenotyping [3,8]. In 2018, a seminal study from Sweden, identified common comorbidities of diabetes were also available. ATC codes of five novel clusters of subjects [3]. In a large Danish study of newly the antidiabetic treatments grouped per pharmacological class and ICD- diagnosed T2DM, authors identified several distinct pathophysiological 10 codes of the comorbidities are reported in Supplementary Tables phenotypes according to beta cell function and insulin sensitivity [8]. S1 and S2 respectively. We included only patients that had more than However, these attempts to study different phenotype clusters in the 18 years at the moment of the T2Dm onset, with at least three measures diabetic population have been done with cross sectional studies, rather for a minimum of two variables in the period under analysis, and who than focusing on the clinical evolution of the disease. Even if some had at least one measure of HbA1C and BMI. The T2DM onset date attempts to study typical T2DM trajectories exists [9–11], these are was estimated for each subject as the minimum date between the date limited to few factors that influence the evolution of the diseases, of the first T2DM diagnosis, the first prescription of the antidiabetic e.g. a limited subset of common comorbidities [9], or the evolution drugs, or the first measure of HbA1c>7% . The latter value was adopted of few variables as the glycated hemoglobin (HbA1c) [10] or the as for an important part of the study subjects an initial HbA1c was body mass index (BMI) [11]. To the best of our knowledge, there are determined before 2010, when HbA1c was standardized and accepted no studies focused on the whole evolution of the disease and that as a diagnostic criterion for diabetes. Thus, this algorithm was stricter take into consideration several factors including comorbidities, medical but, at the same time, more precise in identifying the time of diagnosis treatments, HbA1c changes over time and other important variables for of T2DM in our database. the description of the disease. 2.2. Data preprocessing Deep learning (DL) has been proved to be very effective in the design of multivariate patient trajectories, especially for databases with Since routinely collected data extracted from EHR can contain many missings values, which is very common among data collected errors arising from several different factors, we used clinical criteria from EHRs, e.g. [12,13]. Even if the in recent years the number of to remove any extreme values that could potentially be measurement DL models to extract temporal information from EHRs has notably errors. Removed values are reported in Supplementary Tables S3. The increased, several challenges remain to be handled, with the major study subjects diagram is presented in the Supplementary Figures S1. ones being irregularity, sparsity and heterogeneity of the data and the To mitigate the problem of irregularity of the data we represented each opacity of the models due to their black box nature [14]. patient as a multivariate time series (MTS) 𝑃𝑖 ∈ R𝑛𝑡,𝑖 ×𝑛𝑣 where: The proper clusterization of T2DM subjects by their clinical char- acteristics is a major issue. Progress in this field should enable us to • 𝑛𝑣 is the number of variables considered for the analysis. We implement different strategies of precision medicine and include diag- used the following variables: age, age at diagnosis of the sub- nostic algorithms for defining diabetes subtypes in order to implement jects, sex, diabetes duration, five binary signals representing the most appropriate clinical decisions [15]. common diabetes secondary conditions (cardiovascular disease The objective of this study is to exploit the potential of unsupervised (CVD), chronic kidney disease (CKD), neuropathy, high blood deep learning algorithms to create different longitudinal phenotype pressure (HTN), and ophthalmological complications),nine phar- clusters of T2DM patients using routinely collected data extracted from macological classes of antidiabetic drugs and 10 clinical and ana- EHR. Hence, throughout this work, we make the following contri- lytical variables (HbA1c, BMI, SBP, DBP, High-density lipopro- butions: on one hand we adapt a novel DL model, the kernelized tein (HDL), Low-density lipoprotein (LDL), Total cholesterol, autoencoder (K-AE) [16] to work with EHR data, handling the problem triglycerides (TG), Creatinine and albumin to creatinine ratio of sparsity and irregularity of the data with a different kernel matrix (ACR)). and a different imputation process and trying to reduce the opacity of • 𝑛𝑡,𝑖 is the number of time steps. Since every patient had an the model showing the results and the AE latent space in an interactive almost unique combination of number and frequency of measured dashboard available at dm2.b2slab.upc.edu; on the other hand, we variables, we standardized medical examinations dates, creating use this model to calculate longitudinal clusters of the T2DM disease, a temporal grid from January 1st, 2013 to December 31st, 2017, studying different typical trajectories. with a fixed frequency of six months. Then, we re-sampled the date of each examination according to this grid. In the case of 2. Material and methods two or more examinations inside the same time interval of the temporal grid, we used their mean value. Hence, we linearly interpolated these points. For this reason the maximum number 2.1. Data source and study population of time steps is 10, but it can be lower, for example in the cases of patients who died before the end of the study. Data from this study were extracted from the Information System for the Development of Research in Primary Care (SIDIAP) database [17] The combination of measured variables and dates of the measure- from January 1st 2013 to December 31th 2017. This database contains ments changed from patient to patient, implying the possibility of data from Electronic Health Records collected from approximately missing values in the matrix Pi. We imputed such data using a two- 5.6 million patients registered from 287 Primary Care Centers (PCC) step algorithm: firstly, we imputed missing values for a variable 𝑣𝑛 in Catalonia. The SIDIAP database contains data on demographics, when a patient had at least two existing measurements. For this we patient visits with health professionals, diagnoses, clinical variables, defined for each variable a mixed model, using the other variables as 2 E. Manzini et al. Journal of Biomedical Informatics 135 (2022) 104218 fixed effects and the patient as random effect. For each variable the of hyperparameters, such as the dimensions of the different layers and best mixed model was chosen using the Bayesian Information Criterion the importance given to the kernel alignment, can affect the structure (BIC). In the second step, we estimated values of 𝑣𝑛 for patients having of the latent space, changing the distribution of the clusters. Optimal no, or just one, measures of that variable by using the Iterative Robust cluster number and hyperparameters have been chosen maximizing Model-Based Imputation (IRMI) [20] algorithm. mean silhouette width calculated not in the clustering space (i.e. the The Kernelized AutoEncoder To cluster patients considering the whole latent space of the K-AE) but in the space defined by the original data evolution of different variables during the five years of the follow-up of the subjects. Once we obtained the clusters in the latent space using period we used a Kernelized AutoEncoder (K-AE) [16], a deep learning different configurations of the K-AE, clustering algorithm and number model that can learn the representation of MTS mapping them in a of clusters, we divided the study period into ten intervals of 6 months latent space that embeds similarities and differences between patients. each. Then, for each interval we considered only those patients that The MTS inputs that contain data extracted from EHRs of the patients had at least a measure of HbA1c and BMI. We calculated the silhouette are mapped in the latent space by an encoder built on a bidirectional width for the period considering, aside from the two variables just stacked LSTM layer [21]. Then the decoder, composed by another mentioned, also the pharmacological treatment, the comorbidities, the stacked LSTM, maps the vectors in the latent space back in the original age and the diabetes duration in the corresponding period. MTS. Furthermore, vectors in the latent space are aligned to a kernel Clusters were defined separately for men and women. Once we ob- matrix K built from the original MTS. Hence the model is trained with tained the clusters for both sexes, we matched them to pair each cluster a cost function that is the weighted sum of the reconstruction error and of men with a cluster of women. In order to do so, we represented the kernel alignment cost, namely: each cluster with a vector with the mean of all values of the patients in ‖ 𝑍𝑍 𝑇 that cluster, and we calculated the euclidean distance between men and 𝐾 ‖ 𝑛𝑝 1 ∑ ‖ ‖ 𝐿 = 𝛼𝑑 𝑀𝑆𝐸(𝑃𝑖 , 𝑃̂𝑖 ) + 𝛼𝑘 ‖ − ‖ (1) women clusters, pairing clusters according to the minimum distance. 𝑛𝑝 𝑖=1 ‖ ‖𝑍𝑍 𝑇 ‖𝐹 ‖𝐾‖ ‖ ‖ 𝐹 ‖𝐹 To validate our clustering method, we also clusterized patients using Where: the same variables used to create the latent space considering not their evolution but rather a baseline value, similarly to what was done in • 𝛼𝑑 and 𝛼𝑘 are the weights of the reconstruction loss and kernel other works [3,8]. All results were also compared with a random cluster alignment cost functions respectively. as a reference. • 𝑃̂𝑖 is the estimation of the MTS 𝑃𝑖 made by the decoder. The clinical characteristics in each cluster were described for mean • 𝑛𝑝 is the number of patients. and standard deviation for continuous variables and frequencies and • 𝑀𝑆𝐸(𝑃𝑖 , 𝑃̂𝑖 ) is the Mean Square Error between the reconstructed percentages for categorical variables. The evolution of different vari- and the original MTS calculated for each patient for 𝑛𝑡,𝑖 time steps ables for each cluster has been tested using linear mixed models count- available for that patient. ing for sex and age of the patient as a random effect and studying the • 𝑍𝑍 𝑇 ∈ R𝑛𝑝 ×𝑛𝑝 is a matrix containing the dot products of the effect of each cluster on the variable evolution respect to the disease representations in the latent space of the MTS. duration. The percentages of patients diagnosed with comorbidities or • ‖ ⋅ ‖𝐹 represents the Frobenius norm. that received a given drug prescription for each cluster was tested • K is the kernel matrix for the alignment. with the two-proportion z-test (𝛼 = 0.05). The statistical analyses were performed using R3.6.3 software [27]. 2.3. Kernel matrix 3. Results In the original K-AE proposed by Bianchi et al. the kernel matrix was based on the time series cluster kernel [22]. Since this algorithm is 3.1. Clusters validation results based on some assumptions on the data (that are described as Gaussian Mixture Models with a time dependent mean and a constant covariance In order to better visualize results of this work we developed a web and missing values that must be Missing At Random), we decided to app available at dm2.b2slab.upc.edu. use a different kernel matrix based on the less restrictive Fast Global After applying the inclusion criteria, we obtained 11,028 patients, Alignment kernel (GAK) [23]. of which 5226 women and 5802 men. Patients’ age (calculated at Practically, our kernel matrix K was obtained by iterating the origi- 01.01.2013) ranged from 24 to 100 years, with a mean of 69.8±10.6 nal GAK algorithm several times using at each iteration a different set of years. Age at T2DM onset ranged from 18 to 91 years, with a mean of variables to mitigate the possible effects of errors due to the imputation. 59.4±10.9 years. Of these patients, 8790 were still active at the end At each iteration the unnormalized GAK was computed using only a of the study, 1793 died during the analyzed period and the remain- random subset of variables chosen among the 10 clinical variables, ing 445 were transferred from the database. Mean age at death was the diagnosis and the pharmacological treatment (except for the first 81.4±8.7 years. iteration in which all these variables were used) and it is added to the Clusterization and validation have been carried out separately for precedent GAK value. At the end of the iterations K was normalized to sexes, obtaining in both cases an optimal number of clusters equal have 1 on the principal diagonal. This kernel alignment should help the to 7. The relationship found between men and women clusters was Encoder create a latent space that better represents patients’ evolutions bijective, i.e. each cluster of one sex was paired with exactly one cluster and maintains distances between patients. of the other sex and the other way around. For both sexes, the best All these steps have been implemented in Python 3.8.10, using configuration of parameters of the K-AE generated a latent space that Tslearn 0.4.1 [24] for the creation of the kernel matrix and Tensorflow can be used to cluster patients better with respect to using a baseline 2.4.3 [25] to create and define the deep learning model. value for each variable, as reported in Fig. 1. 2.4. Clusters definition and validation 3.2. Clusters description Once we obtained the vectors in the latent space, we divided pa- The main statistics of the seven clusters are reported in Table 1, tients into different clusters testing different state of art algorithms for while Fig. 2 contains the mean values of the different variables such clustering (k-mean, clara, and hierarchical clustering) in this space [26]. as drug usage and comorbidities for the different clusters. Numerical The definition of the optimal number of clusters in a vector space can values for these variables are reported in Supplementary Tables S4, S5 be a delicate issue. This is true especially in this case where a series and S6. 3 E. Manzini et al. Journal of Biomedical Informatics 135 (2022) 104218 Fig. 1. Mean silhouette width evolution. Colors represent different clustering methods: green lines for clusters obtained with our K-AE model; blue lines for clusters obtained with a single basal value; and red lines for random clusterization. Dot and continuous lines represent respectively women and men clusters. K-AE used to cluster women had the following parameters: encoder dim. = 200; decoder dim. = 100; Latent space dim. = 200; 𝛼𝑘 = 0.4. K-AE used to cluster men had the following parameters: encoder dim. = 150; decoder dim. = 200; Latent space dim. = 150; 𝛼𝑘 = 0.3. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 2. Graphical representation of variables, comorbidities and drug usage for the seven clusters (Neuropathic Cluster-NC, Hypercholesteraemic Cluster-HCC, Multiple Complications Cluster-MCC, Vascular Disease Cluster-VDC, Hypertensive Cluster-HTC, Retinopathy Cluster-RC, and Metabolic Cluster-MC) for women and men. Variables plots represent the mean value of all the measures of the corresponding cluster over the analyzed period, comorbidities plots represent percentage of patients in the cluster that have been diagnosed for the corresponding diseases and drugs plots represent the percentage of patients in the cluster that used at least once one medication in the corresponding class. Note that CKD-EPI represents 100-eGFR in order to be coherent with other variables where higher values indicate worse values. (Used abbreviation: Chol. Tot.-Total Cholesterol, Opth.Comp. = Ophthalmological Complications; Neuropt. = Neuropathies; Sulf. = Sulfonylureas, Metfor. = Metformin.). 4 E. Manzini et al. Journal of Biomedical Informatics 135 (2022) 104218 Table 1 Main statistics for the seven clusters (Neuropathic Cluster-NC, Hypercholesteraemic Cluster-HCC, Multiple Complications Cluster-MCC, Vascular Disease Cluster-VDC, Hypertensive Cluster-HTC, Retinopathy Cluster-RC, and Metabolic Cluster-MC). NC (A) HCC (B) MCC (C) VDC (D) HTC (E) RC (F) MC (G) All 832 1472 562 1305 4135 1232 1490 # patients Male 315 831 250 893 2065 724 724 Fem. 517 641 312 412 2070 508 766 All 69.4 ± 10.2 65.3 ± 11.6 69.3 ± 10.3 72.8 ± 9.4 71.1 ± 9.9 68.9 ± 10.5 69.2 ± 10.9 Age1 [y] Male 68 ± 10.4 64.7 ± 11.6 69 ± 9.5 71.4 ± 9.4 69.2 ± 9.9 66.8 ± 10.1 66.7 ± 11.1 Mean±SD2 Fem. 70.3 ± 10.1 66.2 ± 11.5 69.6 ± 11 75.7 ± 8.9 73.0 ± 49.7 72 ± 10.3 71.5 ± 10.1 All 58.7 ± 10.5 56.1 ± 11.5 57.4 ± 11.6 62.2 ± 9.8 61.6 ± 10 56.5 ± 11.2 57.8 ± 11.2 Age at diagnosis [y] Male 57.6 ± 10.4 55.8 ± 11.4 58.1 ± 10.4 61 ± 9.6 60 ± 9.9 54.8 ± 10.5 55.8 ± 11.3 Mean±SD Fem. 59.4 ± 10.5 56.6 ± 11.6 56.9 ± 12.5 64.7 ± 9.9 63.2 ± 9.8 59 ± 11.9 59.8 ± 10.8 All 10.7 ± 5.1 9.2 ± 3.9 11.9 ± 6.5 10.6 ± 4.8 9.5 ± 3.9 12.4 ± 6.7 11.3 ± 5.7 T2DM duration [y] Male 10.4 ± 4.5 8.9 ± 3.4 10.8 ± 5.4 10.4 ± 4.5 9.2 ± 3.6 12 ± 6.2 10.9 ± 5.4 Mean±SD Fem. 10.9 ± 5.4 9.5 ± 4.4 12.7 ± 7.2 11 ± 5.3 9.8 ± 4.1 13 ± 7.2 11.7 ± 5.9 All 14.66 10.87 21 26.9 14.07 17.61 16.31 Death [%] Male 20.32 13.48 22.4 25.98 14.67 17.27 15.33 Fem. 11.22 7.49 19.87 28.88 13.48 18.11 17.23 1 Age calculated on January 1st, 2013. 2 SD: Standard Deviation. Similar mean ages were observed for clusters A, C, F, and G. Patients a metformin prescription at some time during the analysis period. This in cluster B were the youngest in terms of age and age at diagnosis, with percentage is significantly higher than in all other clusters. In contrast, mean values that were statistically lower than all other clusters. The cluster C was the one with the lower percentage of patients treated with mean age and age at diagnosis of clusters D and E were statistically metformin (77.4%). Sulphonylureas were mainly prescribed in cluster E higher than all other clusters but the difference between these two and B (44,6% and 44.4%) with percentages that are statistically similar clusters was not statistically significant. The subjects in cluster F had a to each other, and higher than in all other clusters. In cluster G, the longer period of diabetes duration; however, this value was statistically highest percentages of prescriptions of all other classes of antidiabetic similar to cluster C. In the clusters B and E, the subjects had a shorter drugs (SGLT-2i, DPP-4i and arGLP1) occurred. Subjects in this cluster period of diabetes duration, with values that are statistically similar were also more frequently insulin-treated (97,5%), followed by cluster (pairwise t-test with Benjamini–Hochberg adjustment, 𝛼 = 0.05 [28]). C (70.28%), while subjects in cluster E were the least frequently treated The percentage of deaths was highest in cluster D, followed by cluster with insulin (8.6%). In all cases, the percentages of insulin-treated C, while cluster B was the one with the lowest death percentage. patients were significantly different among all clusters. The results of The clusters showed different clinical profiles, as reported in Figs. 2, the usage of antidiabetic drugs among the clusters are presented in 3 and in Supplementary Figures S2–S10 and Table S4 . Regarding Supplementary Table S5. glycaemic control, the mean HbA1c ranged from 7.0% to 8.0%. Cluster The distribution of common comorbidities among the clusters was G had the worst glycaemic control, while cluster E the best glycaemic also different. Cardiovascular complications were almost entirely rep- control. The variability of HbA1c was different among the clusters. In resented in cluster D (the oldest), significantly higher than in all other cluster C, we observed stable level of HbA1c, while in clusters A, B, D, E clusters. For instance, cardiovascular disease was only found in 2.2% and F we observed upward trends, with markedly different slopes, that of the subjects from clusters B and E. Neuropathy occurred mostly in of cluster B being more than double that of other clusters. In cluster G, cluster A, while in cluster E, the complication only occurred in 4.1% we observed a downward slope over the analyzed period. The majority of the subjects. Renal complications were primarily present in cluster of the subjects in the clusters A, C, E, F and G were obese (BMI>30), C (69.4%), while the percentage in other clusters was much lower, while in clusters B and D subjects were overweight (BMI>25). Cluster ranging from 1.88% in cluster E to 0.2% in cluster B. High blood C had the highest, and cluster B the lowest, mean BMI. This variable pressure was very common in all clusters, except for cluster B: in cluster tended to decrease over time in the clusters A, B, D, E and G; on the E it was present in all patients, while in other clusters it ranged from contrary, it tended to increase for clusters C and F. 91% to 81%, except for cluster B where it was 9.24%. In contrast, Regarding the lipid profile, subjects in the cluster A, B and E were ophthalmological complications were primarily prevalent in cluster characterized by high levels for total cholesterol, LDL and HDL, being F, followed by cluster C, while in other clusters, this complication cluster B the one with the highest values, while triglycerides were occurred less frequently. higher in cluster C and G than the rest of the clusters. Conversely, Due to differential clinical factors, we named cluster A as Neuro- cluster D was characterized by the lowest levels for HDL, LDL and total pathic Cluster (NC), cluster B as Hypercholesteraemic Cluster (HCC), cholesterol. For all clusters LDL and total cholesterol tend to decrease cluster C as Multiple Complications Cluster (MCC), cluster D as Vascular over time, as well as triglycerides, even if in this case this trend is Disease Cluster (VDC), cluster E as Hypertensive Cluster (HTC), cluster significant only for clusters B, C, E and G, with the slope of cluster F as Retinopathy Cluster (RC) and cluster G as Metabolic Cluster (MC). G being more than double of the one of the other clusters. On the contrary, HDL remains stable over time, except for clusters B and G 4. Discussion where it tends to increase. The renal profile among the clusters was also different. The mean glomerular filtration rate (eGFR) in all clusters Results of our unsupervised deep learning cluster analysis suggest ranged from 56.3 to 78.1 ml/min/1.73m2 . The lowest (eGFR) was seven longitudinal phenotype clusters of T2DM persons with different observed in cluster C (which also had the highest ACR and creatinine clinical evolutions. levels); furthermore, in this cluster we observed a higher decrease of the Health care professionals often need to identify, quantify, and in- eGFR. On the contrary, cluster B had the highest eGFR and the lowest terpret relationships among different health variables [29]. Access to ACR and Creatinine, with values that remained stable over time. big real-world clinical data and the advances in artificial intelligence Metformin was the most frequently prescribed drug in all of the applications in recent years has made it possible to develop algorithms clusters, except in cluster G subjects in whom the most frequently that generally are more accurate for the prediction and classification prescribed drug was insulin. In cluster B, 86.6% of the subjects had of patients [29]. Increasingly, deep learning is changing the analysis 5 E. Manzini et al. Journal of Biomedical Informatics 135 (2022) 104218 Fig. 3. HbA1c and BMI evolution calculated over a period of 10 years for the seven clusters (Neuropathic Cluster-NC, Hypercholesteraemic Cluster-HCC, Multiple Complications Cluster-MCC, Vascular Disease Cluster-VDC, Hypertensive Cluster-HTC, Retinopathy Cluster-RC, and Metabolic Cluster-MC). Dot lines represent the mean values while intervals represent lower and upper quantiles computed with 0.95 confidence intervals. of databases extracted from EHRs. It requires less manual feature field is represented by the stratification of T2DM patients in different engineering, while the high volume and coverage of healthcare datasets groups on the basis of different clinical characteristics [33]. To date, enable successful training of complex deep learning models [30]. several studies have been published regarding the possible phenotyping In general, architectures based on AEs have proven to be optimal of T2DM, although none of them clustered patients based on the evo- when dealing with data extracted from EHR. In DeepPatient [13], the lution of the disease. They rather stratified them based on observations authors suggest that a vector of features learned without supervision at a single time (e.g. at diagnosis), possibly analyzing the evolution can be used effectively for several tasks, from disease prediction to of the diseases for a limited period of time after the observations, clustering. In [12], the authors showed that an architecture based on a but only once the patients were stratified. This limits the affidability variational AE can be used to cluster longitudinal data from patients, of the clusters and also their applicability in clinical practice, since resulting in purer clusters than in methods not based on unsupervised these studies are often based on the measures of variables that are deep learning. Similar results were achieved in [31], where a denoising not routinely collected in clinical practice [15]. To the best of our AE was used to derive patients representations to be used for prediction knowledge, this is the first time that the evolution of T2DM has been and study of patient similarity. As for these studies, also our results investigated through longitudinal clusters. Our algorithm, in fact, used suggest that patient representation obtained through our K-AE model a K-AE to map 5 years of routinely collected data from EHR in a latent can be used to obtain better clusters with respect to classical techniques space that embedded the different characteristics and trajectories of not based on DL. patient phenotypes. The diagnosis of T2DM is made on chronically elevated blood One of the most recent works related to T2DM phenotype clustering glucose concentrations and on exclusion of other forms of diabetes was realized by Ahlqvist et al. [3]. Five different clusters were identi- (autoimmunity, pregnancy, pancreatic disease. . . ) [32]. The fact that fied based on six different variables measured at diagnosis, these being this diagnosis process is based only on one consequence of the diseases glutamic acid decarboxylase (GAD) antibodies, age at diagnosis, BMI, (high HbA1c) and on exclusion criteria should admonish us to the lack HbA1c, and homoeostasis model estimates of 𝛽-cell function (HOMA-B) of understanding we have of the diseases causes and on the urge to and insulin resistance (HOMA-IR). According to clustering of subjects develop precision medicine for diabetes [32]. One step further is this based on these phenotypic variables, five subgroups were defined: SAID 6 E. Manzini et al. Journal of Biomedical Informatics 135 (2022) 104218 severe autoimmune diabetes; SIDD, severe insulin-deficient diabetes; 𝛽-cell function. Secondly, data came from the routine clinical practice SIRD, severe insulin-resistant diabetes; MOD, mild obesity-related di- of primary care attended T2DM subjects, hence there was missing data. abetes; MARD, mild age-related diabetes. Even if this study has been To minimize the effect of this, we only filtered subjects who had a confirmed by subsequent works, although with some differences due minimum of two variables measured during the analyzed period and to the variables used [34] or the population [35], a recent study sug- had at least one measure of HbA1C and BMI. It should also be noted gested that in five years from diagnosis, approximately a quarter of the that recent studies have questioned the utility and existence of different patients changed cluster with respect to their initial classification [36]. clusters in the T2DM population [39,40]. Besides the clusters described This result suggests that T2DM is a complex disease that shows different in current work, the main advantage of its technique is to map patient phenotype trajectories that can hardly be described by a cross sectional data into a single vector of the latent space. Distances between different study. vectors embed differences between the patients they represent in terms Our deep learning algorithm suggested 7 typical trajectories inside of the evolution of the disease, obtaining a metric to evaluate patients’ the diabetic population with different evolutions of the disease. The similarities. Such representations could represent a further step toward Multiple Complications and Metabolic clusters have the highest BMI, a precision medicine in diabetes but also in the management of all the high percentage of insulin drugs use (70.3% and 97.5%, respectively) chronic diseases, since it relies only on routinely collected data and and a percentage of Metformin usage lower than in other clusters. can potentially be applied to the study of all chronic diseases. Moreover, in the Multiple Complications Cluster, we observed higher Another limitation of the methodology is the interpretability. Deep percentages of renal complications (69.4%) and the highest levels of learning models can produce accurate predictions; however, the pres- ACR, with signs of macroalbuminuria. In the Metabolic Cluster the ence of the black-box models decreases the interpretability and trans- percentage of renal complications is much lower but still it is the parency of their inner workings. This may be an important issue since second highest one. This is also the only cluster in with both BMI clinicians often are unwilling to accept machine recommendations and HbA1c tend to decrease over time. It is interesting to note that which lack clarity as to their underlying reasoning [30]. In order those two clusters share some similarities with the SIRD cluster (high to make our results clearer, we developed an interactive dashboard BMI, renal complications,..) and that they could represent two pos- showing the evolution of the clusters but also a representation of the sible evolutions of patients that at diagnosis were assigned to this latent space obtained by the K-AE. cluster, with the Metabolic Cluster characterized by a more intense In conclusion, using unsupervised deep learning algorithms on rou- pharmacological treatment. The Vascular Disease and Hypertensive tinely collected health care data from primary care centers in Catalonia, clusters were composed by older subjects (compared with the rest of seven different clusters were obtained in terms of evolution of clinical the clusters at the moment of diagnosis) and both of them showed characteristics, comorbidities, and important outcomes. Although our better glycaemic control with respect to other clusters (mean HbA1c finding confirms the heterogeneity of the T2DM subjects, further stud- 7.3%±1.3 and 7.0%±1.1, respectively). Also in this case, these two clus- ies are needed to provide more clinical and genetically enriched data ters share similarities with the MARD cluster and they show different to perform more precise type 2 diabetes subtype clustering. The de- treatments, with the Vascular Disease cluster characterized by a higher velopment of unsupervised deep machine learning might be a valuable use of insulin. The Hypercholesteraemic cluster was characterized by tool for future applicability in decision-support to help clinicians in the an early age at diagnosis of T2DM, and relatively low BMI, that tends management of T2DM or other chronic diseases. Ultimately these inno- to decrease over time, sharing some similarities with Ahlqvist SAID vative methodologies will help develop precise strategies to implement and SIDD clusters. Their SIDD cluster can also be compared with our and improve disease management, using a more individualized target Retinopathy Cluster, characterized by a similar age at diagnosis but a approach. higher BMI that slightly increases over time with respect to the HC. Moreover, both SIDD and our RC clusters had the highest percentage CRediT authorship contribution statement of ophthalmological complications. Finally, the Neuropathic cluster is similar to the RC in terms of Hba1C and BMI but with a higher age at Enrico Manzini: Conceptualization, Data curation, Methodology, diagnosis and a shorter diabetes duration. Software, Formal analysis, Writing – original draft, Writing – review & Recent guidelines on treatment of diabetes highlighted the im- editing. Bogdan Vlacho: Data curation, Methodology, Writing – origi- portance of personalized treatments for risk management of patients nal draft, Writing – review & editing. Josep Franch-Nadal: Data cura- with T2DM [37]. From this point of view, subclassification of T2DM patients could be a powerful tool for the personalization of diabetic tion, Methodology, Writing – review & editing. Joan Escudero: Writing management [38]. The treatment of different risk factors and compli- – review & editing, Resources. Ana Génova: Writing – review & editing, cations (such as CVD, obesity, CKD, and others) is in fact often decided Resources. Elisenda Reixach: Conceptualization, Resources. Erik An- separately, without discussions on the whole clinical phenotypic evo- drés: Conceptualization, Resources. Israel Pizarro: Conceptualization, lution of diabetes in each individual patient [38]. Our unsupervised Resources. José-Luis Portero: Conceptualization, Resources. Dídac algorithm suggested 7 different longitudinal clusters of patients that Mauricio: Data curation, Methodology, Writing – review & editing. showed different clinical profiles, but that also had different treatments Alexandre Perera-Lluna: Conceptualization, Resources, Methodology, and comorbidities. As mentioned before, for instance, the Metabolic Writing – review & editing. and the Multiple Complication Clusters share common phenotypic char- acteristics such as age at diagnosis, diabetes duration and high levels Declaration of competing interest of HbA1c. However, the pharmacological treatment in the Metabolic Cluster is more intensive, while the rate of patients developing comor- The authors declare the following financial interests/personal rela- bidities (such as CVD and ophthalmological) were much lower. Further tionships which may be considered as potential competing interests: studies will be necessary to evaluate the future clinical applications Enrico Manzini reports was provided by Pulso. of the proposed stratification. The proposed approach, once validated Enrico Manzini, Alexandre Perera-Lluna reports a relationship with and improved, might be a useful instrument to support the clinical Spanish Ministry of Economy and Competitiveness that includes: fund- management of type 2 diabetes. ing grants. The present study has some limitations. Firstly, we only have access Enrico Manzini, Alexandre Perera-Lluna reports a relationship with to a limited amount of health care data and for a limited period of time Networking Biomedical Research Centre in the subject area of Bio- (5 years). Therefore, the clustering would be more precise if possible engineering, Biomaterials and Nanomedicine that includes: funding linkage with genetic data or specific laboratory parameters related to grants. 7 E. Manzini et al. Journal of Biomedical Informatics 135 (2022) 104218 Enrico Manzini, Alexandre-Perera-Lluna reports a relationship with [5] E. Blind, H. Janssen, K. Dunder, P. Graeff, The European Medicines Agency’s Instituto de Investigación Carlos III that includes: funding grants. approval of new medicines for type 2 diabetes, Diabetes Obes. Metab. 20 (2017) http://dx.doi.org/10.1111/dom.13349. Enrico Manzini, Alexandre Perera-Lluna reports a relationship with [6] M. Mata-Cases, D. Mauricio, I. Vinagre, R. Morros, E. Hermosilla, F. Fina, M. Share4Rare project that includes: funding grants. Josep Franch-Nadal, Rosell Murphy, C. Castell, J. Franch-Nadal, B. Bolíbar, Treatment of hyper- Didac Mauricio reports a relationship with AstraZeneca Pharmaceuti- glycaemia in type 2 diabetic patients in a primary care population database cals LP that includes: funding grants and speaking and lecture fees. in a Mediterranean Area (Catalonia, Spain), Diabetes Metab. 5 (2014) 338, Josep Franch-Nadal, Didac Mauricio reports a relationship with http://dx.doi.org/10.4172/2155-6156.1000338. Ascensia Diabetes Care pain L that includes: speaking and lecture fees. [7] S.G. Nicholls, S.M. Langan, E.I. Benchimol, Routinely collected data: The im- portance of high-quality diagnostic coding to research, CMAJ 189 (33) (2017) Josep Franch-Nadal, Didac Mauricio reports a relationship with E1054–E1055, http://dx.doi.org/10.1503/cmaj.170807, URL https://www.cmaj. Boehringer Ingelheim GmbH that includes: funding grants and speaking ca/content/189/33/E1054. and lecture fees. [8] J.V. Stidsen, J.E. Henriksen, M.H. Olsen, R.W. Thomsen, J.S. Nielsen, J. Josep Franch-Nadal, Didac Mauricio reports a relationship with G K Rungby, S.P. Ulrichsen, K. Berencsi, J.A. Kahlert, S.G. Friborg, I. Brandslund, that includes: funding grants and speaking and lecture fees. A.A. Nielsen, J.S. Christiansen, H.T. Sørensen, T.B. Olesen, H. Beck-Nielsen, Pathophysiology-based phenotyping in type 2 diabetes: A clinical classifi- Josep Franch-Nadal, Didac Mauricio reports a relationship with Lilly cation tool, Diabetes/Metab. Res. Rev. 34 (5) (2018) e3005, http://dx.doi. pain that includes: funding grants and speaking and lecture fees. org/10.1002/dmrr.3005, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/ Josep Franch-Nadal, Didac Mauricio reports a relationship with M dmrr.3005, e3005 DMRR-17-RES-100.R2. D that includes: funding grants and speaking and lecture fees. [9] W. Oh, E. Kim, M.R. Castro, P.J. Caraballo, V. Kumar, M.S. Steinbach, G.J. Simon, Josep Franch-Nadal, Didac Mauricio reports a relationship with Type 2 diabetes mellitus trajectories and associated risks, Big Data 4 (1) (2016) 25–30, http://dx.doi.org/10.1089/big.2015.0029. Novartis Pharmaceuticals Corporation that includes: funding grants and [10] I. Walraven, M.R. Mast, T. Hoekstra, A. Jansen, A.A. van der Heijden, S.P. speaking and lecture fees. Rauh, F. Rutters, E. van’t Riet, P.J. Elders, A.C. Moll, et al., Distinct HbA1c Josep Franch-Nadal, Didac Mauricio reports a relationship with trajectories in a type 2 diabetes cohort, Acta Diabetol. 52 (2) (2015) 267–275, Novo Nordisk Inc that includes: funding grants and speaking and http://dx.doi.org/10.1007/s00592-014-0633-8. lecture fees. [11] R.A. Hubbard, J. Xu, R. Siegel, Y. Chen, I. Eneli, Studying pediatric health Josep Franch-Nadal, Didac Mauricio reports a relationship with outcomes with electronic health records using Bayesian clustering and trajectory analysis, J. Biomed. Inform. 113 (2021) 103654, http://dx.doi.org/10.1016/j.jbi. Sanofi that includes: funding grants and speaking and lecture fees. 2020.103654. [12] J. de Jong, M.A. Emon, P. Wu, R. Karki, M. Sood, P. Godard, A. Ahmad, Acknowledgments H. Vrooman, M. Hofmann-Apitius, H. Fröhlich, Deep learning for clustering of multivariate clinical patient trajectories with missing values, GigaScience 8 (11) This work was supported by the Spanish Ministry of Economy and (2019) 1–14, http://dx.doi.org/10.1093/gigascience/giz134. Competitiveness (www.mineco.gob.es) TEC2014-60337-R, DPI2017- [13] R. Miotto, L. Li, B.A. Kidd, J.T. Dudley, Deep patient: An unsupervised represen- tation to predict the future of patients from the electronic health records, Sci. 89827-R, Networking Biomedical Research Centre in the subject area of Rep. 6 (May) (2016) 1–10, http://dx.doi.org/10.1038/srep26094. Bioengineering, Biomaterials and Nanomedicine (CIBER-BBN), Spain, [14] F. Xie, H. Yuan, Y. Ning, M.E.H. Ong, M. Feng, W. Hsu, B. Chakraborty, N. Liu, initiatives of Instituto de Investigación Carlos III (ISCIII), and Share4Rare Deep learning for temporal data representation in electronic health records: A project (Grant Agreement 780262). This study was possible thanks to systematic review of challenges and methodologies, J. Biomed. Inform. (2021) the commitment of physicians and nurses working in the Catalan Health 103980, http://dx.doi.org/10.1016/j.jbi.2021.103980. Institute to provide optimal care to patients with diabetes. CIBER of [15] W. Chung, K. Erion, J. Florez, A. Hattersley, m.-f. Hivert, C. Lee, M. McCarthy, J. Nolan, J. Norris, E. Pearson, L. Philipson, A. McElvaine, W. Cefalu, S. Rich, Diabetes and Associated Metabolic Diseases (CIBERDEM) is an initiative P. Franks, Precision medicine in diabetes: A consensus report from the American from Instituto de Salud Carlos III, Madrid, Spain. Diabetes Association (ADA) and the European Association for the Study of This analysis is part of the DiaCare Project of Novo Nordisk and the Diabetes (EASD), Diabetes Care 43 (7) (2020) 1617–1635, http://dx.doi.org/10. Fundació TicSalut (Departament de Salut, Generalitat de Catalunya), in 2337/dci20-0022. collaboration with Grupo Pulso, for the benefit of people with type 2 [16] F.M. Bianchi, L. Livi, K. Mikalsen, M. Kampffmeyer, R. Jenssen, Learning representations of multivariate time series with missing data, Pattern Recognit. diabetes. 96 (2019) 106973, http://dx.doi.org/10.1016/j.patcog.2019.106973. [17] B. Bolíbar, F. Fina Avilés, R. Morros, M.d.M. Garcia-Gil, E. Hermosilla, R. Ramos, Appendix A. Supplementary data M. Rosell, J. Rodríguez, M. Medina, S. Calero, D. Prieto-Alhambra, Grupo SIDIAP, [SIDIAP database: Electronic clinical records in primary care as a source of Supplementary material related to this article can be found online information for epidemiologic research], Med. Clin. 138 (14) (2012) 617–621, http://dx.doi.org/10.1016/j.medcli.2012.01.020, URL http://www.ncbi.nlm.nih. at https://doi.org/10.1016/j.jbi.2022.104218. gov/pubmed/22444996. [18] G. Brämer, International statistical classification of diseases and related health References problems. Tenth revision, World Health Stat. Q. Rapport Trimestriel de Statis- tiques Sanitaires Mondiales 41 (1) (1988) 32—36, URL http://europepmc.org/ [1] R. Busse, M. Blümel, D. Scheller-Kreinsen, A. Zentner, Tackling chronic disease abstract/MED/3376487. in Europe. Strategies, interventions and challenges, Political Science, Vol. 2009, [19] ATC classification index with DDDs, 2021, WHO Collaborating Centre for Drug 2010. Statistics Methodology, Oslo, Norway. [2] P. Saeedi, I. Petersohn, P. Salpea, B. Malanda, S. Karuranga, N. Unwin, S. [20] M. Templ, A. Kowarik, P. Filzmoser, Iterative stepwise regression imputation Colagiuri, L. Guariguata, A.A. Motala, K. Ogurtsova, J.E. Shaw, D. Bright, using standard and robust methods, Comput. Statist. Data Anal. 55 (10) (2011) R. Williams, Global and regional diabetes prevalence estimates for 2019 and 2793–2806, http://dx.doi.org/10.1016/j.csda.2011.04.012. projections for 2030 and 2045: Results from the international diabetes federation [21] S. Hochreiter, J. Schmidhuber, Long short-term memory, Neural Comput. 9 diabetes atlas, 9th edition, Diabetes Res. Clin. Pract. 157 (2019) http://dx.doi. (1997) 1735–1780, http://dx.doi.org/10.1162/neco.1997.9.8.1735. org/10.1016/j.diabres.2019.107843. [3] E. Ahlqvist, P. Storm, A. Käräjämäki, M. Martinell, M. Dorkhan, A. Carlsson, P. [22] K.Ø. Mikalsen, F.M. Bianchi, C. Soguero-Ruiz, R. Jenssen, Time series cluster Vikman, R.B. Prasad, D.M. Aly, P. Almgren, Y. Wessman, N. Shaat, P. Spégel, kernel for learning similarities between multivariate time series with missing H. Mulder, E. Lindholm, O. Melander, O. Hansson, U. Malmqvist, A. Lernmark, data, Pattern Recognit. 76 (2018) 569–581, http://dx.doi.org/10.1016/j.patcog. K. Lahti, T. Forsén, T. Tuomi, A.H. Rosengren, L. Groop, Novel subgroups of 2017.11.030. adult-onset diabetes and their association with outcomes: A data-driven cluster [23] M. Cuturi, Fast global alignment kernels, in: Proceedings of the 28th International analysis of six variables, Lancet Diabetes Endocrinol. 6 (5) (2018) 361–369, Conference on Machine Learning, ICML 2011, 2011. http://dx.doi.org/10.1016/S2213-8587(18)30051-2. [24] R. Tavenard, J. Faouzi, G. Vandewiele, F. Divo, G. Androz, C. Holtz, M. Payne, [4] A. Bonnefond, P. Froguel, Clustering for a better prediction of type 2 diabetes R. Yurchak, M. Rußwurm, K. Kolar, E. Woods, Tslearn, A machine learning mellitus, Nat. Rev. Endocrinol. 17 (4) (2021) http://dx.doi.org/10.1038/s41574- toolkit for time series data, J. Mach. Learn. Res. 21 (118) (2020) 1–6, URL 021-00475-4. http://jmlr.org/papers/v21/20-091.html. 8 E. Manzini et al. Journal of Biomedical Informatics 135 (2022) 104218 [25] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G.S. Corrado, [35] X. fen Xiong, Y. Yang, L. Wei, Y. Xiao, L. Li, L. Sun, Identification of two A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, novel subgroups in patients with diabetes mellitus and their association with M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, clinical outcomes: A two-step cluster analysis, J. Diabetes Investig. 7 (2021) R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, http://dx.doi.org/10.1126/scitranslmed.aaa9364. I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. [36] O.P. Zaharia, K. Strassburger, A. Strom, G.J. Bönhof, Y. Karusheva, S. Antoniou, Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, X. Zheng, TensorFlow: K. Bódis, D.F. Markgraf, V. Burkart, K. Müssig, J.H. Hwang, O. Asplund, L. Groop, Large-scale machine learning on heterogeneous systems, 2015, URL https://www. E. Ahlqvist, J. Seissler, P. Nawroth, S. Kopf, S.M. Schmid, M. Stumvoll, A.F. tensorflow.org/, Software available from tensorflow.org. Pfeiffer, S. Kabisch, S. Tselmin, H.U. Häring, D. Ziegler, O. Kuss, J. Szendroedi, [26] M. Maechler, P. Rousseeuw, A. Struyf, M. Hubert, K. Hornik, Cluster: Cluster M. Roden, B.F. Belgardt, A. Buyken, J. Eckel, G. Geerling, H. Al-Hasani, C. analysis basics and extensions, 2021, URL https://CRAN.R-project.org/package= Herder, A. Icks, J. Kotzka, E. Lammert, D. Markgraf, W. Rathmann, Risk of cluster, R package version 2.1.2. diabetes-associated diseases in subgroups of patients with recent-onset diabetes: [27] R. Core Team, R: A language and environment for statistical computing, 2020, A 5-year follow-up study, Lancet Diabetes Endocrinol. 7 (9) (2019) 684–694, R Foundation for Statistical Computing, Vienna, Austria, URL https://www.R- http://dx.doi.org/10.1016/S2213-8587(19)30187-1. project.org/. [37] J.B. Buse, D.J. Wexler, A. Tsapas, P. Rossing, G. Mingrone, C. Mathieu, D.A. [28] Y. Benjamini, Y. Hochberg, Controlling the false discovery rate: A practical and D’Alessio, M.J. Davies, 2019 Update to: management of hyperglycemia in type 2 powerful approach to multiple testing, J. R. Stat. Soc. Ser. B Stat. Methodol. 57 diabetes, 2018. a consensus report by the American diabetes association (ADA) (1) (1995) http://dx.doi.org/10.1111/j.2517-6161.1995.tb02031.x. and the European Association for the Study of Diabetes (EASD), Diabetes Care [29] K. Johnson, J. Soto, B. Glicksberg, S. Khader, R. Miotto, M. Ali, E. Ashley, 43 (2) (2020) 487–493. Artificial intelligence in cardiology, J. Am. Coll. Cardiol. 71 (2018) 2668–2679, [38] H. Tanabe, H. Masuzaki, M. Shimabukuro, Novel strategies for glycaemic control http://dx.doi.org/10.1016/j.jacc.2018.03.521. and preventing diabetic complications applying the clustering-based classification [30] C. Xiao, E. Choi, J. Sun, Opportunities and challenges in developing deep of adult-onset diabetes mellitus: A perspective, Diabetes Res. Clin. Pract. 180 learning models using electronic health records data: A systematic review, J. (2021) 109067, http://dx.doi.org/10.1016/j.diabres.2021.109067. Am. Med. Inform. Assoc. 25 (10) (2018) 1419–1428, http://dx.doi.org/10.1093/ [39] M. Lugner, S. Gudbjörnsdottir, N. Sattar, A. Svensson, M. Miftaraj, K. Eeg- jamia/ocy068. Olofsson, B. Eliasson, S. Franzén, Comparison between data-driven clusters and [31] L. Lei, Y. Zhou, J. Zhai, L. Zhang, Z. Fang, P. He, J. Gao, An effective patient models based on clinical features to predict outcomes in type 2 diabetes: representation learning for time-series prediction tasks based on EHRs, in: 2018 Nationwide observational study, Diabetologia (2021) 1–9, http://dx.doi.org/10. IEEE International Conference on Bioinformatics and Biomedicine, BIBM, 2019, 1007/s00125-021-05485-5. http://dx.doi.org/10.1109/BIBM.2018.8621542. [40] J.M. Dennis, B.M. Shields, W.E. Henley, A.G. Jones, A.T. Hattersley, Disease [32] H. Fitipaldi, M.I. McCarthy, J.C. Florez, P.W. Franks, A global overview of progression and treatment response in data-driven subgroups of type 2 diabetes precision medicine in type 2 diabetes, Diabetes 67 (10) (2018) 1911–1922, compared with models based on simple clinical features: An analysis using http://dx.doi.org/10.2337/dbi17-0045. clinical trial data, Lancet Diabetes Endocrinol. 7 (6) (2019) 442–451, http: [33] A.L. Gloyn, D.J. Drucker, Precision medicine in the management of type 2 //dx.doi.org/10.1016/S2213-8587(19)30087-7. diabetes, Lancet Diabetes Endocrinol. 6 (11) (2018) 891–900, http://dx.doi.org/ 10.1016/S2213-8587(18)30052-4. [34] R.C. Slieker, L.A. Donnelly, H. Fitipaldi, G.A. Bouland, G.N. Giordano, M. Åkerlund, M.J. Gerl, E. Ahlqvist, A. Ali, I. Dragan, et al., Replication and cross-validation of type 2 diabetes subtypes based on clinical variables: An IMI- RHAPSODY study, Diabetologia 64 (9) (2021) 1982–1989, http://dx.doi.org/10. 1007/s00125-021-05490-8. 9