A Structural FHMM for Interpretable Disease Trajectories in T2DM

# A Structural FHMM for Interpretable Disease Trajectories in T2DM
Source: [https://arxiv.org/html/2608.24328](https://arxiv.org/html/2608.24328)
Alessandro MariAffiliation:Swiss Data Science Center \(SDSC\), Ecole Polytechnique Fédérale de Lausanne \(EPFL\), SwitzerlandGuillaume ObozinskiAffiliation:Swiss Data Science Center \(SDSC\), ETH Zürich, SwitzerlandMaria Luisa Marques de Sa FaquettiAffiliation:Pharmacoepidemiology Group, Institute of Pharmaceutical Sciences, Zurich, SwitzerlandAdrian Martinez de la TorreAffiliation:Pharmacoepidemiology Group, Institute of Pharmaceutical Sciences, Zurich, SwitzerlandAndrea BurdenAffiliation:Pharmacoepidemiology Group, Institute of Pharmaceutical Sciences, Zurich, Switzerland

August 24, 2026

###### Abstract

In this work, we propose a structural variant of the Factorial Hidden Markov Model \(FHMM\) for the analysis of disease trajectories in patients with Type 2 diabetes mellitus \(T2DM\)\. The model represents a patient’s latent health state as a combination of multiple independent, simultaneously evolving components, associated with comorbidities and lab results\. This structured latent representation facilitates the identification of clinically meaningful patient states and clustering of common disease trajectories\. We evaluate the proposed approach using The IQVIA Medical Research Data incorporating data from THIN, a Cegedim database of anonymized electronic health records \(EHR\), identifying patients with a first\-ever prescription for a non\-insulin antidiabetic drug \(NIAD\) between January 2006 and December 2019\. The model identifies multiple clinically coherent latent components corresponding to known patterns of diabetes\-related complications and reveals heterogeneous progression pathways, including distinct microvascular\-dominant and multi\-organ trajectories associated with elevated comorbidity burden and mortality\. These results demonstrate that the proposed framework captures meaningful longitudinal structure in EHR data and provides interpretable insights into the evolution of T2DM and its comorbidities\.

## 1Background

Type 2 diabetes mellitus \(T2DM\) is a chronic metabolic disorder characterized by elevated blood glucose levels resulting from insulin resistance and inadequate insulin secretion\. It represents a major global health concern and is frequently accompanied by a range of comorbidities, which substantially worsen quality of life and increase the risk of morbidity and mortality\[[1](https://arxiv.org/html/2608.24328#bib.bib2)\]\. The presence and progression of such comorbidities contribute to the heterogeneity of T2DM trajectories, making disease management and risk stratification particularly challenging\. The inference of a patient’s underlying health state from longitudinal clinical observations, together with information on T2DM and its associated comorbidities, can support the identification of disease progression patterns and emerging health risks\. Leveraging retrospective electronic health records \(EHRs\) to characterize these evolving health states offers the potential to improve disease understanding, stratify patients into clinically meaningful subgroups, and inform timely interventions\.

In this study, we analyze longitudinal observations from EHRs associated with T2DM and multiple comorbidities in order to infer patients’ latent health states and characterize common patterns of disease evolution\. To this end, we propose a statistical model designed to capture the multi\-dimensional and dynamic nature of patient health, with a particular focus on interpretability of the latent representation and its clinical relevance\.

Specifically, we consider a structural variant of the Factorial Hidden Markov Model \(FHMM\)\[[2](https://arxiv.org/html/2608.24328#bib.bib1)\], which enables modeling the dynamics of a multi\-dimensional patient health state through several independent latent components evolving simultaneously over time\. Given longitudinal clinical observations, the proposed model infers the latent configuration of these components, providing a structured representation of disease progression and associated health risks\. In contrast to a standard Hidden Markov Model \(HMM\), which represents patient health using a single latent state, the FHMM framework allows the patient’s overall condition to be approximated as a combination of a small number of evolving latent processes\. This structure increases modeling flexibility and facilitates interpretation, as individual latent chains can be associated with distinct aspects of disease evolution or comorbidity burden\.

By using retrospective information on patients’ health indicators and disease history, the proposed model enables the identification of clinically meaningful latent states and the analysis of transitions between them\. The resulting latent trajectories can be used to cluster patients according to common progression patterns, characterize typical disease pathways, and assess the evolution of comorbidity\-related risks over time\.

Related work\.Tracking and modeling comorbidities from longitudinal health observations in patients with T2DM is crucial for understanding disease heterogeneity and supporting effective disease management\. A wide range of statistical and machine learning approaches have been proposed to analyze large\-scale health data and uncover temporal patterns associated with disease progression and risk\. HMM and their generalizations are among the most widely used probabilistic approaches for modeling disease dynamics in healthcare applications\[[3](https://arxiv.org/html/2608.24328#bib.bib10),[4](https://arxiv.org/html/2608.24328#bib.bib12),[5](https://arxiv.org/html/2608.24328#bib.bib4),[6](https://arxiv.org/html/2608.24328#bib.bib13)\]\. Extensions such as mixtures of Input–Output Hidden Markov Models have been proposed to capture more complex temporal dependencies\[[7](https://arxiv.org/html/2608.24328#bib.bib11)\]\. Nonparametric latent feature models, including approaches based on the Indian Buffet Process, have also been explored for comorbidity analysis, particularly in the context of psychiatric disorders\[[8](https://arxiv.org/html/2608.24328#bib.bib5)\]\.

More recently, deep generative models have been applied to longitudinal clinical data\. In\[[9](https://arxiv.org/html/2608.24328#bib.bib8)\], the authors propose a variational autoencoder\-based framework to learn latent representations that jointly generate clinical measurements and medical concepts\. While effective, their approach requires manual assignment of latent dimensions to specific medical concepts\. In contrast, our model automatically learns a structured factorization of the latent space by construction\. Moreover, the latent representation in our approach is discrete and subject to monotonicity constraints, which enhances interpretability and facilitates the analysis of disease trajectories\.

Methods addressing the joint evolution of multiple diseases have also been proposed\. For example,\[[10](https://arxiv.org/html/2608.24328#bib.bib3)\]introduces a bivariate copula\-based hidden process to model the interaction between two diseases using population\-level time series\. In the context of T2DM, machine learning approaches have been applied to model disease progression\[[11](https://arxiv.org/html/2608.24328#bib.bib7)\], while clustering\-based analyses have been used to study patterns of comorbidities\[[12](https://arxiv.org/html/2608.24328#bib.bib9)\]\. Our work complements these approaches by providing a structured, interpretable latent\-state model adapted to the analysis of complex, multi\-comorbidity disease trajectories\.

## 2Methods

### 2\.1Data

The IQVIA Medical Research Database UK \(IMRD\-UK\) incorporates data supplied by The Health Improvement Network \(THIN\), a Cegedim database of anonymized electronic health records generated from the daily records of General Practitioners \(GPs\)\. The study protocol and use of the data were reviewed and approved by the Scientific Review Committee \(SRC, reference number 20SR062\)\.

The dataset includes patients from participating primary care practices who have provided general consent for their data to be used for research purposes\. Consent may be withdrawn at any time, in which case the individual’s data are removed from the dataset\. Secondary care providers, such as ambulatory care or specialists, provide additional feedback on diagnoses to the GPs\.

Generally, the database comprises routinely collected data from clinical practice, including patient diagnoses, symptoms, prescriptions, referrals, immunizations, laboratory test results \(including glycated hemoglobin \[HbA1c\]\), selected behavioral factors \(body mass index \[BMI\], alcohol use, and smoking status\), and demographic characteristics\. As with most electronic healthcare databases, data are generated during routine care rather than for research purposes, resulting in unequally spaced observations and incomplete capture of certain variables\. Lifestyle factors such as BMI, smoking status, and alcohol consumption are subject to missingness\. In addition, diagnoses are recorded only when they are among the main reasons for a visit, laboratory tests are captured only when ordered; and while medication prescriptions are identified, prescription fills and adherence are not known\.

### 2\.2Cohort description

We identified all adult patients \(aged 18\+\) with a first\-ever prescription of a non\-insulin antidiabetic drug \(NIAD\) between January 1, 2006, and December 31, 2019\. The index date was defined as the date of the first NIAD prescription\. To ensure identification of new users, patients were required to have at least one year of valid data prior to the index date\. We excluded patients without any laboratory measurement for variables of interest \(see Table[1](https://arxiv.org/html/2608.24328#S2.T1)\) within the 4\-years prior to the index date after data aggregation \(see Section[2\.2\.1](https://arxiv.org/html/2608.24328#S2.SS2.SSS1)\)\.

#### 2\.2\.1Data preprocessing

Figure 1:Graphical representation of the preprocessing steps\.To harmonize the dataset and prepare it for modeling, we applied the following data preparation steps\. A summary of these steps is provided in Figure[1](https://arxiv.org/html/2608.24328#S2.F1)\.

##### Comorbidities and Laboratory Measurements:

The comorbidities and laboratory variables included in the analysis are listed in Table[1](https://arxiv.org/html/2608.24328#S2.T1)\. Comorbidities were identified at or prior to the index date and during follow\-up using read codes recorded at each patient’s clinical encounters\. Only chronic conditions were considered\. Accordingly, a comorbidityccwas defined as present at timettif any corresponding read code appeared during the relevant six\-month interval\. Once identified, comorbidities were assumed to persist over time, reflecting their chronic or recurrent nature; thus, the indicator for comorbidityccremained set to 1 following its first occurrence\. In the base model, drugs were only used in modeling the initial state\. A detailed discussion of the effects of drugs used in the transition matrix is provided in Appendix[C](https://arxiv.org/html/2608.24328#A3)\.

The laboratory values are standardized to consistent units \(see Table[1](https://arxiv.org/html/2608.24328#S2.T1)\), and outliers are removed using predefined thresholds \(see Appendix[A](https://arxiv.org/html/2608.24328#A1)\)\. To address irregular measurement timing and partial missingness inherent to routine clinical data, laboratory measurements were aggregated into six\-month intervals, with the median value used when multiple measurements were available within a given interval\. Patients were excluded if any required laboratory variable had no recorded measurements during the four years preceding NIAD initiation\. This step was necessary to enable forward imputation of missing values during the post\-NIAD period\.

Additionally, to separate the timing of doctors’ decisions regarding diagnoses and medication prescriptions, we used the comorbidity status \(presence or absence\) and lab tests shifted forward in time relative to the prescriptions\. This approach helps avoid issues related to six\-month aggregation, where diagnoses might be influenced by medication prescriptions rather than the patient’s actual health state\.

To further preserve the intended temporal ordering between patient health status and clinical decision\-making, comorbidity indicators and laboratory measurements were shifted forward in time relative to medication prescriptions\. This approach mitigates potential reverse\-causation artifacts introduced by six\-month aggregation, whereby diagnoses or laboratory results recorded shortly before any treatment initiation could otherwise appear contemporaneous with \(or influenced by\) drug exposure\.

##### Clinical Stratification:

Patient observations were stratified by age and sex to account for clinically meaningful heterogeneity\. Age was determined at the index date and categorized into two groups: patients aged 64 years or younger and those aged 65 years or older\. This stratification reflects differences in diabetes progression and the burden of age\-related comorbidities across the life course\[[13](https://arxiv.org/html/2608.24328#bib.bib24)\]\. Additionally, analyses were stratified by sex to capture known differences in disease prevalence and comorbidity profiles, such as osteoporosis, between male and female patients\[[14](https://arxiv.org/html/2608.24328#bib.bib23)\]\.

ComorbiditiesLab ValuesDrugs \(with ATC code\)•Stroke•Myocardial infarction \(MI\)•Congestive heart failure \(CHF\)•Peripheral vascular disease \(PVD\)•Diabetic retinopathy \(DR\)•Hypertension \(HT\)•Chronic kidney disease \(CKD\)•Osteoporosis•Chronic obstructive pulmonary disease \(COPD\)•Chronic liver disease \(CLD\)•Glycated hemoglobin \(Hb A1c\) in %•Low\-density Lipoprotein \(LDL\) cholesterol in mmol/L•High\-density Lipoprotein \(HDL\) cholesterol in mmol/L•Triglycerides \(TG\) in mmol/L•Blood pressure \(Systolic/Diastolic\) in mmHg•Body Mass Index \(BMI\) in kg/m2•Glomerular filtration rate \(GFR\) in mL/min/1\.73m2•Diabetes Drugs111Including Gliclazide, Glimepiride, Glipizide, Metformin, Pioglitazone, Rosiglitazone, Sitagliptin, Dapagliflozin\.•Antithrombotic agents \(B01\)•Cardiac therapy \(C01\)•Antihypertensives \(C02\)•Diuretics \(C03\)•Beta blocking agents \(C07\)•Calcium channels blockers \(C08\)•ACE inhibitors \(C09\)•Statins \(C10\)•Corticosteroids \(H02\)•Drugs for treatment of bone diseases \(M05\)

Table 1:List of variables that are included in the model\.

### 2\.3Structural FHMM

Our modeling approach is based on Factorial Hidden Markov Model \(FHMM\)\[[2](https://arxiv.org/html/2608.24328#bib.bib1)\], which is an extension of the hidden Markov chain, where the output is conditioned not on one hidden Markov chain withMMstates, but onKKindependent hidden Markov chains\. This allows to model the time series depending on several independent latent processes \(see Figure[2](https://arxiv.org/html/2608.24328#S2.F2)\)\. In terms of comorbidities and lab\-values progression modeling, the FHMM approach translates into the assumption that there existKKhidden health state components, which are independent and that are jointly influencing the output\. Although the independence assumption may appear unrealistic, since body systems rarely function in isolation, the FHMM can be understood as an independent basis decomposition of the dynamics, where several hidden processes evolve separately and jointly generate the observed comorbidities and lab values\. This view highlights FHMM as a structured approximation: the latent chains act as abstract basis processes that capture key aspects of disease progression\.

X1X\_\{1\}X2X\_\{2\}X3X\_\{3\}H1\(1\)H\_\{1\}^\{\(1\)\}H2\(1\)H\_\{2\}^\{\(1\)\}H3\(1\)H\_\{3\}^\{\(1\)\}H1\(2\)H\_\{1\}^\{\(2\)\}H2\(2\)H\_\{2\}^\{\(2\)\}H3\(2\)H\_\{3\}^\{\(2\)\}H1\(3\)H\_\{1\}^\{\(3\)\}H2\(3\)H\_\{2\}^\{\(3\)\}H3\(3\)H\_\{3\}^\{\(3\)\}Y1LABY\_\{1\}^\{\\text\{LAB\}\}Y2LABY\_\{2\}^\{\\text\{LAB\}\}Y3LABY\_\{3\}^\{\\text\{LAB\}\}Y1COMY\_\{1\}^\{\\text\{COM\}\}Y2COMY\_\{2\}^\{\\text\{COM\}\}Y3COMY\_\{3\}^\{\\text\{COM\}\}H1H\_\{1\}H2H\_\{2\}H3H\_\{3\}Y1Y\_\{1\}Y2Y\_\{2\}Y3Y\_\{3\}

Figure 2:Directed graphical model of a FHMM withM=3M=3hidden chains\. Arrows represent causal relationships and gray nodes are observed variables\. White nodes are not observed\. The variablesXtX\_\{t\}are the drugs taken at timettthat influence the patient health stateHt=def\(Ht\(1\),…,Ht\(M\)\)H\_\{t\}\\stackrel\{\{\\scriptstyle\\rm def\}\}\{\{=\}\}\(H\_\{t\}^\{\(1\)\},\\dots,H\_\{t\}^\{\(M\)\}\)\. The results of the lab testsYtLABY\_\{t\}^\{\\text\{LAB\}\}and the set of comorbiditiesYtCOMY\_\{t\}^\{\\text\{COM\}\}are direct consequences of the health stateHtH\_\{t\}at timett\.We assume that the evolution of comorbidity states and laboratory test measurements, denoted byYtY\_\{t\}, follows the generative process defined by

H1\(m\)\|X1\\displaystyle H\_\{1\}^\{\(m\)\}\\mid X\_\{1\}∼π\(m\)\(⋅∣X1;θπ\),\\displaystyle\\sim\\pi^\{\(m\)\}\(\\cdot\\mid X\_\{1\};\{\\theta^\{\\pi\}\}\),m=1,…,M,\\displaystyle m=1,\\dots,M,\(1\)Ht\(m\)\|Ht−1\(m\),Xt\\displaystyle H\_\{t\}^\{\(m\)\}\\mid H\_\{t\-1\}^\{\(m\)\},X\_\{t\}∼P\(m\)\(⋅∣Ht−1\(m\),Xt;θP\),\\displaystyle\\sim P^\{\(m\)\}\(\\cdot\\mid H\_\{t\-1\}^\{\(m\)\},X\_\{t\};\{\\theta^\{P\}\}\),m=1,…,M,\\displaystyle m=1,\\dots,M,\(2\)Yt\|Ht\\displaystyle Y\_\{t\}\\mid H\_\{t\}∼p\(⋅∣Ht;η\),\\displaystyle\\sim p\(\\cdot\\mid H\_\{t\};\\eta\),\(3\)whereHt=def\(Ht\(1\),…,Ht\(M\)\)H\_\{t\}\\stackrel\{\{\\scriptstyle\\rm def\}\}\{\{=\}\}\(H\_\{t\}^\{\(1\)\},\\dots,H\_\{t\}^\{\(M\)\}\)denotes the joint latent configuration ofMMparallel hidden Markov chains withKKpossible states each\. At each time point, the combination of theseMMlatent variables, each taking one of K\-states, determines the distribution of the observableYtY\_\{t\}\.

In the presence of the auxillary input variables, such as in our case medications, one can model the output assuming the dependence of the hidden states on the input, which leads to Input\-output FHMM formulation\.

Exact inference for \([3](https://arxiv.org/html/2608.24328#S2.E3)\) can, in principle, be carried out by representing the FHMM as a single HMM withKMK^\{M\}states, where the transition matrix can be written as

P=P\(1\)⊗P\(2\)⊗⋯⊗P\(M\)∈ℝKM×ℝKM,P=P^\{\(1\)\}\\otimes P^\{\(2\)\}\\otimes\\cdots\\otimes P^\{\(M\)\}\\in\\mathbb\{R\}^\{K^\{M\}\}\\times\\mathbb\{R\}^\{K^\{M\}\},where⊗\\otimesis the Kronecker product\. Nevertheless, the introduction of additional chains makes the estimation intractable, as a naive implementation of the Forward\-Backward algorithm has complexityO⁡\(K2​M​T\)O\(K^\{2M\}T\)\. The complexity can be reduced toO⁡\(M​KM\+1​T\)O\(MK^\{M\+1\}T\)by carefully computing the matrix\-vector dot product\[[15](https://arxiv.org/html/2608.24328#bib.bib14)\], but still grows exponentially in the number of chains\. One of the ways to make the estimation tractable is to use variational approximation, which we describe in Section[2\.4](https://arxiv.org/html/2608.24328#S2.SS4)\.

We impose the following problem\-specific structure on the FHMM to improve identifiability and interpretability\.

#### 2\.3\.1Problem\-specific structure

The transition matrix specifies the probabilities of moving between latent health states from one time step to the next\. The general FHMM defined in Section[2\.3](https://arxiv.org/html/2608.24328#S2.SS3)is not easily interpretable as the states of the chaines are not ordered and any permutation of the states labels\{1,…,K\}\\\{1,\.\.\.,K\\\}gives the same model\. In order to make the identification of the states possible, we assume that states are ordered by worsening of the patient’s condition, and that lower state values represent a healthy patient\. Additionally, we assume that the health of a patient can degrade only byddunits \(a jump ofddstates\) and that the last state is an absorbing state, representing the impossibility of improving the patient’s health condition\. These assumptions mean that the transition matrixPθPP\_\{\\theta^\{P\}\}is a band matrix with bandwidthdd, where the last row of the matrixPθPP\_\{\\theta^\{P\}\}contains all zeros except in the last column, which is equal to one\. Figure[3](https://arxiv.org/html/2608.24328#S2.F3)shows all possible transitions whend=1d=1\. The transition matrix is fully parameterized without the inclusion of input variables, such as medications\. A discussion on incorporating medications into the parameterization is provided in Appendix[C](https://arxiv.org/html/2608.24328#A3)\.

12345Figure 3:Possible transitions between states of a hidden chain with55states and jumps of sized=1d=1\.The Emission probabilityp\(⋅\|Ht;η\)p\(\\cdot\|H\_\{t\};\\eta\)should capture that the higher the patients’ state value is, the worse their health conditions are\. In particular, the conditional mean

𝔼⁡\[Yt\|Ht\]=defg⁡\(η0\+∑mH~t\(m\)⊤​η\(m\)\),\\mathbb\{E\}\[Y\_\{t\}\|H\_\{t\}\]\\stackrel\{\{\\scriptstyle\\rm def\}\}\{\{=\}\}g\\left\(\\eta\_\{0\}\+\\sum\_\{m\}\{\\ \\tilde\{H\}\_\{t\}^\{\(m\)\}\}^\{\\top\}\\eta^\{\(m\)\}\\right\),should be a monotone function of each state componentHt\(m\)H\_\{t\}^\{\(m\)\}for allmm, whereH~t\(m\)∈\{0,1\}K\\tilde\{H\}\_\{t\}^\{\(m\)\}\\in\\\{0,1\\\}^\{K\}is a one\-hot encoded version ofHt\(m\)∈\{1,…,K\}H\_\{t\}^\{\(m\)\}\\in\\\{1,\\dots,K\\\}andggis some link function\. Note thatη1\(m\)=0\\eta\_\{1\}^\{\(m\)\}=0for identifiability issues and that the intercept is captured byη0\\eta\_\{0\}\.

As in\[[16](https://arxiv.org/html/2608.24328#bib.bib18)\], to guarantee monotonicity we parametrize the variableη\(m\)\\eta^\{\(m\)\}as the cumulative sum ofη\+\(m\)≥0\\eta\_\{\+\}^\{\(m\)\}\\geq 0andη−\(m\)≥0\\eta\_\{\-\}^\{\(m\)\}\\geq 0, i\.e\.,

η\(m\)=defD​η\+\(m\)−D​η−\(m\),\\eta^\{\(m\)\}\\stackrel\{\{\\scriptstyle\\rm def\}\}\{\{=\}\}D\\eta\_\{\+\}^\{\(m\)\}\-D\\eta\_\{\-\}^\{\(m\)\},whereDDis a lower triangular matrix with all ones\. This means that each difference between successive coefficients is either non\-negative or non\-positive: ifη−\(m\)=0\\eta\_\{\-\}^\{\(m\)\}=0, then

ηi\(m\)−ηi−1\(m\)=η\+,i\(m\)≥0⇔ηi\(m\)≥ηi−1\(m\),\\eta\_\{i\}^\{\(m\)\}\-\\eta\_\{i\-1\}^\{\(m\)\}=\\eta\_\{\+,i\}^\{\(m\)\}\\geq 0\\iff\\eta\_\{i\}^\{\(m\)\}\\geq\\eta\_\{i\-1\}^\{\(m\)\},henceη\(m\)\\eta^\{\(m\)\}is non\-decreasing\. Similarly, ifη\+\(m\)=0\\eta\_\{\+\}^\{\(m\)\}=0, then

ηi\(m\)−ηi−1\(m\)=−η−,i\(m\)≤0⇔ηi\(m\)≤ηi−1\(m\),\\eta\_\{i\}^\{\(m\)\}\-\\eta\_\{i\-1\}^\{\(m\)\}=\-\\eta\_\{\-,i\}^\{\(m\)\}\\leq 0\\iff\\eta\_\{i\}^\{\(m\)\}\\leq\\eta\_\{i\-1\}^\{\(m\)\},henceη\(m\)\\eta^\{\(m\)\}is non\-increasing\. To decide whether the slope should be positive or negative, we add a group\-lasso penalty\[[17](https://arxiv.org/html/2608.24328#bib.bib17)\]:

p⁡\(η\)=λ×∑m\(‖η−\(m\)‖2\+‖η−\(m\)‖2\)\.p\(\\eta\)=\\lambda\\times\\sum\_\{m\}\\left\(\\\|\\eta\_\{\-\}^\{\(m\)\}\\\|\_\{2\}\+\\\|\\eta\_\{\-\}^\{\(m\)\}\\\|\_\{2\}\\right\)\.\(4\)This penalty enforces sparsity at the group level, i\.e\. it encourages eitherη\+\(m\)=0\\eta\_\{\+\}^\{\(m\)\}=0\(purely decreasing slope\) orη−\(m\)=0\\eta\_\{\-\}^\{\(m\)\}=0\(purely increasing slope\)\.

The initial state is modeled as a functionπθπ\\pi\_\{\\theta^\{\\pi\}\}of the list of drugs taken over a period of 4 years prior to NIAD, age, smoking states and alcohol status jointly denoted byX1X\_\{1\}\. The distribution is assumed to be unimodal and is parametrized by a gamma distribution, i\.e\.,

\[πθπ​\(X1\)\]i=Softmax​\{\(α⁡\(X1\)−1\)​log⁡i−βinit​i\}∝iα⁡\(X1\)−1​exp⁡\(−βinit​i\),\\left\[\\pi\_\{\\theta^\{\\pi\}\}\(X\_\{1\}\)\\right\]\_\{i\}=\\text\{Softmax\}\\left\\\{\(\\alpha\(X\_\{1\}\)\-1\)\\log i\-\\beta^\{\\text\{init\}\}i\\right\\\}\\propto i^\{\\alpha\(X\_\{1\}\)\-1\}\\exp\(\-\\beta^\{\\text\{init\}\}i\),\(5\)withα⁡\(X1\)=Winit​X1\+αinit\\alpha\(X\_\{1\}\)=W^\{\\text\{init\}\}X\_\{1\}\+\\alpha^\{\\text\{init\}\}andθπ=\{Winit,αinit,βinit\}\{\\theta^\{\\pi\}\}=\\\{W^\{\\text\{init\}\},\\alpha^\{\\text\{init\}\},\\beta^\{\\text\{init\}\}\\\}\. All parameters are assumed to be non\-negative\. Ifα\\alphais close to zero, the probability mass is concentrated mostly around the lowest state value, whereas larger values ofα\\alphaincrease the probability of starting from a higher state value\.

Assuming that a drug prescription is equivalent to a degradation of the patient’s health, we imposeWinit≥0W^\{\\text\{init\}\}\\geq 0andαinit≥0\\alpha^\{\\text\{init\}\}\\geq 0\. In addition,WinitW^\{\\text\{init\}\},αinit\{\\alpha^\{\\text\{init\}\}\}andβinit\{\\beta^\{\\text\{init\}\}\}follow a priori a \(truncated\) Gaussian distribution centered at 0, 1 and 1 respectively, which is equivalent to add the penalizations‖Winit‖22\\\|W^\{\\text\{init\}\}\\\|\_\{2\}^\{2\},‖αinit−1‖22\\\|\{\\alpha^\{\\text\{init\}\}\}\-1\\\|\_\{2\}^\{2\}and‖βinit−1‖22\\\|\{\\beta^\{\\text\{init\}\}\}\-1\\\|\_\{2\}^\{2\}to the objective function\. The underlying assumptions are

1. 1\.the drug effect is a priori small \(i\.e\.,WinitW^\{\\text\{init\}\}is likely to be small\),
2. 2\.in the absence of drugs, the likelihood of starting from a higher state is small, that is, the prior distribution is\[πθπ​\(X1\)\]i∝exp⁡\(−i\)\.\\left\[\\pi\_\{\\theta^\{\\pi\}\}\(X\_\{1\}\)\\right\]\_\{i\}\\propto\\exp\(\-i\)\.

### 2\.4Model estimation

Parameters of the models are estimated via Expectation\-Maximization \(EM\) following the work of Ghahramani\[[2](https://arxiv.org/html/2608.24328#bib.bib1)\], with the aim to maximize the Evidence Lower Bound \(ELBO\)

maxθ,qF\(q,θ\)=maxθ,q𝔼h1:T∼q\[logpθ\(y1:T,h1:T∣x1:T\)q\(h1:T∣x1:T,y1:T\)\],\\max\_\{\\theta,q\}F\(q,\\theta\)=\\max\_\{\\theta,q\}\\mathbb\{E\}\_\{h\_\{1:T\}\\sim q\}\\left\[\\log\\frac\{p\_\{\\theta\}\(y\_\{1:T\},h\_\{1:T\}\\mid x\_\{1:T\}\)\}\{q\(h\_\{1:T\}\\mid x\_\{1:T\},y\_\{1:T\}\)\}\\right\],\(6\)over a class of approximative distributionsq∈𝒟q\\in\\mathcal\{D\}\(E step\) and the parametersθ=\(θπ,θP,η\)\\theta=\(\\theta^\{\\pi\},\\theta^\{P\},\\eta\)\(M step\), where we used the notationsa:bs\_\{a:b\}to denote the sequencessfromaatobb\. In\[[2](https://arxiv.org/html/2608.24328#bib.bib1)\], the authors proposed to use a structured mean field variational approximation\[[18](https://arxiv.org/html/2608.24328#bib.bib15)\]

q\(h1:T∣x1:T\)=def∏mq\(h1:T\(m\)∣x1:T,y1:T\),q\(h\_\{1:T\}\\mid x\_\{1:T\}\)\\stackrel\{\{\\scriptstyle\\rm def\}\}\{\{=\}\}\\prod\_\{m\}q\(h\_\{1:T\}^\{\(m\)\}\\mid x\_\{1:T\},y\_\{1:T\}\),with

q\(H1\(m\)=i∣x1,y1:T\)∝e1,i\(m\)×πi\(m\),q\(Ht\(m\)=j∣Ht−1\(m\)=i,xt,y1:T\)∝et,j\(m\)×Pi,j\(m\),\\displaystyle\\begin\{split\}q\\left\(H\_\{1\}^\{\(m\)\}=i\\mid x\_\{1\},y\_\{1:T\}\\right\)&\\propto e\_\{1,i\}^\{\(m\)\}\\times\\pi\_\{i\}^\{\(m\)\},\\\\ q\\left\(H\_\{t\}^\{\(m\)\}=j\\mid H\_\{t\-1\}^\{\(m\)\}=i,x\_\{t\},y\_\{1:T\}\\right\)&\\propto e\_\{t,j\}^\{\(m\)\}\\times P\_\{i,j\}^\{\(m\)\},\\end\{split\}\(7\)as an approximation of the optimal distributionq\(h1:T∣x1:T,y1:T\)=p\(h1:T∣y1:T,x1:T\)q\(h\_\{1:T\}\\mid x\_\{1:T\},y\_\{1:T\}\)=p\(h\_\{1:T\}\\mid y\_\{1:T\},x\_\{1:T\}\), which reduces the inference complexity toO⁡\(K2​M​T\)O\(K^\{2\}MT\)in the case of Gaussian emissions\. The parametersϕ=\{et,j\(m\)\}t,j,m\\phi=\\\{e\_\{t,j\}^\{\(m\)\}\\\}\_\{t,j,m\}are learned by optimizing the ELBO \([6](https://arxiv.org/html/2608.24328#S2.E6)\)\.

In the case of binary outputs, such as comorbidities, the complexity still grows exponentially in the number of chainsO⁡\(M​KM\+1​T\)O\(MK^\{M\+1\}T\)\[[15](https://arxiv.org/html/2608.24328#bib.bib14)\]\. However, inspired by the Gaussian case, a second order approximation of the emission log\-likelihoodlogpη\(yt∣ht\(1:M\)\)\\log p\_\{\\eta\}\(y\_\{t\}\\mid\{h\_\{t\}^\{\(1:M\)\}\}\)is derived using the inequality\[[19](https://arxiv.org/html/2608.24328#bib.bib16)\]

−log⁡\(1\+exp⁡y\)≥−log⁡\{1\+exp⁡\(−z\)\}−\(y\+z\)/2−λ⁡\(z\)2​\(y2−z2\),\-\\log\(1\+\\exp y\)\\geq\-\\log\\left\\\{1\+\\exp\(\-z\)\\right\\\}\-\(y\+z\)/2\-\\frac\{\\lambda\(z\)\}\{2\}\(y^\{2\}\-z^\{2\}\),\(8\)withλ⁡\(z\)=tanh​\(z/2\)2​z\\lambda\(z\)=\\frac\{\\text\{tanh\}\(z/2\)\}\{2z\}and some additional variablez∈ℝz\\in\\mathbb\{R\}\. Exploiting this inequality reduces the complexity of the model training toO⁡\(K2​M​T\)O\(K^\{2\}MT\)\(see Appendix[B](https://arxiv.org/html/2608.24328#A2)for more details\)\. After training, parametersθ\\thetaare kept fixed and the parametersϕ\\phiare learned when presented with new data\. We refer this step as ”fine\-tuning”\.

### 2\.5Trajectory analysis

The emission parametersη\(m\)\\eta^\{\(m\)\}model the risk of developing a comorbidity for a given state, but do not capture the progression of comorbidities over time\. On the other hand, the sequence of hidden states\{Ht\(m\)\}t\>0\\\{H\_\{t\}^\{\(m\)\}\\\}\_\{t\>0\}encodes the temporal dynamics of disease evolution, allowing us to infer how patients transition between latent health states and how these transitions influence the development and accumulation of comorbidities across time\.

In particular, the variational posterior distributionqϕ\(Ht\(m\)∣x1:T,y1:T\)q\_\{\\phi\}\(H\_\{t\}^\{\(m\)\}\\mid x\_\{1:T\},y\_\{1:T\}\)is given by \([7](https://arxiv.org/html/2608.24328#S2.E7)\), after fitting the parametersϕ\\phifor each patient\. The latter is computed using the Baum–Welch algorithm for HMMs \(see\[[2](https://arxiv.org/html/2608.24328#bib.bib1)\]\)\. From the posterior distribution, we calculate the expected state of each latent chain at each time step as

𝔼q\[ht\(m\)\]=∑hthtqϕ\(m\)\(ht∣x1:T,y1:T\),\\mathbb\{E\}\_\{q\}\[h\_\{t\}^\{\(m\)\}\]=\\sum\_\{h\_\{t\}\}h\_\{t\}q\_\{\\phi\}^\{\(m\)\}\(h\_\{t\}\\mid x\_\{1:T\},y\_\{1:T\}\),yielding anMM\-dimensional representation per time point\. We then apply K\-means clustering to all time steps, selecting the number of clusters based on the Davies–Bouldin index and prediction strength\[[20](https://arxiv.org/html/2608.24328#bib.bib21),[21](https://arxiv.org/html/2608.24328#bib.bib22)\]\.

### 2\.6Experimental setup

Patient data were partitioned into three distinct subsets: training \(80%\), validation \(10%\), and test \(10%\)\. Stratification was applied to ensure a consistent distribution of comorbidities across these subsets \(see Table[6](https://arxiv.org/html/2608.24328#A1.T6)in the Appendix\)\.

All models were trained for a maximum of 300 EM iterations, with early stopping applied if the objective functionFFshowed no improvement\. Model parameters were initialized to zero, except forαinit\\alpha^\{\\text\{init\}\},βinit\\beta^\{\\text\{init\}\}, which were initialized to one\. Multiple FHMMs with varying numbers of chains and states were trained, and model selection was performed using the Akaike Information Criterion \(AIC\)\[[22](https://arxiv.org/html/2608.24328#bib.bib19)\]and the Bayesian Information Criterion \(BIC\)\[[23](https://arxiv.org/html/2608.24328#bib.bib20)\]\.

## 3Results

### 3\.1Patient cohort and characteristics

After applying the preprocessing workflow described in Section[2\.2\.1](https://arxiv.org/html/2608.24328#S2.SS2.SSS1), the final cohort comprised 88201 patients, including 47983 men and 40218 women\. Figure[4](https://arxiv.org/html/2608.24328#S3.F4)illustrates the filtering process of patients based on missing laboratory data during the four years preceding NIAD initiation\. In total, 100325 patients were excluded due to missing data for key laboratory values\.

![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/cohort_selection_process.png)Figure 4:Filtering of patients based on missing laboratory data during the four years preceding NIAD initiation\. A total of 100325 patients were excluded\. Abbreviations are defined in Table[1](https://arxiv.org/html/2608.24328#S2.T1)\.Table 2:Patient characteristics at the index date \(NIAD initiation\) for the cohort after preprocessing\. Comorbidities are reported as percentages of the cohort population\. Laboratory measurements are summarized as mean values, with the proportion of missing observations reported in parentheses\. Abbreviations are defined in Table[1](https://arxiv.org/html/2608.24328#S2.T1)\. SD denotes standard deviation and IQR denotes interquartile range\.The median follow\-up duration after NIAD initiation was approximately 4\.0 years across all groups\. Baseline patient characteristics at the index date \(first NIAD prescription\) are summarized in Table[2](https://arxiv.org/html/2608.24328#S3.T2)\. Continuous variables, including age, body mass index \(BMI\), and glycated hemoglobin \(HbA1c\), are reported as medians with interquartile ranges \(IQRs\), while categorical variables—such as sex, smoking status, and comorbidities—are presented as counts and percentages\.

Smoking prevalence was higher among men than women\. Hypertension and diabetic retinopathy were the most prevalent comorbidities at baseline, with diabetic retinopathy showing the greatest increase in prevalence following NIAD initiation\. Comorbidities that developed after NIAD initiation are indicated in parentheses in Table[2](https://arxiv.org/html/2608.24328#S3.T2)\.

### 3\.2Model selection

Table[3](https://arxiv.org/html/2608.24328#S3.T3)presents the AIC and BIC scores for various models evaluated on the validation set\. We use the notationFHMM​\(M,K\)d\\text\{FHMM\}\(M,K\)\_\{d\}to denote a factorial hidden Markov model \(FHMM\) withMMchains,KKstates per chain, and a transition matrix bandwidth ofdd\. Similarly,HMMd\\text\{HMM\}\_\{d\}refers to a standard hidden Markov model \(HMM\) withKKstates and a transition matrixPθPP\_\{\\theta^\{P\}\}of bandwidthdd\.

Table 3:Comparison between models with various number of chains\. Metrics are evaluated on the validation set and drugs are not used as input\.For comparison, an HMM withK=27=33K=27=3^\{3\}states is evaluated against an FHMM of equivalent size\. In the FHMM, each chain has 2 possible transitions \(d=1d=1\), which corresponds to an HMM bandwidth ofd=8=23d=8=2^\{3\}\. Because the relationship between HMM state values and worsening health conditions is less direct \(where a single disease progression may be too restrictive\), the HMM emission probabilities are parameterized without monotonicity constraints, and drug information is excluded from the transition matrix input\.

Memory constraints limit our comparison to a state dimension ofK~=K3=27\\tilde\{K\}=K^\{3\}=27, withK=3K=3\. In the HMM, the full matrix

ξt=𝔼⁡\[H~t​H~t−1⊤\]∈ℝK3×K3\\xi\_\{t\}=\\mathbb\{E\}\\left\[\\tilde\{H\}\_\{t\}\\tilde\{H\}\_\{t\-1\}^\{\\top\}\\right\]\\in\\mathbb\{R\}^\{K^\{3\}\\times K^\{3\}\}must be computed for all patients and timesteps, whereH~t\\tilde\{H\}\_\{t\}is the one\-hot encoding ofHtH\_\{t\}\. In contrast, the FHMM primarily requires computing

ξt\(m\)=𝔼⁡\[H~t\(m\)​H~t−1\(m\)⊤\]∈ℝK×K\\xi\_\{t\}^\{\(m\)\}=\\mathbb\{E\}\\left\[\\tilde\{H\}\_\{t\}^\{\(m\)\}\{\\tilde\{H\}\_\{t\-1\}^\{\(m\)\}\\ \}^\{\\top\}\\right\]\\in\\mathbb\{R\}^\{K\\times K\}for all patients, timestepstt, and chainsmm\. Since memory complexity scales quadratically withKK, we keptKKsmall\. In particular, increasing the number of states from 4 to 6 did not necessarily improve the BIC on the validation set for larger models \(see Table[3](https://arxiv.org/html/2608.24328#S3.T3)\)\.

The top\-performing model wasFHMM​\(8,4\)2\\text\{FHMM\}\(8,4\)\_\{2\}across all sex and age groups, except among*Men<65<65*, for whomFHMM​\(6,4\)2\\text\{FHMM\}\(6,4\)\_\{2\}performed better according to both BIC and AIC\. The modelFHMM​\(8,4\)2\\text\{FHMM\}\(8,4\)\_\{2\}corresponds to the largest latent space, capable of encoding up to48=65,5364^\{8\}=65,536possible states\. In what follows, we will focus on the best performing model for each age and sex group\. Generally, FHMM models consistently outperform HMM models in terms of BIC\.

### 3\.3Model analysis

Next we provide the analysis of the results for*Women<65<65*subgroup\. Other subgroups can be found in the Appendix\. The model selected on validation was ”fine\-tuned” on the testing set, meaning that the variational parametersϕ\\phiwere newly inferred based on the latter\.

![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female_m65_8_4_bin_no_p.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female_m65_8_4_cont_no_p.png)

Figure 5:Top:Emission parametersη\\etamodeling the risk of developing a comorbidity for a given state\. Larger values indicate a higher probability of having the comorbidity\. The values are reported on a logarithmic scale\.Bottom:Mean value of lab tests with respect to the intercept \(i\.e\., state 0\)\. Red \(resp\. blue\) indicates an increase \(resp\. decrease\) in𝔼⁡\[Yt\|Ht\]\\mathbb\{E\}\[Y\_\{t\}\|H\_\{t\}\]relative to the intercept\. The model includes 8 chains with 4 states per chain and is fitted on onWomen<65<65\.![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female_m65_8_4_init_no_p.png)Figure 6:Estimated contribution of drug prescriptions to initial latent state allocation \(WinitW^\{\\text\{init\}\}\), based on a model with 8 chains and 4 states per chain fitted to the*Women<65<65*subgroup\. Higher values indicate a greater probability of starting in a more severe latent state\. These estimates capture baseline correlations, not causal effects, since they describe state assignment rather than counterfactual deviation\.##### Emission parameters\.

The emission parameters for each latent state are shown in Figure[5](https://arxiv.org/html/2608.24328#S3.F5)\. By construction \(Section[2\.3\.1](https://arxiv.org/html/2608.24328#S2.SS3.SSS1)\), higher states correspond to greater probabilities of specific comorbidities, revealing clusters of conditions that tend to co‐occur:

- •Chain 0, termed ”Osteoporosis progression,” combines elevated risks of osteoporosis, chronic kidney disease and fatty liver disease with a high‐HDL/LDL and low‐triglyceride profile\.
- •Chain 1, ”Cardiovascular progression,” captures macrovascular outcomes \(peripheral vascular disease, heart failure, death\) and renal failure \(low GFR\), along with osteoporosis\.
- •Chain 2, ”Infarction progression,” is marked by myocardial infarction risk accompanied by BMI reduction\.
- •Chain 3, ”Metabolic syndrome progression,” features obesity, dyslipidemia, chronic obstructive pulmonary disease, liver disease and death risk\.
- •Chain 4, ”Hypertension progression,” reflects isolated elevation of systolic and diastolic blood pressure\.
- •Chain 5, ”Retinopathy progression,” shows increased diabetic retinopathy probability alongside poor glycemic and lipid control \(elevated HbA1c, triglycerides and BMI, reduced HDL\)\.
- •Chain 6, ”Vascular progression,” combines stroke, peripheral vascular disease and COPD with a comparatively favorable lipid profile\.
- •Chain 7, ”Cardiorenal progression,” integrates myocardial infarction, stroke, chronic heart failure, chronic kidney disease and dyslipidemia with obesity\. Hypertension and diabetic retinopathy tend to develop independently of other comorbidities in our model\.

##### Initial latent state\.

The relationship between the initial latent state and the input variables is shown in Figure[6](https://arxiv.org/html/2608.24328#S3.F6)\. Positive values \(in red\) correspond to a higher probability of beginning in a more severe state\. Smoking shows a marked positive value in chain 3, reflecting its association with increased COPD risk\. The age parameter is positive across several chains, indicating that older patients tend to start in worse health states\. This pattern is not observed in chain 5, where age has no significant effect on the initial state\. Receiving treatment for bone diseases is associated with an increased likelihood of being in the osteoporosis state at the start\.

Table 4:Patient characteristics in each cluster for theWomen<65<65cohort with average values and standard deviation where applicable\. Comorbidities are presented as absolute counts at NIAD initiation, along with the number of new comorbidities developed post\-NIAD\. Proportions, expressed as percentages of the cohort population, are reported in parentheses\.
##### State Clustering\.

Following Section[2\.5](https://arxiv.org/html/2608.24328#S2.SS5), we identified 13 clusters in the latent representation according to the selection criterion \(see Figure[12](https://arxiv.org/html/2608.24328#A4.F12)in the Appendix\)\. Table[4](https://arxiv.org/html/2608.24328#S3.T4)details the baseline characteristics of each cluster at NIAD initiation, and Figure[7](https://arxiv.org/html/2608.24328#S3.F7)presents their time\-averaged clinical profiles\.

![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female_m65/proba_avg.png)Figure 7:Average cluster\-specific probabilities and lab values forWomen<65<65\. Values have been centered by subtracting the average at the time of NIAD initiation to improve readability\. Red bars indicate increases, and blue bars indicate decreases relative to the average patient at NIAD\. Comorbidity risks are shown as changes in probability \(left axis\), and lab values as percentage changes relative to the reference value \(right axis\)\. Annotated numbers indicate absolute changes in the original units \(see Table[1](https://arxiv.org/html/2608.24328#S2.T1)\)\. Error bars represent the standard error of the mean\.Several clusters demonstrate similarities at NIAD initiation\. Clusters 0, 1, and 11 share characteristics of younger patients \(age<50<50\) with low comorbidity burden\. However, their clinical trajectories differ significantly \- particularly cluster 1, which despite its younger profile shows higher smoking prevalence, elevated BMI, and increased mortality risk \(see Table[4](https://arxiv.org/html/2608.24328#S3.T4)\)\.

The clusters can be divided based on their mortality risk\. At one end of the spectrum, clusters 0, 2, 4, 5, and 11 form a lower\-risk group characterized by relatively low death rates\. An intermediate category includes clusters 1, 9, and 10, which show progressing disease with moderate mortality risk\. The highest\-risk group comprises clusters 3, 6, 7, 8, and 12, showing substantially elevated mortality rates that indicate advanced disease states\. Within this high\-risk category, cluster 7 is particularly notable: it is defined by the highest BMI \(morbid obesity with BMI\>40\>40\) and carries an exceptionally high comorbidity burden, indicating particularly severe disease progression\.

The healthiest subgroup, cluster 0 \(”young low\-risk”\), is characterized by preserved metabolic and renal function and low cardiovascular risk\. Cluster 3 \(”pulmonary\-metabolic”\) combines metabolic dysfunction and pulmonary disease with high smoking prevalence and renal impairment\. This cluster’s profile, with increased risk in nearly all comorbidities, positions it as a transition toward severe multimorbidity\. On the other hand, cluster 10 \(”bone\-vascular”\) features high osteoporosis prevalence and increased peripheral vascular disease, suggesting shared pathophysiological pathways between skeletal and vascular systems\.

Among the most severe phenotypes, cluster 8 \(”advanced nephropathy”\) shows severe GFR loss, persistent hypertension, and high mortality\. Cluster 12 \(”cardio\-pulmonary multimorbidity”\) is defined by COPD, CHF, and elevated stroke risk\. Cluster 6 \(”hepatorenal syndrome”\) represents a very severe situation, linking chronic liver disease, renal dysfunction \(low GFR\), and poor glycemic control, consistent with fatty liver disease trends\[[24](https://arxiv.org/html/2608.24328#bib.bib28)\]\.

![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female_m65/distribution.png)Figure 8:Distribution of the clusters over time\.
##### State progression\.

The temporal distribution of the clusters is shown in Figure[8](https://arxiv.org/html/2608.24328#S3.F8)\. Clusters associated with advanced or progressing disease \(e\.g\., cluster 7 and transitional clusters 2, 6, 10, 12\) gain members over time, while healthier clusters \(0, 1, 4\) decline\.

![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female_m65/cumulative_incidence.png)Figure 9:Cumulative Incidence Function \(CIF\) for the time to first cluster transition inWomen<65<65, stratified by initial cluster assignment at 6 months post\-NIAD\. The estimates account for censoring due to loss to follow\-up and study termination, computed using the Aalen nonparametric estimator\[[25](https://arxiv.org/html/2608.24328#bib.bib6)\]\.Figure[9](https://arxiv.org/html/2608.24328#S3.F9)presents the cumulative incidence functions for time to first cluster transition, stratified by initial cluster assignment\. The transition dynamics are further illustrated in the appendix \(Figure[13](https://arxiv.org/html/2608.24328#A4.F13)\)\. Patients in cluster 0 often progress into early microvascular states, particularly cluster 11 or cluster 4, defined by higher prevalence of diabetic retinopathy and hypertension, respectively\. After approximately four years, these states may progress to early multimorbidity \(cluster 5\) or a combination of both diseases \(cluster 2\), as confirmed by Figure[9](https://arxiv.org/html/2608.24328#S3.F9)showing most cluster 4 patients transitioning to cluster 2\. In contrast, cluster 1, despite representing younger patients, exhibits higher mortality risk and demonstrates progression to advanced disease states\. These patients typically evolve into clusters 6 or 7 in later stages, indicating more severe disease trajectories\.

Cluster 10 \(bone\-vascular\) represents an intermediate severity profile\. Patients in this cluster tend to progress toward broader multimorbidity through two primary pathways: either via clusters 2 or 5, or occasionally involving hepatic complications \(cluster 6\)\. The advanced disease profiles represented by cluster 12 \(cardio\-pulmonary multimorbidity\) typically transition to cluster 7, the most severe multimorbid state, that serves as a frequent absorbing state\. Similarly, clusters 2, 5, and 11 demonstrate low transition probabilities, functioning as stable states where patients tend to remain over time \(at least when starting in those clusters\)\.

The transition graph in Figure[10](https://arxiv.org/html/2608.24328#S3.F10)\(showing only transitions with less than 1% probability\) reveals that clusters 0, 2, 4, 5, 10, and 11 form a transient subgraph\. This configuration indicates that once patients transition out of these clusters, they rarely return, suggesting a unidirectional progression toward more severe disease states\. The observed pattern reflects the clinical reality that as patients develop serious comorbidities, their condition tends to evolve toward increasingly complex and severe multimorbidity profiles rather than reverting to less severe states\.

![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female_m65/graph_avg.png)Figure 10:Transition graph\. Red arrows \(resp\. blue\) means an increase \(resp\. decrease\) of the comorbidity risk \(Women<65<65\)\. Arrow thickness is proportional to the number of observed transitions\. Transitions with less than 1% probability are not shown\.

## 4Discussion

We identified eight latent components that capture distinct patterns of comorbidity in type 2 diabetes in women under 65 years of age \(with additional population analyses provided in the Appendix\)\. At the interpretation level, the inferred latent chains correspond to recognizable diabetes\-related complications, spanning microvascular, macrovascular, metabolic, skeletal, pulmonary, and cardiorenal involvement\. Several components reflect well\-established comorbidity patterns in T2DM; for example, chains dominated by macrovascular and renal outcomes align with the known coupling between cardiovascular disease and chronic kidney disease in diabetes, often described as a cardiorenal syndrome\[[26](https://arxiv.org/html/2608.24328#bib.bib29)\]\. Similarly, components emphasizing microvascular complications, such as retinopathy and hypertension, appear to develop independently of other complications and follow expected clinical patterns, including associations with elevated blood pressure and poor glycemic and lipid control, respectively\[[27](https://arxiv.org/html/2608.24328#bib.bib30),[28](https://arxiv.org/html/2608.24328#bib.bib25)\]\. The skeletal involvement in our osteoporosis\-related component co\-occurred with chronic kidney and liver disease, reflecting the clustering of multisystem complications in T2DM\. This pattern is consistent with the understanding that diabetes affects multiple organ systems through interconnected metabolic and vascular pathways\[[29](https://arxiv.org/html/2608.24328#bib.bib26),[30](https://arxiv.org/html/2608.24328#bib.bib27)\]\.

While many of these associations are individually well documented, our model learns them jointly from longitudinal data without prespecifying disease groupings or assuming independence between outcomes\. This yields an interpretable decomposition of multimorbidity into parallel disease processes that may co\-occur and progress at different rates within individuals\.

At the patient level, clustering of latent\-state trajectories highlights how these components combine into progression pathways\. Patients often transition from relatively healthy or low\-risk metabolic states \(clusters 0, 2, 4, 5, 10, 11\) toward advanced multimorbid states, with clusters 6 and 7 representing terminal disease profiles associated with elevated mortality\. These trajectories broadly separate into at least two modes of progression: one dominated by microvascular complications \(clusters 2, 4, 5, 10, 11\), and another driven by pulmonary, hepatic, renal, and cardio\-pulmonary involvement \(clusters 3, 6, 8, 9, 12\)\.

From a methodological perspective, our model reveals patterns that are difficult to extract with alternative approaches\. While static clustering approaches summarize patient similarity but do not encode temporal dynamics\[[12](https://arxiv.org/html/2608.24328#bib.bib9)\], and deep generative models can learn flexible latent representations from longitudinal clinical data but often require additional assumptions or post hoc alignment to relate latent dimensions to clinical concepts\[[9](https://arxiv.org/html/2608.24328#bib.bib8)\], our approach provides a structured factorization that enables direct clinical interpretation\. Classical HMM\-based models capture temporal dynamics but often rely on a single latent process, potentially mixing distinct disease mechanisms\[[3](https://arxiv.org/html/2608.24328#bib.bib10),[4](https://arxiv.org/html/2608.24328#bib.bib12)\]\. The factorial structure adopted here enables multimorbidity to be represented as the superposition of multiple latent chains, providing a natural framework for chronic diseases characterized by overlapping and asynchronous complications\.

This study has several limitations\. Real\-world electronic health records are collected based on routine clinical care, resulting in unequally spaced measurements that often required forward imputation to maintain complete information at each time\-point\. Some variables, particularly BMI and other lifestyle variables, exhibit missingness that is likely not missing at random\. While forward imputation avoids using future information, alternative imputation schemes might be more appropriate for handling missing values\. The current model uses Gaussian emissions for scalability, though extensions to other distributions \(e\.g\., count\-based likelihoods\) are possible\. The data were grouped into binned time steps, imposing an assumption of regular observation intervals that does not reflect the irregular nature of real\-world measurements\. The analysis was restricted to a relatively small set of commonly occurring chronic comorbidities and laboratory values, excluding less frequent conditions that might contribute to multimorbidity patterns\. Issues of data imbalance persist, as many comorbidities represent rare events, and treatment adherence was not explicitly modeled, adding uncertainty to treatment\-related effects\. Finally, the requirement for HbA1c measurements prior to the index date may introduce selection bias, as tested individuals are more likely to have suspected or established comorbidities\. As a result, the study cohort may represent a systematically sicker subset of the broader population, potentially limiting generalizability\.

## 5Conclusion

Our study demonstrates that comorbidities in T2DM can be effectively modeled as the interaction of multiple latent disease processes evolving over time\. The structural FHMM approach identified eight coherent latent components that explain the different disease trajectories, providing insights into the complex progression patterns of T2DM and its complications\.

The advantages of our approach include its ability to learn meaningful patterns of comorbidities jointly from longitudinal data, without requiring pre\-specified disease groupings\. This makes it broadly applicable, particularly in contexts where prior knowledge is unavailable or cannot be encoded\. The model captures temporal dynamics often missed by static clustering approaches and delivers more interpretable results compared to deep generative models\. The factorial structure allows for the representation of comorbidities as a superposition of multiple concurrent disease processes, offering a better understanding of disease progression than single\-process HMMs\.

Looking forward, several directions could improve the clinical utility of this approach\. While our model effectively characterizes T2DM progression patterns, real\-time clinical implementation would require incorporation of additional comorbidities beyond the core set analyzed\. Although prescribed medications were included in our initial analysis, they did not significantly improve model performance, likely due to challenges in capturing complex treatment patterns from EHR data\. Future work could explore more sophisticated feature engineering of medication data\. Nevertheless the interpretability demonstrated by our approach suggests strong potential for integration into real\-world clinical workflows, where it could provide valuable insights into disease progression patterns and inform personalized care strategies for T2DM patients\.

## 6Funding

This research was funded by a Swiss Data Science Centre Collaboration Grant \(C19\-09\)\.

## References

- \[1\]D\. Tomic, J\. E\. Shaw, and D\. J\. Magliano\(2022\)The burden and risks of emerging complications of diabetes mellitus\.Nature Reviews Endocrinology18\(9\),pp\. 525–539\.External Links:ISSN 1759\-5037,[Document](https://dx.doi.org/10.1038/s41574-022-00690-7),[Link](https://doi.org/10.1038/s41574-022-00690-7)Cited by:[§1](https://arxiv.org/html/2608.24328#S1.p1.1)\.
- \[2\]Z\. Ghahramani and M\. Jordan\(1995\)Factorial hidden markov models\.Advances in Neural Information Processing Systems8\.Cited by:[Appendix B](https://arxiv.org/html/2608.24328#A2.p1.1),[Appendix B](https://arxiv.org/html/2608.24328#A2.p3.3),[Appendix B](https://arxiv.org/html/2608.24328#A2.p4.1),[§1](https://arxiv.org/html/2608.24328#S1.p3.1),[§2\.3](https://arxiv.org/html/2608.24328#S2.SS3.p1.1),[§2\.4](https://arxiv.org/html/2608.24328#S2.SS4.p1.1),[§2\.4](https://arxiv.org/html/2608.24328#S2.SS4.p1.2),[§2\.5](https://arxiv.org/html/2608.24328#S2.SS5.p2.1)\.
- \[3\]A\. M\. Alaa, S\. Hu, and M\. Schaar\(2017\)Learning from clinical judgments: semi\-markov\-modulated marked hawkes processes for risk prognosis\.InInternational Conference on Machine Learning,pp\. 60–69\.Cited by:[§1](https://arxiv.org/html/2608.24328#S1.p5.1),[§4](https://arxiv.org/html/2608.24328#S4.p4.1)\.
- \[4\]X\. Wang, D\. Sontag, and F\. Wang\(2014\)Unsupervised learning of disease progression models\.InProceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining,pp\. 85–94\.Cited by:[§1](https://arxiv.org/html/2608.24328#S1.p5.1),[§4](https://arxiv.org/html/2608.24328#S4.p4.1)\.
- \[5\]O\. Ben\-Assuli, T\. Heart, J\. R\. Vest, R\. Ramon\-Gonen, N\. Shlomo, and R\. Klempfner\(2021\)Profiling readmissions using hidden markov model\-the case of congestive heart failure\.Information Systems Management38\(3\),pp\. 237–249\.Cited by:[§1](https://arxiv.org/html/2608.24328#S1.p5.1)\.
- \[6\]Y\. Liu, S\. Li, F\. Li, L\. Song, and J\. M\. Rehg\(2015\)Efficient learning of continuous\-time hidden markov models for disease progression\.Advances in neural information processing systems28\.Cited by:[§1](https://arxiv.org/html/2608.24328#S1.p5.1)\.
- \[7\]T\. Ceritli, A\. P\. Creagh, and D\. A\. Clifton\(2022\)Mixture of input\-output hidden markov models for heterogeneous disease progression modeling\.InWorkshop on Healthcare AI and COVID\-19,pp\. 41–53\.Cited by:[§1](https://arxiv.org/html/2608.24328#S1.p5.1)\.
- \[8\]F\. J\. Ruiz, I\. Valera, C\. Blanco, and F\. Perez\-Cruz\(2014\)Bayesian nonparametric comorbidity analysis of psychiatric disorders\.The Journal of Machine Learning Research15\(1\),pp\. 1215–1247\.Cited by:[§1](https://arxiv.org/html/2608.24328#S1.p5.1)\.
- \[9\]C\. Trottet, M\. Schürch, A\. Mollaysa, A\. Allam, and M\. Krauthammer\(2023\)Generative time series models with interpretable latent processes for complex disease trajectories\.InDeep Generative Models for Health Workshop NeurIPS 2023,External Links:[Link](https://openreview.net/forum?id=tiqs7trqcC)Cited by:[§1](https://arxiv.org/html/2608.24328#S1.p6.1),[§4](https://arxiv.org/html/2608.24328#S4.p4.1)\.
- \[10\]Z\. Oflaz, C\. Yozgatligil, and A\. S\. Selcuk\-Kestel\(2023\)Modeling comorbidity of chronic diseases using coupled hidden markov model with bivariate discrete copula\.Statistical Methods in Medical Research,pp\. 09622802231155100\.Cited by:[§1](https://arxiv.org/html/2608.24328#S1.p7.1)\.
- \[11\]L\. Fregoso\-Aparicio, J\. Noguez, L\. Montesinos, and J\. A\. García\-García\(2021\)Machine learning and deep learning predictive models for type 2 diabetes: a systematic review\.Diabetology & Metabolic Syndrome13\(1\),pp\. 1–22\.Cited by:[§1](https://arxiv.org/html/2608.24328#S1.p7.1)\.
- \[12\]A\. Martinez\-De la Torre, F\. Perez\-Cruz, S\. Weiler, and A\. M\. Burden\(2022\)Comorbidity clusters associated with newly treated type 2 diabetes mellitus: a bayesian nonparametric analysis\.Scientific Reports12\(1\),pp\. 20653\.Cited by:[§1](https://arxiv.org/html/2608.24328#S1.p7.1),[§4](https://arxiv.org/html/2608.24328#S4.p4.1)\.
- \[13\]R\. R\. Kalyani, S\. H\. Golden, and W\. T\. Cefalu\(2017\)Diabetes and aging: unique considerations and goals of care\.Diabetes Care40\(4\),pp\. 440–443\.External Links:[Document](https://dx.doi.org/10.2337/dci17-0005),[Link](https://doi.org/10.2337/dci17-0005)Cited by:[§2\.2\.1](https://arxiv.org/html/2608.24328#S2.SS2.SSS1.Px2.p1.1)\.
- \[14\]G\. T\. Russo, A\. Giandalia, E\. L\. Romeo, M\. Nunziata, M\. Muscianisi, M\. C\. Ruffo, A\. Catalano, and D\. Cucinotta\(2016\)Fracture risk in type 2 diabetes: current perspectives and gender differences\.International Journal of Endocrinology2016,pp\. 1615735\.External Links:[Document](https://dx.doi.org/10.1155/2016/1615735),[Link](https://doi.org/10.1155/2016/1615735)Cited by:[§2\.2\.1](https://arxiv.org/html/2608.24328#S2.SS2.SSS1.Px2.p1.1)\.
- \[15\]R\. Schweiger, Y\. Erlich, and S\. Carmi\(2018\)FactorialHMM: fast and exact inference in factorial hidden Markov models\.Bioinformatics35\(12\),pp\. 2162–2164\.External Links:ISSN 1367\-4803,[Document](https://dx.doi.org/10.1093/bioinformatics/bty944),[Link](https://doi.org/10.1093/bioinformatics/bty944)Cited by:[§2\.3](https://arxiv.org/html/2608.24328#S2.SS3.p4.2),[§2\.4](https://arxiv.org/html/2608.24328#S2.SS4.p2.1)\.
- \[16\]J\. Chiquet, Y\. Grandvalet, and C\. Charbonnier\(2012\)Sparsity with sign\-coherent groups of variables via the cooperative\-lasso\.The Annals of Applied Statistics6\(2\),pp\. 795–830\.External Links:ISSN 19326157, 19417330,[Link](http://www.jstor.org/stable/41713474)Cited by:[§2\.3\.1](https://arxiv.org/html/2608.24328#S2.SS3.SSS1.p4.1)\.
- \[17\]M\. Yuan and Y\. Lin\(2006\)Model selection and estimation in regression with grouped variables\.Journal of the Royal Statistical Society: Series B \(Statistical Methodology\)68\(1\),pp\. 49–67\.Cited by:[§2\.3\.1](https://arxiv.org/html/2608.24328#S2.SS3.SSS1.p4.4)\.
- \[18\]M\. J\. Wainwright and M\. I\. Jordan\(2008\)Graphical Models, Exponential Families, and Variational Inference\.Foundations and Trends in Machine Learning1\(1\-2\),pp\. 1–305\.External Links:ISSN 1935\-8237,[Link](https://doi.org/10.1561/2200000001),[Document](https://dx.doi.org/10.1561/2200000001)Cited by:[Appendix B](https://arxiv.org/html/2608.24328#A2.p1.1),[§2\.4](https://arxiv.org/html/2608.24328#S2.SS4.p1.2)\.
- \[19\]T\. S\. Jaakkola and M\. I\. Jordan\(1997\)A variational approach to Bayesian logistic regression models and their extensions\.InProceedings of the Sixth International Workshop on Artificial Intelligence and Statistics,D\. Madigan and P\. Smyth \(Eds\.\),Proceedings of Machine Learning Research, Vol\.R1,pp\. 283–294\.Note:Reissued by PMLR on 30 March 2021\.External Links:[Link](https://proceedings.mlr.press/r1/jaakkola97a.html)Cited by:[§2\.4](https://arxiv.org/html/2608.24328#S2.SS4.p2.1)\.
- \[20\]D\. L\. Davies and D\. W\. Bouldin\(1979\)A cluster separation measure\.IEEE Transactions on Pattern Analysis and Machine IntelligencePAMI\-1\(2\),pp\. 224–227\.External Links:[Document](https://dx.doi.org/10.1109/TPAMI.1979.4766909)Cited by:[Appendix D](https://arxiv.org/html/2608.24328#A4.p1.1),[§2\.5](https://arxiv.org/html/2608.24328#S2.SS5.p2.2)\.
- \[21\]R\. Tibshirani and G\. Walther\(2005\)Cluster validation by prediction strength\.Journal of Computational and Graphical Statistics14\(3\),pp\. 511–528\.External Links:[Document](https://dx.doi.org/10.1198/106186005X59243),[Link](https://doi.org/10.1198/106186005X59243)Cited by:[Appendix D](https://arxiv.org/html/2608.24328#A4.p1.1),[§2\.5](https://arxiv.org/html/2608.24328#S2.SS5.p2.2)\.
- \[22\]H\. Akaike\(1974\)A new look at the statistical model identification\.IEEE Transactions on Automatic Control19\(6\),pp\. 716–723\.External Links:[Document](https://dx.doi.org/10.1109/TAC.1974.1100705)Cited by:[§2\.6](https://arxiv.org/html/2608.24328#S2.SS6.p2.1)\.
- \[23\]G\. Schwarz\(1978\)Estimating the dimension of a model\.The Annals of Statistics6\(2\),pp\. 461–464\.External Links:ISSN 00905364, 21688966,[Link](http://www.jstor.org/stable/2958889)Cited by:[§2\.6](https://arxiv.org/html/2608.24328#S2.SS6.p2.1)\.
- \[24\]Z\. M\. Younossi, A\. B\. Koenig, D\. Abdelatif, Y\. Fazel, L\. Henry, and M\. Wymer\(2016\)Global epidemiology of nonalcoholic fatty liver disease—meta‐analytic assessment of prevalence, incidence, and outcomes\.Hepatology64\(1\),pp\. 73–84\.External Links:[Document](https://dx.doi.org/10.1002/hep.28431)Cited by:[§3\.3](https://arxiv.org/html/2608.24328#S3.SS3.SSS0.Px3.p5.1)\.
- \[25\]O\. Aalen\(1978\)Nonparametric estimation of partial transition probabilities in multiple decrement models\.The Annals of Statistics6\(3\),pp\. 534–545\.External Links:ISSN 00905364, 21688966,[Link](http://www.jstor.org/stable/2958557)Cited by:[Figure 16](https://arxiv.org/html/2608.24328#A5.F16),[Figure 16](https://arxiv.org/html/2608.24328#A5.F16.5),[Figure 20](https://arxiv.org/html/2608.24328#A6.F20),[Figure 20](https://arxiv.org/html/2608.24328#A6.F20.5),[Figure 24](https://arxiv.org/html/2608.24328#A7.F24),[Figure 24](https://arxiv.org/html/2608.24328#A7.F24.5),[Figure 9](https://arxiv.org/html/2608.24328#S3.F9),[Figure 9](https://arxiv.org/html/2608.24328#S3.F9.5)\.
- \[26\]C\. Ronco, P\. McCullough, S\. D\. Anker, I\. Anand, N\. Aspromonte, S\. M\. Bagshaw, R\. Bellomo, T\. Berl, I\. Bobek, D\. Cruz,et al\.\(2008\)Cardiorenal syndrome\.Journal of the American College of Cardiology52\(19\),pp\. 1527–1539\.External Links:[Document](https://dx.doi.org/10.1016/j.jacc.2008.07.051)Cited by:[§4](https://arxiv.org/html/2608.24328#S4.p1.1)\.
- \[27\]N\. Cheung, P\. Mitchell, and T\. Y\. Wong\(2010\)Diabetic retinopathy\.Lancet376\(9735\),pp\. 124–136\.Note:Epub 2010 Jun 26External Links:[Document](https://dx.doi.org/10.1016/S0140-6736%2809%2962124-3)Cited by:[§4](https://arxiv.org/html/2608.24328#S4.p1.1)\.
- \[28\]R\. Lee, T\. Y\. Wong, and C\. Sabanayagam\(2015\)Epidemiology of diabetic retinopathy, diabetic macular edema and related vision loss\.Eye and Vision \(London, England\)2,pp\. 17\.External Links:[Document](https://dx.doi.org/10.1186/s40662-015-0026-2),[Link](https://doi.org/10.1186/s40662-015-0026-2)Cited by:[§4](https://arxiv.org/html/2608.24328#S4.p1.1)\.
- \[29\]P\. Vestergaard\(2007\)Discrepancies in bone mineral density and fracture risk in patients with type 1 and type 2 diabetes—a meta\-analysis\.Osteoporosis International18\(4\),pp\. 427–444\.External Links:[Document](https://dx.doi.org/10.1007/s00198-006-0253-4)Cited by:[§4](https://arxiv.org/html/2608.24328#S4.p1.1)\.
- \[30\]N\. Napoli, M\. Chandran, D\. D\. Pierroz, B\. Abrahamsen, A\. V\. Schwartz, and S\. L\. Ferrari\(2017\)Mechanisms of diabetes mellitus\-induced bone fragility\.Nature Reviews Endocrinology13\(4\),pp\. 208–219\.External Links:[Document](https://dx.doi.org/10.1038/nrendo.2016.153)Cited by:[§4](https://arxiv.org/html/2608.24328#S4.p1.1)\.

## Appendix AData preparation

Table 5:Predefined thresholds used to remove outliers\. Values are removed if they are not included in the lower and upper limits\.Table 6:Number of patient diagnosed with the comorbidity for each data split\. Numbers are given in percentage\.
## Appendix BVariational approximation

This section provides details on the fitting algorithm used to train the FHMM\. The fitting algorithm mainly follows\[[2](https://arxiv.org/html/2608.24328#bib.bib1)\], where the authors reduce the inference complexity toO⁡\(K2​M​T\)O\(K^\{2\}MT\)by exploiting a structured mean field variational approximation\[[18](https://arxiv.org/html/2608.24328#bib.bib15)\]for the posterior distributionqϕq\_\{\\phi\}used in the E step of their EM procedure\. In particular, due to their specific choice, \([6](https://arxiv.org/html/2608.24328#S2.E6)\) can be written as

F\(qϕ,θ\)=𝔼qϕ\[∑t=1Tlogpθ\(yt\|H~t\(1:M\)\)−∑t=1T∑m=1M\(H~t\(m\)\)⊤loget\(m\)\]−Z\+Zqϕ,F\(q\_\{\\phi\},\\theta\)=\\mathbb\{E\}\_\{q\_\{\\phi\}\}\\left\[\\sum\_\{t=1\}^\{T\}\\log p\_\{\\theta\}\\left\(y\_\{t\}\|\\tilde\{H\}\_\{t\}^\{\(1:M\)\}\\right\)\-\\sum\_\{t=1\}^\{T\}\\sum\_\{m=1\}^\{M\}\(\\tilde\{H\}\_\{t\}^\{\(m\)\}\)^\{\\top\}\\log e\_\{t\}^\{\(m\)\}\\right\]\-Z\+Z\_\{q\_\{\\phi\}\},whereZZandZqϕZ\_\{q\_\{\\phi\}\}are respectively normalizing constants ofpθ\(y1:T\|h1:T,x1:T\)p\_\{\\theta\}\(y\_\{1:T\}\|h\_\{1:T\},x\_\{1:T\}\)andqϕ\(h1:T∣x1:T\)q\_\{\\phi\}\(h\_\{1:T\}\\mid x\_\{1:T\}\)and whereH~t\(m\)\\tilde\{H\}\_\{t\}^\{\(m\)\}is a one\-hot encoded version ofHt\(m\)H\_\{t\}^\{\(m\)\}\.

The E step requires the optimization of the parametersϕ\\phi, which can be computing by solving

∂F∂log⁡eτ\(n\)\|qϕ^,θ=0,\\left\.\\frac\{\\partial F\}\{\\partial\\log e\_\{\\tau\}^\{\(n\)\}\}\\right\|\_\{q\_\{\\hat\{\\phi\}\},\\theta\}=0,where the gradient is

∂F∂log⁡eτ\(n\)\|qϕ,θ=∑t=1T𝔼qϕ\[\{logpθ\(yt\|H~t\(1:M\)\)−∑m=1M\(H~t\(m\)\)⊤loget\(m\)\}×\{H~τ\(n\)−𝔼qϕ\[H~τ\(n\)\]\}\]\.\\left\.\\frac\{\\partial F\}\{\\partial\\log e\_\{\\tau\}^\{\(n\)\}\}\\right\|\_\{q\_\{\{\\phi\}\},\\theta\}=\\sum\_\{t=1\}^\{T\}\\mathbb\{E\}\_\{q\_\{\\phi\}\}\\left\[\\left\\\{\\log p\_\{\\theta\}\\left\(y\_\{t\}\|\\tilde\{H\}\_\{t\}^\{\(1:M\)\}\\right\)\-\\sum\_\{m=1\}^\{M\}\\left\(\\tilde\{H\}\_\{t\}^\{\(m\)\}\\right\)^\{\\top\}\\log e\_\{t\}^\{\(m\)\}\\right\\\}\\right\.\\left\.\\times\\left\\\{\{\\tilde\{H\}\_\{\\tau\}^\{\(n\)\}\}\-\\mathbb\{E\}\_\{q\_\{\\phi\}\}\\left\[\\tilde\{H\}\_\{\\tau\}^\{\(n\)\}\\right\]\\right\\\}\\right\]\.\(9\)
The M step requires the optimization of𝔼qϕ\[logpθ\(y1:T,h1:T\)\]\\mathbb\{E\}\_\{q\_\{\\phi\}\}\\left\[\\log p\_\{\\theta\}\(y\_\{1:T\},h\_\{1:T\}\)\\right\]with respect toθπ,θP,η\\theta^\{\\pi\},\\theta^\{P\},\\eta\. Each parameter can be optimized individually by solving

θ^π\\displaystyle\\hat\{\\theta\}^\{\\pi\}=arg​maxθπ∑mγ1\(m\)⊤logπθπ\(m\),\\displaystyle=\\argmax\_\{\\theta^\{\\pi\}\}\\sum\_\{m\}\{\\gamma\_\{1\}^\{\(m\)\}\}^\{\\top\}\\log\\pi\_\{\\theta^\{\\pi\}\}^\{\(m\)\},θ^P\\displaystyle\\hat\{\\theta\}^\{P\}=arg​maxθP∑m∑tlogPθP\(m\)Tr\(ξt\(m\)\),\\displaystyle=\\argmax\_\{\\theta^\{P\}\}\\sum\_\{m\}\\sum\_\{t\}\\log P\_\{\\theta^\{P\}\}^\{\(m\)\}\\Tr\(\\xi\_\{t\}^\{\(m\)\}\),η^\\displaystyle\\hat\{\\eta\}=arg​maxη∑t𝔼qϕ\[logpη\(yt∣ht\(1:M\)\)\],\\displaystyle=\\argmax\_\{\\eta\}\\sum\_\{t\}\\mathbb\{E\}\_\{q\_\{\\phi\}\}\\left\[\\log p\_\{\\eta\}\(y\_\{t\}\\mid\{h\_\{t\}^\{\(1:M\)\}\)\}\\right\],\(10\)where

γt\(m\)\\displaystyle\\gamma\_\{t\}^\{\(m\)\}=def𝔼qϕ​\[H~t\(m\)\]∈ℝK,\\displaystyle\\stackrel\{\{\\scriptstyle\\rm def\}\}\{\{=\}\}\\mathbb\{E\}\_\{q\_\{\\phi\}\}\\left\[\\tilde\{H\}\_\{t\}^\{\(m\)\}\\right\]\\in\\mathbb\{R\}^\{K\},ξt\(m\)\\displaystyle\\xi\_\{t\}^\{\(m\)\}=def𝔼qϕ​\[H~t\(m\)​H~t−1\(m\)⊤\]∈ℝK×K,\\displaystyle\\stackrel\{\{\\scriptstyle\\rm def\}\}\{\{=\}\}\\mathbb\{E\}\_\{q\_\{\\phi\}\}\\left\[\\tilde\{H\}\_\{t\}^\{\(m\)\}\{\\tilde\{H\}\_\{t\-1\}^\{\(m\)\}\\ \}^\{\\top\}\\right\]\\in\\mathbb\{R\}^\{K\\times K\},are computed using the Baum–Welch algorithm for HMMs \(see\[[2](https://arxiv.org/html/2608.24328#bib.bib1)\]\), whose complexity isO⁡\(K2​T\)O\(K^\{2\}T\)\.

Equations \([9](https://arxiv.org/html/2608.24328#A2.E9)\) and \([10](https://arxiv.org/html/2608.24328#A2.E10)\) involves expectations that are expensive to compute\. In particular, it requires the computation ofq\(Ht1\(m\),Ht2\(m\)∣x1:T\)q\(H\_\{t\_\{1\}\}^\{\(m\)\},H\_\{t\_\{2\}\}^\{\(m\)\}\\mid x\_\{1:T\}\)for anyt1,t2t\_\{1\},t\_\{2\}, which has a memory and computational complexity ofO⁡\(M​T2​K2\)O\(MT^\{2\}K^\{2\}\)\. Inspired by the Gaussian case, for which analytical expressions exists \(see\[[2](https://arxiv.org/html/2608.24328#bib.bib1)\]\), a second order approximation of the emission log\-likelihoodlogpη\(yt∣ht\(1:M\)\)\\log p\_\{\\eta\}\(y\_\{t\}\\mid\{h\_\{t\}^\{\(1:M\)\}\}\)is derived using the inequality \([8](https://arxiv.org/html/2608.24328#S2.E8)\)\. In that case, the gradient is reduced to

∂F∂log⁡eτ\(n\)\|q,θ≈∑tCov\(H~τ\(n\),H~t\(n\)\)×\{ytη\(n\)−η\(n\)\(12\+λ\(z\)η\(n\)η⊤γt\)−λ⁡\(z\)2\(Diag\(η\(n\)η\(n\)⊤\)−2η\(n\)η\(n\)⊤γ\(n\)\)−loget\(n\)\}\.\\left\.\\frac\{\\partial F\}\{\\partial\\log e\_\{\\tau\}^\{\(n\)\}\}\\right\|\_\{q,\\theta\}\\approx\\sum\_\{t\}\\cov\\left\(\\tilde\{H\}\_\{\\tau\}^\{\(n\)\},\\tilde\{H\}\_\{t\}^\{\(n\)\}\\right\)\\times\\left\\\{y\_\{t\}\\eta^\{\(n\)\}\-\\eta^\{\(n\)\}\\left\(\\frac\{1\}\{2\}\+\\lambda\(z\)\\eta^\{\(n\)\}\\eta^\{\\top\}\\gamma\_\{t\}\\right\)\\right\.\\\\ \\left\.\-\\frac\{\\lambda\(z\)\}\{2\}\\left\(\\diag\(\\eta^\{\(n\)\}\{\\eta^\{\(n\)\}\}^\{\\top\}\)\-2\\eta^\{\(n\)\}\{\\eta^\{\(n\)\}\}^\{\\top\}\\gamma^\{\(n\)\}\\right\)\-\\log e\_\{t\}^\{\(n\)\}\\right\\\}\.\(11\)Finding the zeros of \([11](https://arxiv.org/html/2608.24328#A2.E11)\) \(E step\) can then be easily done without computing the expensive termsCov\(H~τ\(n\),H~t\(n\)\)\\cov\\left\(\\tilde\{H\}\_\{\\tau\}^\{\(n\)\},\\tilde\{H\}\_\{t\}^\{\(n\)\}\\right\)\. On the other hand, the M step for binary outputs becomes

maxη𝔼qϕ\[logpη\(yt∣ht\(1:M\)\)\]=maxη𝔼qϕ\[ytH~t⊤η−log\(1\+exp\(H~t⊤η\)\)\]≥maxη⁡\(yt−12\)​γt⊤​η−log⁡\{1\+exp⁡\(−z\)\}−z/2−λ⁡\(z\)2​\(η⊤​𝔼​\[H~t​H~t⊤\]​η−z2\),\\max\_\{\\eta\}\\mathbb\{E\}\_\{q\_\{\\phi\}\}\\left\[\\log p\_\{\\eta\}\(y\_\{t\}\\mid h\_\{t\}^\{\(1:M\)\}\)\\right\]=\\max\_\{\\eta\}\\mathbb\{E\}\_\{q\_\{\\phi\}\}\\left\[y\_\{t\}\\tilde\{H\}\_\{t\}^\{\\top\}\\eta\-\\log\(1\+\\exp\(\\tilde\{H\}\_\{t\}^\{\\top\}\\eta\)\)\\right\]\\geq\\\\ \\max\_\{\\eta\}\\left\(y\_\{t\}\-\\frac\{1\}\{2\}\\right\)\\gamma\_\{t\}^\{\\top\}\\eta\-\\log\\left\\\{1\+\\exp\(\-z\)\\right\\\}\-z/2\-\\frac\{\\lambda\(z\)\}\{2\}\\left\(\\eta^\{\\top\}\\mathbb\{E\}\\left\[\\tilde\{H\}\_\{t\}\{\\tilde\{H\}\_\{t\}\}^\{\\top\}\\right\]\\eta\-z^\{2\}\\right\),\(12\)whereH~t=\(H~t\(1\),…,H~t\(M\)\)\\tilde\{H\}\_\{t\}=\\left\(\\tilde\{H\}\_\{t\}^\{\(1\)\},\\dots,\\tilde\{H\}\_\{t\}^\{\(M\)\}\\right\)andγt=\(γt\(1\),…,γt\(M\)\)\\gamma\_\{t\}=\\left\(\\gamma\_\{t\}^\{\(1\)\},\\dots,\\gamma\_\{t\}^\{\(M\)\}\\right\)\. Note that Equation[12](https://arxiv.org/html/2608.24328#A2.E12)is a lower bound of the ELBO, which can be optimized with resepect tozz, which admits a closed form solution

zt=η⊤​𝔼​\[H~t​H~t⊤\]​η\.z\_\{t\}=\\sqrt\{\\eta^\{\\top\}\\mathbb\{E\}\\left\[\\tilde\{H\}\_\{t\}\{\\tilde\{H\}\_\{t\}\}^\{\\top\}\\right\]\\eta\}\.

## Appendix CInput medications

When drugs are used as input, the entries of the matrixPθPP\_\{\\theta^\{P\}\}are given

\[PθP\]i​j∝exp\(αi​jdrug\+fl\(Xt,i\)𝟙\{i\>j\}\+fr\(Xt,i\)𝟙\{j\>i\}\),\[P\_\{\\theta^\{P\}\}\]\_\{ij\}\\propto\\exp\(\\alpha\_\{ij\}^\{\\text\{drug\}\}\+f\_\{l\}\(X\_\{t\},i\)\\mathbbm\{1\}\_\{\\\{i\>j\\\}\}\+f\_\{r\}\(X\_\{t\},i\)\\mathbbm\{1\}\_\{\\\{j\>i\\\}\}\),wherefl​\(Xt,i\)=Wl,idrug⋅Xt≥0f\_\{l\}\(X\_\{t\},i\)=W\_\{l,i\}^\{\\text\{drug\}\}\\cdot X\_\{t\}\\geq 0andfr​\(Xt,i\)=Wr,idrug⋅Xtf\_\{r\}\(X\_\{t\},i\)=W\_\{r,i\}^\{\\text\{drug\}\}\\cdot X\_\{t\}\. This parametrization captures the fact that drugs can either add probability mass on lower state values \(Wl,idrug\>0W\_\{l,i\}^\{\\text\{drug\}\}\>0\) or add probability mass on higher state values \(Wr,idrug\>0W\_\{r,i\}^\{\\text\{drug\}\}\>0\)\. It is important to note that adding probability mass to lower state values will automatically reduce the probability mass on higher state values, and vice versa\. In addition,Wl,idrug,Wr,idrugW\_\{l,i\}^\{\\text\{drug\}\},W\_\{r,i\}^\{\\text\{drug\}\}, follow a priori a Gaussian distribution centered at 0, which is equivalent to add the penalizations‖Wl,idrug‖22\\\|W\_\{l,i\}^\{\\text\{drug\}\}\\\|\_\{2\}^\{2\}to the objective function\. The underlying assumptions is that the drug effect is a priori small \(i\.e\.,Wl,idrugW\_\{l,i\}^\{\\text\{drug\}\}is likely to be small\)\.

Table[7](https://arxiv.org/html/2608.24328#A3.T7)shows the performance of the model when drugs are used as input\. The use of drugs in the model, is marginal and it does not improve the BIC score\. The AIC is improved in few cases\. Figure[11](https://arxiv.org/html/2608.24328#A3.F11)show how the input variables affect the evolution of the latent state for the modelFHMM​\(8,4\)2\\text\{FHMM\}\(8,4\)\_\{2\}fitted on the male population younger than 65\. Age generally has a negative effect on chain transitions, increasing the probability of moving to a worsening state\. Smoking also accelerates progression along chain 4, which is associated with a higher risk of COPD\. The use of calcium channel blockers is linked to chain 5, which models hypertension\. Here we observe an apparent worsening effect, although this may simply reflect that the model associates higher states with greater use of calcium channel blockers\. By contrast, the addition of a new antidiabetic drug typically produces a positive effect on the patient’s state\.

Table 7:Comparison between models with drugs \(XtX\_\{t\}\) and without drugs\. Metrics are evaluated on the validation set\.![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male_m65_8_4_input.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male_m65_8_4_bin.png)

Figure 11:Top:Average effect of drugs on the transition matrix, estimated as1/K∑i=1K\(Wr,idrug−Wl,idrug\)1/K\\sum\_\{i=1\}^\{K\}\\left\(W\_\{r,i\}^\{\\text\{drug\}\}\-W\_\{l,i\}^\{\\text\{drug\}\}\\right\)\. Red \(resp\. blue\) means an increase of the probability to transition to a higher \(resp\. lower\) state value\.Bottom:Emission parametersη\\etamodeling the risk of developing a comorbidity for a given state\. Larger values indicate a higher probability of having the comorbidity\. The model includes 8 chains with 4 states per chain and is fitted on*Men<65<65*\.
## Appendix DAdditional Figures

The number of clusters is determined using the Davies\-Bouldin index\[[20](https://arxiv.org/html/2608.24328#bib.bib21)\]and the prediction strength introduced by\[[21](https://arxiv.org/html/2608.24328#bib.bib22)\]\. The prediction strength is estimated using a 3\-fold cross\-validation setup, where two\-thirds of the data is used for training and one\-third for testing\.

![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female_m65/selection.png)Figure 12:Cluster selection using Davies\-Bouldin index and the prediction strength\. The optimal number of cluster is indicated by the black dashed line\. This is fitted on the population ofWomen<65<65\.![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female_m65/flow_chart_cluster_0.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female_m65/flow_chart_cluster_4.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female_m65/flow_chart_cluster_11.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female_m65/flow_chart_cluster_2.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female_m65/flow_chart_cluster_10.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female_m65/flow_chart_cluster_3.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female_m65/flow_chart_cluster_6.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female_m65/flow_chart_cluster_8.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female_m65/flow_chart_cluster_9.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female_m65/flow_chart_cluster_12.png)

Figure 13:Evolution of patient clusters 0 \(healthy\), 2 \(pathway to multimorbidity\), 3 \(pulmonary\-metabolic\), 4 \(isolated hypertension\), 6 \(hepatic\), 8 \(advanced nephropathy\), 9 \(obese hypertension\), 10 \(bone\-vascular\), 11 \(isolated DR\), and 12 \(cardio\-pulmonary multimorbidity\) over time \(Women<65<65\), conditional on the initial cluster assignment \(6 months post\-NIAD\)\. The number of patients in each cluster is indicated alongside the corresponding bars\. For clarity, flows representing less than 1% of the initial cluster size, as well as transitions to no available follow\-up data, are not displayed\. Clusters 1, 5, and 7 are excluded due to the low number of transitions observed from these clusters\.
## Appendix EMale<65<65

Table 8:Patient characteristics in each cluster for theMen<65<65cohort with average values and standard deviation where applicable\. Comorbidities are presented as absolute counts at NIAD initiation, along with the number of new comorbidities developed post\-NIAD\. Proportions, expressed as percentages of the cohort population, are reported in parentheses\.![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male_m65/proba_avg.png)Figure 14:Average cluster\-specific probabilities and lab values forMen<65<65\. Values have been centered by subtracting the average at the time of NIAD initiation to improve readability\. Red bars indicate increases, and blue bars indicate decreases relative to the average patient at NIAD\. Comorbidity risks are shown as changes in probability \(left axis\), and lab values as percentage changes relative to the reference value \(right axis\)\. Annotated numbers indicate absolute changes in the original units \(see Table[1](https://arxiv.org/html/2608.24328#S2.T1)\)\. Error bars represent the standard error of the mean\.![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male_m65/flow_chart_cluster_0.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male_m65/flow_chart_cluster_1.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male_m65/flow_chart_cluster_2.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male_m65/flow_chart_cluster_3.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male_m65/flow_chart_cluster_4.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male_m65/flow_chart_cluster_5.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male_m65/flow_chart_cluster_7.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male_m65/flow_chart_cluster_8.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male_m65/flow_chart_cluster_9.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male_m65/flow_chart_cluster_10.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male_m65/flow_chart_cluster_11.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male_m65/flow_chart_cluster_12.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male_m65/flow_chart_cluster_13.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male_m65/flow_chart_cluster_14.png)

Figure 15:Evolution of patient clusters over time \(Men<65<65\), conditional on the initial cluster assignment \(6 months post\-NIAD\)\. The number of patients in each cluster is indicated alongside the corresponding bars\. For clarity, only the most frequent transitions are shown \(see Figure[16](https://arxiv.org/html/2608.24328#A5.F16)for more details\)\. Cluster 6 is excluded due to the low number of transitions observed from this cluster\.![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male_m65/cumulative_incidence.png)Figure 16:Cumulative Incidence Function \(CIF\) for the time to first cluster transition inMen<65<65, stratified by initial cluster assignment at 6 months post\-NIAD\. The estimates account for censoring due to loss to follow\-up and study termination, computed using the Aalen nonparametric estimator\[[25](https://arxiv.org/html/2608.24328#bib.bib6)\]\.![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male_m65/graph_avg.png)Figure 17:Transition graph\. Red arrows \(resp\. blue\) means an increase \(resp\. decrease\) of the comorbidity risk \(Men<65<65\)\. Arrow thickness is proportional to the number of observed transitions\. Transitions with less than 1% probability are not shown\.
## Appendix FMale65\+65\+

Table 9:Patient characteristics in each cluster for theMen65\+65\+cohort with average values and standard deviation where applicable\. Comorbidities are presented as absolute counts at NIAD initiation, along with the number of new comorbidities developed post\-NIAD\. Proportions, expressed as percentages of the cohort population, are reported in parentheses\.![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male65+/proba_avg.png)Figure 18:Average cluster\-specific probabilities and lab values forMen65\+65\+\. Values have been centered by subtracting the average at the time of NIAD initiation to improve readability\. Red bars indicate increases, and blue bars indicate decreases relative to the average patient at NIAD\. Comorbidity risks are shown as changes in probability \(left axis\), and lab values as percentage changes relative to the reference value \(right axis\)\. Annotated numbers indicate absolute changes in the original units \(see Table[1](https://arxiv.org/html/2608.24328#S2.T1)\)\. Error bars represent the standard error of the mean\.![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male65+/flow_chart_cluster_0.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male65+/flow_chart_cluster_1.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male65+/flow_chart_cluster_4.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male65+/flow_chart_cluster_5.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male65+/flow_chart_cluster_7.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male65+/flow_chart_cluster_9.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male65+/flow_chart_cluster_12.png)

Figure 19:Evolution of patient clusters over time \(Men65\+65\+\), conditional on the initial cluster assignment \(6 months post\-NIAD\)\. The number of patients in each cluster is indicated alongside the corresponding bars\. For clarity, only the most frequent transitions are shown \(see Figure[20](https://arxiv.org/html/2608.24328#A6.F20)for more details\)\. Clusters 2,3,6,8,10,11 are excluded due to the low number of transitions observed from these clusters\.![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male65+/cumulative_incidence.png)Figure 20:Cumulative Incidence Function \(CIF\) for the time to first cluster transition inMen65\+65\+, stratified by initial cluster assignment at 6 months post\-NIAD\. The estimates account for censoring due to loss to follow\-up and study termination, computed using the Aalen nonparametric estimator\[[25](https://arxiv.org/html/2608.24328#bib.bib6)\]\.![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/male65+/graph_avg.png)Figure 21:Transition graph\. Red arrows \(resp\. blue\) means an increase \(resp\. decrease\) of the comorbidity risk \(Men65\+65\+\)\. Arrow thickness is proportional to the number of observed transitions\. Transitions with less than 1% probability are not shown\.
## Appendix GFemale65\+65\+

Table 10:Patient characteristics in each cluster for theWomen65\+65\+cohort with average values and standard deviation where applicable\. Comorbidities are presented as absolute counts at NIAD initiation, along with the number of new comorbidities developed post\-NIAD\. Proportions, expressed as percentages of the cohort population, are reported in parentheses\.![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female65+/proba_avg.png)Figure 22:Average cluster\-specific probabilities and lab values forWomen65\+65\+\. Values have been centered by subtracting the average at the time of NIAD initiation to improve readability\. Red bars indicate increases, and blue bars indicate decreases relative to the average patient at NIAD\. Comorbidity risks are shown as changes in probability \(left axis\), and lab values as percentage changes relative to the reference value \(right axis\)\. Annotated numbers indicate absolute changes in the original units \(see Table[1](https://arxiv.org/html/2608.24328#S2.T1)\)\. Error bars represent the standard error of the mean\.![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female65+/flow_chart_cluster_0.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female65+/flow_chart_cluster_2.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female65+/flow_chart_cluster_3.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female65+/flow_chart_cluster_4.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female65+/flow_chart_cluster_5.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female65+/flow_chart_cluster_6.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female65+/flow_chart_cluster_7.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female65+/flow_chart_cluster_9.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female65+/flow_chart_cluster_10.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female65+/flow_chart_cluster_11.png)
![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female65+/flow_chart_cluster_12.png)

Figure 23:Evolution of patient clusters over time \(Women65\+65\+\), conditional on the initial cluster assignment \(6 months post\-NIAD\)\. The number of patients in each cluster is indicated alongside the corresponding bars\. For clarity, only the most frequent transitions are shown \(see Figure[24](https://arxiv.org/html/2608.24328#A7.F24)for more details\)\. Clusters 1,8 are excluded due to the low number of transitions observed from these clusters\.![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female65+/cumulative_incidence.png)Figure 24:Cumulative Incidence Function \(CIF\) for the time to first cluster transition inWomen65\+65\+, stratified by initial cluster assignment at 6 months post\-NIAD\. The estimates account for censoring due to loss to follow\-up and study termination, computed using the Aalen nonparametric estimator\[[25](https://arxiv.org/html/2608.24328#bib.bib6)\]\.![Refer to caption](https://arxiv.org/html/2608.24328v1/figure/female65+/graph_avg.png)Figure 25:Transition graph\. Red arrows \(resp\. blue\) means an increase \(resp\. decrease\) of the comorbidity risk \(Women65\+65\+\)\. Arrow thickness is proportional to the number of observed transitions\. Transitions with less than 1% probability are not shown\.