Cirrhosis is a chronic liver disease that features fibrosis and regenerative nodules [1]. It can be classified clinically as compensated or decompensated [2]. The former entails an early disease stage characterized by a significantly scarred liver that still functions sufficiently to perform its essential responsibilities, without overt symptoms [3]. The latter denotes an advanced stage and is manifested by substantially decreased functional capacity that often leads to serious prognoses [4]. Ascites, varices, hepatic encephalopathy (HE), and hepatorenal syndrome (HRS) are important symptoms that mark the onset of liver decompensation [5]. When patients deteriorate from the compensated to the decompensated stage, they face a much higher mortality rate and need substantial clinical resources. Volk et al [6] stressed the criticality of identifying patients at higher risk of decompensation, so that they can be monitored closely and receive timely interventions [7].

However, accurately estimating the risk that patients with cirrhosis likely deteriorate from the compensated to the decompensated stage is difficult, largely due to the inherently multifactorial and heterogeneous nature of liver cirrhosis progression. Effective predictions are further complicated by interrelated variables that pertain to liver function measures, lifestyle factors, comorbid conditions, and temporal disease dynamics [8]. Physicians tend to depend on invasive hepatic venous pressure gradient measures, which, however, are not conveniently available in many health care organizations [9,10]. Noninvasive tests and serum-based biomarkers are also accessible to physicians. For example, transient elastography (eg, FibroScan; Echosens) and other modalities (eg, acoustic radiation using impulse imaging and magnetic resonance elastography) can assess liver stiffness to estimate decompensation risk indirectly [9,11,12]. Most physicians often use score-based methods, such as the Child-Pugh score, to estimate a patient’s likelihood of deterioration. These methods are usually restricted by their exclusion of essential etiologies and can benefit from being complemented by other methods or additional factors [13,14]. For example, other clinical measures that consider multiple recurrent events of different types in the patient history, such as disease burden score, can be incorporated as an outcome variable in regression models [15].

Data-driven predictive analytics is particularly promising and can complement score-based methods [16]. Toward that end, machine learning–based methods integrate structured, longitudinal clinical variables (eg, laboratory results and elastography measures) to capture complex relationships and interplays among different conditions and medications, which are not attainable by fixed-point, score-based methods [17]. Similarly, deep learning–based methods are suited for estimating patient outcomes and risks because they can model complex relationships and temporal dependencies for individual-level predictions, using clinical, biomedical, and time series patient data available in electronic health records (EHRs) [18]. A recent review of deep learning shows that these methods have become increasingly applied in clinical settings because of their inclusion of temporal intervals, multiple feature types, and larger sample sizes for better predictive performance [19]. However, deep learning–based methods also have their own challenges; for example, predictive models built from EHR data require bias mitigation and fairness evaluation to ensure equitable performance across patient subgroups. These challenges underscore the need for evaluation beyond prediction accuracy [20]. Overall, predictive analytics methods can leverage detailed patient data and harness algorithmic power to complement score-based methods.

Different data-driven predictive methods can be applied to estimate decompensation among patients with cirrhosis [21], but many appear limited in prediction effectiveness [22]. For example, most existing methods target a single binary outcome, such as compensated versus decompensated or presence versus absence of a specific symptom, and often ignore the multifaceted nature of disease progression that involves multiple, distinct decompensation symptoms [23]. This implies the clinical value of multitask learning for modeling multiple, interrelated patient outcomes simultaneously to better support downstream tasks and predictions, as exemplified by methods incorporating attention mechanisms [24]. In addition, many methods overlook demographic similarity, which can inform decompensation symptom estimates by leveraging individual risk profiles across different patient groups [25,26]. Nor do prevalent methods consider symptom interplays; for example, the onset of ascites might be associated with a patient’s subsequent risk of HE or HRS. Modeling such interplays could provide a means to exploit valuable intersymptom relationships for more accurate patient risk estimates and richer clinical insights [27]. In a related sense, modeling symptom interplays can address the “black box” nature of deep learning and provide post hoc interpretability. To illustrate, patient-level comorbidity graphs allow for visualizations of the relationships among diseases to facilitate clinicians’ comprehension, reasoning, and adoption.

The main prediction task was symptom-level decompensation risk. We estimated 4 distinct decompensation symptoms among patients within a 3-year horizon after their cirrhosis diagnosis: ascites, varices, HE, and HRS. A patient’s transition from the compensated to the decompensated stage is manifested by the onset of 1 or more of these symptoms. Guided by clinical knowledge about cirrhosis and its progression, we developed a graph-based multitask deep learning method that considers patient demographics, etiological characteristics, medications, and interdisease interactions. This method also accounts for the relationship of decompensation symptoms with mortality, the inclusion of which in the predictive model construction can augment the main prediction task because mortality is a critical patient outcome resulting from decompensation [28].

The proposed method includes mortality prediction as an auxiliary task during model training, which helps the method learn meaningful information of disease severity and adverse progression. Clinical occurrence of decompensation symptoms is highly associated with worsening prognosis, making mortality a useful complementary signal for training a predictive model, although it is not needed when applying the method to new patients (for decompensation symptom prediction). Incorporating this important patient outcome helps the model encode shared representations of risk across related clinical end points and can strengthen the identification of patients who are likely to develop decompensation symptoms. In this way, the auxiliary task enriches the model learning process by leveraging an essential and closely related outcome to reflect the disease progression more broadly.

Moreover, patients are stratified into cohorts based on their condition, and a comorbidity graph is constructed to encode clinically plausible interdisease relationships. Compared with existing methods, the proposed method is novel in 3 important aspects. First, it adopts a multitask framework and incorporates mortality prediction as an auxiliary task to better estimate distinct symptoms simultaneously using a graph neural network (GNN). Second, this method incorporates patient similarity to create different patient groups and constructs group-level graphs to model comorbidity interplays for predictive efficacy, in line with cirrhosis literature recognizing that similar patients tend to have more comparable outcomes than dissimilar patients [25]. Third, it integrates signals from individual- and group-level embeddings to account for both individual heterogeneity and cohort characteristics, so it can distinguish disease chronicity and further enhance predictions by considering temporal dynamics of disease. Furthermore, the proposed method offers post hoc interpretability to help clinicians understand its predictions for clinical decision-making and use.

This study underscores the advantages of graph-based modeling and multitask learning for capturing comorbidity interplays, patient demographics, and symptoms to better identify patients with cirrhosis at higher risk of decompensation in a timely manner. Evaluation results affirm the proposed method’s superior predictive utilities over many prevalent methods, thus demonstrating the clinical value and feasibility of predicting multiple decompensation symptoms simultaneously. This method adds to predictive analytics research beyond a simplistic, binary view and enables nuanced risk stratification and personalized management of individual patients with a specific decompensation trajectory. The proposed method is generalizable and can be adapted to other chronic diseases featuring similar multifactorial progression and important symptom interplays that challenge analytics methods. Supported by timely and accurate decompensation estimates at the individual level, physicians can identify higher-risk patients for proactive management, therapeutic planning, and efficient resource allocation, with cost containment implications.

Let P={P1,P2,…,Pi} be a set of patients and D={<di,1,δi,1>,<di,2,δi,2>,…,<di,D,δi,D>} be a set of diseases (comorbidities), where di,j is the binary indicator of morbidity j for patient Pi, and δi,j is the temporal difference (time interval) between the onset of that comorbidity and the cirrhosis diagnosis. M={<medi,1,dosagei,1,duri,1,Si,1>,<medi,2,dosagei,2,duri,2,Si,2>,…,<medi,M,dosagei,M,duri,M,Si,M>} represents medication information, where medi,m∈{0,1} indicates whether a medication m is taken by patient Pi. If medi,m=1, we record dosage dosagei,m, duration duri,m, and the starting date Si,m; otherwise, dosagei,m, duri,m, and Si,m are set to 0. AGEi and SEXi denote patient age and gender, respectively. Yi,k,k∈{1…5} are binary indicators of decompensation symptoms (k=1,2,3,4) and death (k=5) for patient Pi.

The main prediction task is formally defined as follows: Given patient comorbidities D, medications M, and demographics (AGE and SEX), a predictive method learns patient embeddings to estimate symptom-specific decompensation risks (ascites, varices, HE, and HRS) within a 3-year window after the cirrhosis diagnosis. For each patient, a prediction is made on the date of cirrhosis diagnosis. That is, all predictor values (demographics, comorbidities, and medication variables) are extracted from patient data on or before the cirrhosis diagnosis date. Outcomes indicate whether distinct decompensation events (ascites, varices, HE, and HRS) occurred within a 3-year window after the cirrhosis diagnosis.

The proposed method makes use of an auxiliary task, mortality estimate, for symptom predictions. Ascites, varices, HE, and HRS mark the onset of decompensation; death is a critical event in a patient’s trajectory. Consideration of mortality as an auxiliary task provides a severity signal that regularizes the shared representation. It also allows the method to learn essential feature patterns and symptom constellations for both timely complication and lethality estimates, thereby increasing its ability to distinguish milder from rapidly progressive deteriorations.

Figure 1 presents the architectural framework of the proposed method, which has 3 important components: group-specific comorbidity graph construction, group- and individual-level graph representation learning, and multitask estimates. The group-specific comorbidity graph construction stratifies patients by demographics, such as age, gender, and liver condition indicators. For each patient group, it builds a disease relationship graph in which edge weights reflect the group’s mutual comorbidity co-inference. Patient-level graph representation learning obtains graph embeddings from each group-specific graph and personalizes them at the individual level by integrating each patient’s temporal comorbidity. The multitask estimates component then uses an attention mechanism to feed the personalized embeddings as key and value vectors into a multihead attention to perform the predictions. We describe each key component next.

This graph construction component creates distinct graphs for different demographic cohorts to capture how comorbidities co-occur and interact within each patient group. To create group-specific comorbidity graphs, patients are first stratified into different demographic cohorts, which could enhance downstream predictions because demographic similarity relates to decompensation trajectories [36,37]. For each cohort h, we created a comorbidity graph Gh=(Vh,Eh), where Vh indicates nodes (comorbidities), and Eh denotes edges between 2 comorbidities that coexist in at least one patient within the cohort. For a patient group Gh, we quantified the edge weight between 2 comorbidities p and q based on a CIS, calculated as follows:

With graph Gh, we applied Node2Vec [88] to represent nodes and edges in a way that preserves both local and global disease relationships. For each disease lj in graph Gh, Node2Vec outputs an embedding vector ϕh(lj )∈Rv. To account for disease temporal dynamics at the patient level, we constructed individual-level embeddings:

where vi,j is the individual-level embedding of comorbidity j for patient Pi, and δi,j denotes the time since diagnosis of the comorbidity, which serves as a proxy for disease chronicity. The individual-level embeddings were constructed using the group-level embeddings by incorporating disease duration to distinguish long-standing risk factors from acute conditions, thereby enhancing individual-level risk stratification for decompensation predictions.

As broad demographic trends, such as age-gender strata and overall disease impact, modulate decompensation risk beyond any single comorbidity, we endowed each cohort h with a compact “group fingerprint” gemb(Gh)∈Ru, obtained by applying principal component analysis (PCA) to the V-dimensional Node2Vec space:

where U∈Rv×u is the matrix of top-u PCA eigenvectors, ϕh_ is cohort h’s mean node representation, and Bh is the total number of important comorbidities in cohort h. gemb(Gh) offers a low-rank projection that can reduce noise, curb overparameterization, and provide the network with a succinct summary of cohort-wide factors [89], which individual-level embeddings alone cannot encode (eg, shared lifestyle and care access patterns).

We then concatenated patients’ individual-level embedding vi,j with a group fingerprint gemb(Gh) to form a combined context tensor.

where V~i represents the concatenated context tensor for patient Pi, and B denotes the total number of comorbidities in Gh(i). This effort integrates both individual disease dynamics and group-level characteristics, thus providing a more comprehensive representation of each patient’s profile.

As group parameters, PCA was exclusively fit on the training-set embeddings; testing-set embeddings were transformed using the fixed, training-fitted parameters (centering and projection) without refitting. Figure 2 presents the pipeline of patient embedding learning and use, which comprises IF score learning, subgroup graph constructions, and individual-level embedding generation, which jointly support downstream prediction tasks.

A binary mortality indicator was included as an auxiliary task for accurate decompensation symptom predictions because mortality is an important clinical event associated with decompensation [28]. The proposed method adopts a multitask architecture to predict 4 decompensation symptoms and mortality within a 3-year window. Mortality can inform (severe) prognosis representation learning, thereby improving decompensation predictions. Inclusion of this auxiliary task is aligned with clinical knowledge that recognizes the relationship of decompensation progression and mortality. By incorporating mortality in the modeling, we have an additional means to distinguish patient trajectories that vary in diseases [90,91]. We construct an inputXi=[AGEi,SEXi,{di,j,δi,j}j=1D,{medi,m,dosagei,m,duri,m,Si,m}m=1M,v^i] for the multitask prediction component, which contains patient Pi’s demographics, comorbidities, medications, and the learned patient context vector. Xi is passed to a fully connected layer and converted to a latent vector E:

The multitask prediction component fuses the tabular latent vector with the personalized graph representation into a single context vector ci and then applies multihead self-attention to obtain an attended representation for prediction.

where fgraph(⋅) is the graph representation module mapping the per-disease context V~i to a fixed-length vector, and Ei` is input to a fully connected layer to produce a shared representation (task-specific heads):

A linear head predicts mortality as an auxiliary task. A shared multilayer perceptron head takes the concatenated vector of Ei, refined and predicted death probability y^i,DEADT as input and projects it to logits for ascites, variceal bleeding, HE, and HRS. For patient i and task k, we computed the following:

ASCT, VART, ENCET, and HEPT represent the 4 decompensation symptoms, ascites, varices, HE, and HRS, respectively, while DEADT represents death. Next, we applied a sigmoid activation on y^i,k to produce an estimated probability for each symptom. The auxiliary task regularizes model learning by steering the network toward features and trajectories associated with severe clinical outcomes. The auxiliary task helps the method learn patterns indicative of disease progression and thereby differentiates high- and low-risk trajectories more effectively for the main (decompensation) prediction task.

The proposed method was trained with binary cross-entropy (BCE) on each target yk:

where lBCE(y^i,k,yi,k) is a BCE loss, IF(lj,yk) is the IF score measuring the influence of patient Pi’s comorbidity to target yk, and di,j is a binary indicator of each comorbidity. λ∑j=1B(di,j⋅IF(.,.)) serves as a scaling factor, adjusted by the influence of comorbidity, with λ controlling the strength of the scaling. Among the targets, the mortality head uses the standard BCE without IF modulation because its risk is already modeled through the symptom heads implicitly. In addition, we applied standard positive-class weights [92] to address data imbalances (ie, assigning a higher weight to the infrequent HRS label) and improve the model performance among less frequent classes. We present key hyperparameters in Multimedia Appendix 1.

From a major health care organization in Taiwan, we obtained demographic and clinical data for 29,915 patients diagnosed with cirrhosis between 2007 and 2013. Each patient was described by variables pertaining to demographics, comorbidity, and medication (Multimedia Appendix 2). We considered 4 outcome classes (ascites, varices, HE, and HRS) and an auxiliary task (mortality). To ensure comparability across patients, we applied z-score standardization to all disease durations [93]. Table 2 presents patient demographics and descriptive statistics, as well as the outcome class distribution. As shown, patients with HRS were slightly older (mean 57.62, SD 13.017 years) than those with other decompensation symptoms, with the varices group being the youngest (mean 55.88, SD 12.58 years). In all 4 outcome classes, approximately two-thirds of patients were men, suggesting a relative gender imbalance in the sample. Among the symptoms, patients with HRS had the most lingering conditions in their deterioration path to decompensation (239/273, 87.56%) and patients with varices the least (4010/6339, 63.26%). In the sample, 17,464 (58.38%) patients were deceased, and 12,451 (41.62%) were alive.

For each patient, variable values on or before the cirrhosis diagnosis date were used for model training, fine-tuning, and testing. Each outcome class was coded as a binary variable, with 1 indicating its occurrence within the 3-year window and 0 otherwise. To evaluate the proposed method, we split the sample into training (80%) and held-out testing (20%), with stratified datasets preserving class proportions across different outcome classes [94]. We repeated this procedure 10 times, using a different random seed each time. We report each investigated method’s performance averaged across the 10 trials. We also conducted 2-tailed t tests to examine whether the observed performance differentials were statistically significant. The evaluation relies on recall, precision, F-measure, and area under the curve (AUC). We report macro-averaged values because they reflect overall performance across different outcome classes, with the consideration of the imbalanced distributions in the sample. By and large, F-measure and AUC are more indicative of a method’s predictive utilities than recall and precision [95]. For probability-to-label conversion, we adopted appropriate symptom-specific decision thresholds. Because HRS is highly imbalanced, we increased the weight of the positive class 10-fold for the proposed method and all benchmarks. For this outcome class, the decision threshold was selected separately for each method that used the 5-fold cross-validation on the training data. We evaluated thresholds from 0.05 to 0.95 in increments of 0.05 and selected the threshold that produces the highest F₁-score. This procedure suggests a threshold of 0.85 for the proposed method. For ascites, varices, and HE, we applied a threshold of 0.50. The held-out test data were not used for threshold selection or class calibration.

We considered 12 prevalent methods for benchmarks. Some methods can be applied for both single-task and multitask predictions, so we used “multimethod” to refer to the multitask setting. These benchmarks include 4 conventional machine learning methods (XGB, ADAB, RF, and linear regression), 2 deep learning–based methods (single-task Single-fully connected neural network [FCN] and multitask Multi-FCN), 2 CNNs (single-task Single-CNN and multitask Multi-CNN), 2 TCNs (single-task Single-TCN and multitask Multi-TCN), and 2 LSTM networks (single-task Single-LSTM and multitask Multi-LSTM). These methods represent different analytics approaches and exemplify prevalent methods. All benchmarks used the same input variables and adopted the identical training (80%) and testing (20%) split to ensure consistency in the comparison. For each method, 5-fold cross-validation was applied on training data for parameter tuning. We describe model configurations, key hyperparameter settings, tuning ranges, and procedure in Multimedia Appendix 3.

For the proposed method, we created different demographic patient cohorts. In line with Liu et al [37], we stratified patients by gender and age (≤40, >40 and ≤60, >60 years) to reveal age distributions in the sample. We also designed a binary liver condition indicator that incorporates clinical and lifestyle factors (eg, biliary disease, Wilson disease, lack of medication use, overweight, smoking, and prior mild or severe liver disease or transplantation). This binary indicator consolidates important risk factors into a single feature to reduce sparsity and improve statistical power [96]. Figure 3 presents an exemplary graph to illustrate different diseases and their co-occurrence for a specific patient cohort (group). As shown, the edge between HBV and mild liver disease (MLD; liver damage and coma from hepatitis and unspecified viral hepatitis) has a weight of 2.85, reflecting a strong relationship between these comorbidities for this patient group, which influences decompensation risk. Parameters of the cohort-specific graph were learned from training data and were fixed at testing time. We summarize the different cohorts together with their size in Multimedia Appendix 2. To reduce the effects of varying cohort sizes, the cohort-level embeddings were combined with shared multitask representations and individual-level features to enhance prediction.

This study was approved by the Chang Gung Medical Foundation Institutional Review Board (201800270B0). Data were deidentified prior to any analyses, in compliance with established ethical standards for health care research. All research procedures were performed in accordance with relevant guidelines and regulations. No consent was required from patients in the sample. Physician participation in the interview was voluntary; each physician was explicitly informed of the voluntary nature, with his or her right to withdraw at any time preserved. Physicians received no financial compensation for their participation in the interview.

Table 3 presents the average performance of each method across 10 trials, together with the 95% CIs. As shown, the proposed method outperformed all benchmarks in most performance measures, exhibiting superior capability to estimate patients’ risk of decompensation symptoms. It attained an AUC value of 0.871 (SD 0.004), showing a 6.09% improvement over the best-performing benchmark (XGB). The proposed method also achieved an F-measure value of 0.678 (SD 0.010), which exhibits a 13.95% improvement over the best-performing benchmark (Single-FCN). It had a recall value significantly higher than that of the best-performing benchmark, a 15.42% improvement over Single-FCN, although its precision was marginally lower than that of the best-performing benchmark (0.663, SD 0.005 vs 0.675, SD 0.131 by ADAB). These comparative results affirm the utilities of this graph-based, multitask, deep learning method for predicting patients’ deterioration from the compensated to the decompensated stage.

We present the respective methods’ predictive performance for each decompensation symptom in Multimedia Appendix 4. As shown, the proposed method achieved the highest AUC values across all symptoms and had the highest recall and F₁-scores for all symptoms except HRS, manifesting its superior predictive performance. For HRS, it had an AUC value of 0.871 (SD 0.012), which was 3.81% higher than that of the best-performing benchmark (0.839, SD 0.013, Single-TCN). Its F-measure value was 0.518 (SD 0.018), comparable to that of the best-performing benchmark (0.533, SD 0.011, Single-LSTM). For varices, the proposed method had an F-measure value of 0.744 (SD 0.006), which was 14.29% higher than that of the best-performing benchmark (0.651, SD 0.017, Single-FCN). It had an AUC value of 0.877 (SD 0.003), exhibiting a 9.90% improvement over the best-performing benchmark (0.798, SD 0.006, RF). Its recall was 22.87% higher than the best-performing benchmark (0.656, SD 0.006, Multi-TCN), and its precision was 2.67% higher than the best-performing benchmark (0.712, SD 0.019, RF). For ascites, the proposed method attained an F-measure of 0.704 (SD 0.006), showing a 7.98% improvement over the best-performing benchmark (0.652, SD 0.019, Single-FCN). It had an AUC value of 0.843 (SD 0.005), which was 2.55% higher than that of the best-performing benchmark (0.822, SD 0.004, XGB). The proposed method had a recall value of 0.762 (SD 0.006), which was 18.14% higher than that of the best-performing benchmark (0.645, SD 0.030, Single-FCN), although its precision was lower than that of the best-performing benchmark (ADAB). For HE, the proposed method attained an F-measure value of 0.744 (SD 0.006), which was 29.62% higher than the best-performing benchmark (0.574, SD 0.070, Single-FCN). Its AUC value was 0.894 (SD 0.004), showing a 7.32% improvement over the best-performing benchmark (0.833, SD 0.003, XGB). The proposed method’s recall (0.818, SD 0.011) was notably 42.51% higher than that of the best-performing benchmark (0.574, SD 0.060, Single-FCN), although its precision was lower than that of the best-performing benchmark (ADAB).

These results further indicated that the proposed method can estimate decompensation symptoms among patients with cirrhosis more effectively than many prevalent methods. It attained the highest F-measure values for ascites (0.704, SD 0.006), varices (0.744, SD 0.006), and HE (0.744, SD 0.006). For HRS, which is particularly difficult to predict, the proposed method achieved the best AUC value of 0.871 (SD 0.012) while maintaining an F-measure value comparable to the best-performing benchmark. The prediction utilities, overall or for each decompensation symptom, affirm the proposed method’s greater efficacy for identifying patients at higher risk of decompensated cirrhosis than existing methods.

We conducted an ablation analysis to examine each important component’s contribution to the proposed method’s predictive utilities. The ablation considered 5 models: M4T, M4t+AUX, M4t+AUX+Graph, M4T+AUX+Graph+Personal, and M4T+AUX+Graph+Personal+Loss (the complete proposed method). All models followed the same procedure as the main evaluation for training and evaluation (ie, identical patient cohort, input variables, data preprocessing, training-testing split, and performance measures). M4T served as a baseline, with a multitask deep learning framework for all 4 decompensation symptoms (ascites, varices, HE, and HRS). M4T+AUX extends M4T by including an auxiliary mortality task (DEADT) to provide shared severity supervision. M4T+AUX+Graph further incorporates group-specific comorbidity graph representation (ie, Node2Vec-derived disease embeddings from different demographic patient groups). M4T+AUX+Graph+Personal integrates individual-level comorbidity temporal patterns into the graph embeddings. M4T+AUX+Graph+Personal+Loss denotes the proposed method in its entirety, containing all the key components.

Table 4 presents the respective models’ average AUC values across the 4 decompensation symptoms. As shown, the multitask baseline (M4T) had an AUC value of 0.803. With the inclusion of the auxiliary mortality task, M4T+AUX attained 0.811, showing a 1.00% improvement over the baseline model. Adding group-level graph representations, M4T+AUX+Graph increased the AUC value to 0.839, exhibiting a 3.45% improvement over M4T+AUX. Adding patient-level temporal personalization, M4T+AUX+Graph+Personal further increased the AUC value to 0.861, exhibiting a marginal increment of 2.62%. Finally, the complete proposed methods attained an AUC value of 0.871 by incorporating IF-based loss reweighting, showing a marginal increment of 1.16%, relative to M4T+AUX+Graph+Personal. The results indicate that each component contributed to the proposed method’s performance in a comparable way, with graph representation and temporal personalization offering greater incremental performance gains, while the auxiliary prediction task and IF-weighted optimization yielded additional improvements.

Beyond prediction performance, the proposed method’s estimates provide post hoc interpretability, unfolding essential patient-, disease-, and task-level factors that explain its estimates. We present 4 post hoc analyses: a contrastive comorbidity graph analysis, a Shapley Additive Explanations (SHAP) analysis [97], a counterfactual disease risk analysis, and a decompensation symptoms interdependence analysis. To enhance probability interpretability, we applied post hoc probability calibration per outcome using Platt scaling (logistic calibration). For each target variable, a calibration model was adjusted to fit a held-out (20%) training subset, thereby mapping raw logits to calibrated probabilities. The calibrated probabilities were then used in the subsequent counterfactual disease risk analysis as well as the interdependence analysis of decompensation symptoms.

Figure 4 presents exemplary visualizations of a patient’s comorbidity graph, together with her group-specific comorbidity graph as a reference. As shown, the left panel displays the group-specific comorbidity graph for the reference patient cohort (female, age >40 and ≤60 years, and having prior liver issues), where edges are weighted by group-level disease co-occurrence (IF×group multiplier). The right panel depicts an individual patient’s network adjusted for chronicity, according to the temporal order and persistence of each condition. The edge weights reveal how different clinical conditions (hepatitis B, MLD, and nonchronic diabetes) of the focal patient have interacted over time and thereby unfold essential comorbidities. The focal patient’s network has fewer nodes and edges, reflecting that this patient has fewer comorbidities than the cohort. These illustrated networks allow physicians to compare a patient with a reference cohort and assess the chronicity and relationship of comorbidities to interpret the method’s estimates for clinical decision-making and patient management.

We also performed a SHAP analysis to scrutinize feature importance for each outcome class (ascites, varices, HE, and HRS). SHAP values were calculated for individual test patients using a training-derived background sample, and features were then ranked by the absolute magnitude of SHAP to identify the most influential contributors for each outcome. That is, SHAP was applied for explanation only instead of modifying model training or performance estimates. Figure 5 illustrates the proposed method’s predictions via SHAP-enabled visualizations. For consistency of interpretation, we used the same illustrative patient in Figures 4-7, who was randomly selected from testing data and had at least one positive decompensation symptom. As shown, top predictors pertain to the patient’s disease network, which reinforces the consistency between graph-derived variables and clinically important variables. In this example, ascites, varices, and encephalopathy predictions were mainly influenced by chronicity of hepatic deterioration and viral hepatitis, as reflected by MLD, HBV, HCV, diabetes, and hypertension features in Figure 5. The visualization of feature importance reveals influential factors and helps physicians better interpret the proposed method’s estimates at the patient level.

To understand the association between comorbidities and decompensation symptom risks, we performed a probabilistic counterfactual risk analysis using patient comorbidity profiles. We first computed baseline predicted probabilities for each decompensation symptom using each patient’s observed comorbidities and comorbidity durations. We then generated counterfactual inputs by removing one comorbidity at a time while keeping all other covariates fixed. Removal was proceeded by setting the focal comorbidity’s binary indicator, duration feature, and disease-specific node embedding to 0, with all remaining features and embeddings intact. For each removed comorbidity d and outcome yk, we recomputed risks and calculated the change of risk Δd,k as follows.

A Δd,k>0 means that removal of d increases the estimated probability of decompensation symptom k; Δd,k<0 indicates a decrease in risk. Figure 6 presents patient-level Δd,k values in a comorbidity-by-outcome matrix in which each row refers to a removed comorbidity and each column denotes a decompensation symptom. For the illustrative patient in Figure 6, removing MLD yields negative shifts across all outcomes (ascites: −0.036; varices: −0.008; HE: −0.002; and HRS: −0.073), with the largest magnitude for HRS, implying that, for this patient, the proposed method predicts positive associations of MLD with the estimated risks for different symptoms, particularly HRS. This outcome-specific decomposition is useful and can help physicians understand the predictions for each patient, although it is an associative interpretation of the learned predictive function and should not be interpreted as direct physiological causality or treatment effect. All percentage-point shifts (Δ) in Figure 6 were computed from calibrated probabilities.

We also examined the relationships among patient outcomes, mortality, and decompensation symptoms. Figure 7 presents a patient-level intertask association map that depicts how model-estimated risk for one outcome is linked to the estimated risks of others. These relationships can help clinicians identify outcomes strongly associated with others at the patient level. This analysis perturbs the probability of one outcome and quantifies its directional effect on other outcomes. For each outcome B∈{ASCT,VART,ENCET,HEPT,DEADT}, we fit a destination-specific logistic regression stacker with estimated probabilities of all other outcomes, XB={p^A:A∈{ASCT,VART,ENCET,HEPT,DEADT}∖{B}}, and response yB. Each stacker estimates P(yB=1|XB) as a calibration and association layer over the proposed method’s outputs. For each patient, we used XB to compute baseline stacked probabilities p^Bbase, then replaced an p^A with its training-cohort median p~A and kept all remaining predictors fixed to recompute p^Bcf. For each directed perturbation A→B, we defined the edge effect as ΔpA→B=p^Bcf-p^Bbase (in percentage points [pp]). A positive ΔpA→B reveals that neutralizing source signal A increases destination risk B, whereas a negative ΔpA→B indicates a decrease. In Figure 7, the logistic stackers were fitted on the training data and then applied to the illustrated case drawn from the test data to demonstrate their capabilities under deployment conditions.

Figure 7 illustrates how changes in the prediction probability of source outcome can lead to a reversal of the estimated destination outcome. For this patient, all 5 displayed effects are negative, indicating that replacing the source-outcome probability with its training cohort median decreases the estimated probability of the corresponding destination outcome: ΔpASCT→DEADT=−50.12 pp, ΔpASCT→VART=−42.60 pp, ΔpASCT→ENCET=−45.59 pp, ΔpVART→ASCT=−32.50 pp, and ΔpVART→ENCET=−45.91 pp. For example, ΔpVART→ENCET=−45.91 pp means that replacing the variceal-bleeding probability with its cohort median lowers the stacked encephalopathy probability by 0.459 for this patient. The illustrative patient-level dependency map can support triage and highlight destination risks most sensitive to other outcome signals, with an emphasis on the associations among patient outcomes. All reported percentage-point shifts (Δ) in Figure 7 were computed from calibrated probabilities.

Compared with many prevalent methods, the proposed method exhibited significantly better prediction efficacy, reflected by an average AUC value of 0.871 (SD 0.004) and F-measure of 0.678 (SD 0.010). This method is effective for symptom-level decompensation risk stratification beyond a dichotomous compensation versus decompensation view. An ablation analysis revealed the importance of each key component of the proposed method, perhaps with greater performance gains by graph-based representation and patient-level temporal patterns, which underscores the criticality of capturing comorbidity structure and disease chronicity for downstream predictions. Moreover, the proposed method can provide interpretable analyses and visualizations (such as patient-level comorbidity graphs, SHAP attributions, comorbidity-removal counterfactual maps, and intertask assist maps), which should increase user confidence, trust, and adoption.

The proposed method addresses several important limitations of prior work on liver-related risk predictions. For example, most existing methods target a single outcome, such as compensation versus decompensation, without symptom-level consideration. In addition, patient subgroup structure and comorbidity interplays are largely overlooked. The proposed method adds to extant literature by considering distinct decompensation symptoms and their interplays, together with an auxiliary task. It incorporates demographic cohort structure using group-specific comorbidity graphs, as well as patient-level temporal disease information through personalized embeddings. Relative to existing regression-based, tree-based, and deep learning–based methods, our method makes methodological contributions by integrating clinically guided structure into learning via IF- and CIS-informed graph construction and IF-weighted multitask loss. Furthermore, this study provides multiple post hoc analyses for interpretability and visualizations of symptom-specific associative patterns, so the proposed method can serve as a structured and clinically interpretable means to complement score-based methods for assessing decompensation risks among patients with cirrhosis.

This study has several limitations that should be considered to interpret the results with caution. First, we used patients from a health care system, who were diagnosed with cirrhosis between 2007 and 2013. This cohort reflects the pre–direct-acting antiviral for HCV, before a major practice change around 2014. Clinical care of patients with chronic liver diseases has evolved over time, such as antiviral access, vaccination coverage, and preventive measures that might affect patients’ decompensation risks, disease trajectories, and outcomes. The illustrative post hoc interpretability analyses may offer a partial depiction regarding how etiology shifts may affect predictions. Physicians can perform counterfactual risk analyses to assess how removal or attenuation of specific etiologic factors, such as viral hepatitis–related condition reduction through vaccination and modern antiviral treatment, might alter estimated probabilities for ascites, varices, HE, and HRS. Although such analyses are associative rather than causal, they still can help assess model sensitivity under plausible epidemiologic changes. Future work should enhance robustness and external validation with recalibration using additional patient data, desirable cohorts after 2020. Second, the proposed method formulates outcomes as fixed-horizon binary labels and optimizing them with task-wise binary cross-entropy, rather than using an explicit competing-risks or survival framework. Its inclusion of mortality prediction as an auxiliary task helps the main prediction task but without considering the temporal order of decompensation and mortality. Therefore, estimates should be interpreted as risk within the 3-year window, instead of cause-specific hazard or cumulative-incidence estimates. Future work can extend this method with time-to-event and competing-risks modeling, such as using a cause-specific or subdistribution hazard approach or a multistate survival model. Third, the proposed method was developed and evaluated using available etiological data in EHRs; it lacks laboratory and blood test data. Future research could extend the method by incorporating additional data for further improved, finer-grained predictions, such as modeling the temporal order of clinical events and adopting a stricter temporal view to formulate the prediction problem.

These limitations in turn suggest promising research directions. We used a particular patient cohort to develop and evaluate the proposed method, and continued research could consider multiple patient cohorts that differ in demographics, ethnic backgrounds, and clinical systems to increase both method generalizability and result robustness. Expanding the scope of inputs by including additional cirrhosis-related conditions such as jaundice [38] or laboratory biomarkers and multimodal data (eg, medical imaging and clinical notes) could further enhance clinical values [98]. Estimating the first-appearing symptom and its timing also is important for clinical decision support and patient management. Moreover, future research can extend the proposed method’s interpretability by incorporating intrinsically interpretable methods for increased reasoning and explanations, which can enhance transparency and physicians’ trust and clinical adoption [99].

Accurately predicting decompensation among patients with cirrhosis is clinically important and methodologically challenging. We developed and evaluated a graph-based multitask deep learning framework to estimate 4 decompensation symptoms (ascites, varices, HE, and HRS) within a 3-year prediction horizon after the cirrhosis diagnosis. This method incorporates demographic subgroup structure, comorbidity interplays, patient-level temporal information, and IF-weighted multitask optimization for downstream predictions. The proposed method provides more effective predictions than many existing methods, across different performance measures, thereby enabling early at-risk patient identification and risk stratification. The illustrated post hoc interpretability offers patient-level, associative explanations of its estimates for user comprehension and confidence.