This retrospective multicenter cohort study included patients who initiated dialysis between January 2014 and December 2020 in the nephrology and hemodialysis departments of four clinical centers in China.

Inclusion criteria included: (1) patients aged 18–75 years who newly started maintenance dialysis, and (2) patients who underwent a non-contrast chest or abdominal multidetector CT scan that included the L1 level within 1 month of dialysis initiation. Exclusion criteria were: (1) acute infection during the peridialytic period, (2) malignant tumors, (3) hepatic failure, (4) conditions that impair intestinal nutrient absorption such as inflammatory bowel disease, chronic diarrhea, or short bowel syndrome, (5) kidney transplantation or withdrawal from dialysis, (6) a change in dialysis modality, and (7) lack of valid follow-up records.

As shown in Figure 1, a total of 3929 initial dialysis patients were screened from the nephrology and dialysis units of four tertiary hospitals in China. Among them, 1820 patients met the inclusion criteria. After applying the exclusion criteria, 1051 patients were finally included in this study.

This study was approved by the Ethics Committee of Zhongda Hospital (approval number 2022ZDSYLL003-P01) and was registered in the Chinese Clinical Trial Registry (registration number ChiCTR2300068453). The study complied with the principles of the Declaration of Helsinki. Because of the retrospective study design, the ethics committee waived the requirement for written informed consent.

All study data were collected by trained research personnel. Data collection procedures and equipment were standardized across the four study sites. Demographic and clinical information at enrollment was obtained from patients’ medical records, including age, sex, height, weight, smoking history, alcohol use, dialysis modality (hemodialysis or peritoneal dialysis), history of diabetes, hypertension, coronary artery disease, chronic heart failure, stroke, cardiac intervention, CVD, hyperlipidemia, and anemia, as well as the use of beta blockers, angiotensin converting enzyme inhibitors or angiotensin receptor blockers, calcium channel blockers, diuretics, erythropoietin, iron supplements, antiplatelet agents, Compound α-keto acid, and glucocorticoids. Body mass index (BMI) was calculated as weight in kilograms divided by height in meters squared.

White blood cell count (WBC), hemoglobin (Hb), platelet count, albumin (ALB), fasting plasma glucose (FPG), uric acid (UA), triglycerides (TG), total cholesterol (TC), low density lipoprotein cholesterol (LDL-C), serum creatinine (SCr), cystatin C (CysC), the ratio of serum creatinine to cystatin C (SCr/CysC), and blood urea nitrogen were measured within 1 week before the initiation of dialysis using standard laboratory methods.

CT examinations were performed using Discovery CT750, Revolution CT, and Optima CT 660 scanners from GE Healthcare in Milwaukee, Wisconsin, the SOMATOM Sensation scanner from Siemens Healthineers in Erlangen, Germany, or the Ingenuity CT system from Philips in Amsterdam, the Netherlands. All scans were acquired with standard settings, including 120 kVp, automatic dose modulation systems (automA and smartmA for GE Healthcare scanners, CareDose 4D for Siemens Healthineers scanners, and DoseRight for Philips scanners), a 512 × 512 matrix, a collimation of 0.625 mm, and a slice thickness of 5 mm. CT images were imported into ImageJ software (version 1.46, developed by the National Institutes of Health). The L1 level was first identified, and the axial slice showing the largest transverse diameter of the L1 transverse processes was selected for analysis. Tissue areas were then identified and quantified according to preset Hounsfield Unit ranges. The radiodensity of each tissue type was expressed as the mean CT attenuation within its corresponding HU range. The HU thresholds used to classify different tissues followed commonly accepted values in published studies. Skeletal muscle was defined as −29 to +150 HU (Golder et al., 2022). Low attenuation muscle was defined as −29 to +29 HU (Kim et al., 2020). Subcutaneous fat was identified using a threshold of −190 to −30 HU, visceral fat using −150 to −50 HU (Weinberg et al., 2018), and total fat tissue using −190 to −30 HU (Crawford et al., 2020).

The primary endpoint of this study was CVD death. CVD death was defined as death caused by acute coronary syndrome, sudden cardiac death, life-threatening arrhythmias, congestive heart failure, or ischemic/hemorrhagic stroke (Tsai et al., 2020; Sato et al., 2016). Sudden cardiac death was defined as cardiac arrest occurring within 1 hour after the onset of symptoms (Ramesh et al., 2016). Life-threatening arrhythmias were defined as documented ventricular tachycardia or ventricular fibrillation. Congestive heart failure was confirmed by electrocardiography, chest radiography, or echocardiography together with symptoms such as dyspnea or edema. Stroke was diagnosed based on characteristic imaging findings and clinical examination. All enrolled patients were followed from the date of dialysis initiation. Follow up ended at the earliest occurrence of the study endpoint, death, withdrawal from the study, loss to follow up, or the end of the follow up period on 31 December 2022. To ensure consistency across centers, a standardized definition of CVD death and a unified data extraction protocol were applied at all participating sites. This research adhered to the Transparent Reporting of a Multivariable Prediction Model for Individual Prognosis or Diagnosis (TRIPOD) Statement (Moons et al., 2015).

In the training set, univariate logistic regression was first used to identify variables that might be considered for model development. Variables with a P value less than 0.05 in the univariate analysis for CVD death were regarded as significant and entered the subsequent selection process. To reduce overfitting and improve predictive performance, the Least Absolute Shrinkage and Selection Operator (LASSO) method was then applied for further variable selection (Tibshirani, 1996). This method applies a penalty to the regression coefficients, which helps address multicollinearity among variables and improves model simplicity by automatically removing predictors with weak statistical contribution or high collinearity (He et al., 2023). In addition, multicollinearity among the selected variables was assessed using the variance inflation factor, and a value less than 2 indicated no evident multicollinearity.

In this study, eight machine learning models were used to predict the risk of cardiovascular mortality in patients starting dialysis. These models included Support Vector Machine (SVM), Gradient Boosting Machine (GBM), Neural Network, Extreme Gradient Boosting (XGBoost), Adaptive Boosting (AdaBoost), Light Gradient Boosting Machine (LightGBM), Categorical Boosting (CatBoost), and Logistic Regression. In this study, the number of patients who experienced cardiovascular death was much smaller than the number who did not, resulting in an imbalance between the positive and negative samples that could affect the accuracy of the prediction models. To address this issue, the Synthetic Minority Over sampling Technique was used to balance the dataset (Chawla et al., 2002).

Model performance was evaluated in both the internal and external validation sets. Discrimination was assessed using the AUC. Classification performance was evaluated using accuracy, precision, sensitivity, specificity, and the F1 score. Calibration was examined with calibration curves, and clinical usefulness at different threshold probabilities was assessed with decision curve analysis. These complementary measures were used to provide a comprehensive evaluation of the predictive ability and clinical applicability of each model. Based on performance across these metrics in the training and testing sets, we selected the model with the best overall predictive performance.

Model interpretability was evaluated using the Shapley Additive Explanations (SHAP) method. SHAP is a model independent approach based on game theory that can be used to assess the overall importance of features and to explain predictions for individual samples. This method helps reduce the black box nature of machine learning models, improves transparency, and provides insight into how each input variable contributes to the predicted outcome. It also allows a deeper understanding of the model’s decision process and the interactions among features (Lundberg et al., 2018).

To enhance the clinical applicability of the prediction model, we developed an online calculator based on the final model. This web-based tool allows users to enter relevant clinical information and obtain real time estimates of the predicted probability of the outcome. The calculator was implemented using the Shiny framework in R, with the user interface built through the ui.R file and the server logic handled by the server.R file. The trained model object was loaded to generate predictions. The final application was deployed on the shinyapps.io platform and can be accessed from any device with an internet connection.

The extent of missing data for each variable in the original dataset is summarized in Supplementary Table S1. Variables with missing rates exceeding 25% were excluded, and the remaining missing values were handled using multiple imputation. Normally distributed continuous variables were summarized as means with standard deviations, whereas skewed continuous variables were presented as medians with interquartile ranges (IQR). Categorical variables were reported as frequencies and percentages. Between group comparisons for continuous variables were performed using the independent samples t-test or the Mann-Whitney U test, depending on distribution. Categorical variables were compared using the chi-square test. Given the limitations of P values in detecting differences between groups in large sample studies, the standardized mean difference was included as an additional measure to evaluate the magnitude of group differences. As a standardized effect size, the standardized mean difference is not influenced by sample size, measurement scale, or variance, and therefore provides a more objective assessment of group imbalance. An absolute standardized mean difference below 0.20 was considered small, and a value below 0.10 indicated that the difference was negligible (Cohen, 1988). All statistical analyses were performed using R software version 4.3.1 (R Foundation for Statistical Computing, Vienna, Austria) or STATA version 16.0 (StataCorp LLC, College Station, Texas, United States). A two-sided P value less than 0.05 was regarded as statistically significant.

The training and internal validation cohorts consisted of 645 patients, including 452 patients in the training set and 193 patients in the internal validation set. The median age was 55 years (IQR, 45–65). Overall, the median follow-up time was 44.6 months (IQR, 29.2–63.0). During follow up, 94 patients (14.6%) experienced CVD death. Compared with patients who did not experience CVD death, those who did were older and had lower levels of Hb, SCr, SCr/CysC, and SMD. They also had higher levels of low attenuation muscle area, the low attenuation muscle to skeletal muscle area ratio, subcutaneous fat area, and total fat area (Table 1; Supplementary Figure S1). In addition, the prevalence of diabetes, coronary artery disease, chronic heart failure, and CVD was higher in patients who experienced CVD death, as were the proportions of those with a history of cardiac intervention and those using antiplatelet agents (Table 1; Supplementary Figure S2).

The external validation cohort consisted of 406 patients, with a median age of 52 years (IQR, 42–63). During follow up, 40 patients (9.9%) experienced CVD death. Detailed patient characteristics are shown in Supplementary Table S2. Differences between the development dataset and the external validation dataset are summarized in Supplementary Table S3. Several variables differed significantly between the two datasets, such as age, sex, BMI, smoking and alcohol history, WBC, ALB, UA, LDL-C, CysC, SCr/CysC, SMI, SMD, low attenuation muscle area, low attenuation muscle density, the low attenuation muscle to skeletal muscle area ratio, subcutaneous fat area, subcutaneous fat density, visceral fat area, visceral fat density, total fat area, total fat density, diabetes, coronary artery disease, history of cardiac intervention, CVD, and the use of iron agents, antiplatelet agents, Compound α-keto acid, and glucocorticoids.

Using univariable logistic regression analysis, 17 variables associated with cardiovascular death in initial dialysis patients were identified, all with P values less than 0.05 (Table 2). LASSO regression was then applied to further select key predictors, resulting in eight variables: age, diabetes, CVD, history of cardiac intervention, dialysis modality, SMD, Hb and SCr (Figure 2). Multicollinearity was evaluated for these eight variables. Spearman correlation analysis showed correlation coefficients below 0.4, and all variables had tolerance values greater than 0.6 and variance inflation factors below 1.5, indicating no significant multicollinearity (Supplementary Figure S3; Supplementary Table S4). These eight variables were therefore used to construct the prediction model for CVD-related mortality in initial dialysis patients.

After the predictive features were selected, eight machine learning models were developed, including SVM, GBM, Neural Network, XGBoost, AdaBoost, LightGBM, CatBoost, and Logistic Regression. In the internal validation set, the AUC of the eight models was compared, as shown in Figure 3A. The CatBoost model demonstrated the best discrimination, with an AUC of 0.843 (95% CI, 0.766–0.919), followed by SVM (AUC 0.833; 95% CI, 0.756–0.910) and Logistic Regression (AUC 0.832; 95% CI, 0.755–0.909). Performance metrics for the internal validation set are presented in Table 3. Among all models, CatBoost showed the most favorable overall performance, achieving the highest accuracy (0.812), precision (0.778), and F1 score (0.824), with sensitivity (0.875) and specificity (0.750) also at comparatively high levels. The calibration curves for the internal validation set are shown in Figure 4. The predicted probabilities of the CatBoost model aligned closely with the reference line, indicating good calibration and suggesting that the model captured the individual risk levels with higher accuracy. Figure 5A presents the decision curves for predicting CVD mortality in the internal validation set. Across most threshold probabilities, the CatBoost model provided greater net benefit than the other models, demonstrating stronger clinical usefulness, particularly in settings requiring individualized risk assessment.

The ROC curves of the eight machine learning models in the external validation set are presented in Figure 3B. The CatBoost model showed the best discrimination with an AUC of 0.799 (95% CI, 0.729–0.869). The Neural Network model demonstrated comparable performance, with an AUC of 0.799 (95% CI, 0.731–0.867). Table 4 summarizes the predictive performance of each model in the external validation set. The CatBoost model achieved the highest accuracy (0.775), specificity (0.812), precision (0.797), and F1 score (0.766), suggesting strong generalizability and stable predictive performance on external data. The Neural Network model showed better sensitivity (0.825), although its accuracy (0.750), specificity (0.675), precision (0.717), and F1 score (0.767) were lower than those of the CatBoost model. Calibration curves for the external validation set are presented in Figure 6, showing that the CatBoost model exhibited calibration closer to the ideal diagonal line. Decision curve analysis for all models is shown in Figure 5B, and across most threshold probabilities, CatBoost provided a higher standardized net benefit than the other models. Given its superior performance over other models and consistent results in both validation sets, CatBoost was selected as the optimal model for interpretability analysis.

To further clarify predictive performance across follow-up durations, we evaluated the discrimination of the CatBoost model for CVD mortality at 1, 3, and 5 years in the external validation cohort. As shown in Figure 7, the AUC was 0.700 (95% CI, 0.540–0.861) at 1 year, 0.828 (95% CI, 0.786–0.870) at 3 years, and 0.845 (95% CI, 0.810–0.879) at 5 years. The 1-year estimate had a wider confidence interval, possibly due to fewer events in the first year, while discrimination was more stable at 3 and 5 years.

We applied the SHAP approach to explain the output of the final model by quantifying the contribution of each variable. In the global explanation, the SHAP summary bar plot (Figure 8A) showed that CVD, SMD, and Hb were the three most influential predictors, indicating their substantial contribution to the estimated risk of CVD mortality. In addition, variables such as dialysis modality, diabetes, and a history of cardiac interventions also contributed to model decisions. In Figure 8B, each point represents an individual patient, with its position on the horizontal axis indicating the SHAP value for that variable and the color reflecting the magnitude of the variable (yellow for higher values and purple for lower values). For example, CVD showed a clear positive influence on adverse outcomes, meaning that patients with CVD were more likely to be classified by the model as having a higher risk of CVD mortality. The importance and direction of effect of each variable are visualized through its ranking and SHAP value distribution. Similarly, SMD and Hb had strong effects on the predictions but in the opposite direction. Lower levels of SMD or Hb were generally associated with positive SHAP values, suggesting that reduced muscle density and lower Hb levels may be important contributors to an increased risk of CVD mortality. This SHAP scatter plot illustrates how variable levels relate to predicted risk and enhances the interpretability of the model by providing a visual link between clinical features and prediction outcomes. A SHAP value heatmap (Supplementary Figure S4) was used to visualize individual-level contributions of each predictor across the cohort, with patients ordered by the model-predicted outcome. Each column represents a patient and each row a variable, with color denoting SHAP values. Patients predicted to experience CVD mortality (orange section) exhibited stronger positive SHAP values (yellow tones) for variables such as CVD, SMD, and Hb. This pattern further supports the stable and consistent contribution of these predictors to the model’s predictions.

The relationship between Hb levels and their corresponding SHAP values is shown in Supplementary Figure S5A, with SMD encoded by color to further illustrate its influence on the model output. Low Hb levels were associated with positive SHAP values, with a more pronounced effect observed among patients with lower SMD. To explore this interaction pattern in greater detail, Supplementary Figure S5B plots SMD on the horizontal axis with Hb represented by color. A similar trend was observed: patients with lower SMD and lower Hb were more likely to be classified as high risk by the model. These findings highlight a potential synergistic effect of SMD and Hb in distinguishing patients at higher risk of CVD mortality.

Figure 9 presents the SHAP force plots for four patients in the testing set, illustrating how the CatBoost model generated individualized predictions by showing the specific contribution of each variable. In these plots, yellow bars represent features that pushed the prediction toward higher risk, whereas purple bars indicate features that pushed the prediction toward lower risk. The length of each bar reflects the magnitude of the feature’s effect on the model output. The value f(x) denotes the predicted risk for that patient, and E [f(x)] represents the model’s baseline value. As shown in Figure 9A, the patient’s predicted risk was primarily driven by several risk factors with positive contributions, including a history of CVD (CVD = 1, contribution +0.144), low SMD (16.4, contribution +0.119), low Hb (87, contribution +0.0776), and hemodialysis modality (dialysis modality = 1, contribution +0.0258). Although a few variables had minor negative contributions, they were insufficient to offset the cumulative influence of these risk factors. The final model output for this patient was f(x) = 0.886, which exceeded the baseline value E [f(x)] = 0.501, indicating a high predicted probability of CVD mortality. In contrast, the force plot for another patient shown in Figure 9C indicates that several variables contributed negatively to the model output, including higher SMD (47.8, contribution −0.173) and the absence of CVD (CVD = 0, contribution −0.151). The combined influence of these protective factors substantially lowered the predicted risk. The final model output for this patient was f(x) = 0.100, which was well below the baseline value and suggested a low probability of CVD mortality. Figures 9B,D similarly demonstrate that key variables such as CVD, Hb, and SMD consistently shaped predictions across different patients, showing good discriminative ability and stable interpretability and supporting the reliability of the model for individual level risk assessment.

Figure 10 shows the SHAP waterfall plots for the corresponding individuals, providing a clearer view of the order and direction of each variable’s cumulative contribution to the prediction. The variables are arranged according to their contribution, beginning at the model’s baseline value E [f(x)] and moving to the right for positive contributions and to the left for negative contributions, resulting in the final prediction f(x). For example, in Figure 10A, the presence of CVD, low SMD, and low Hb markedly increased the predicted risk. In Figure 10B, higher Hb and the absence of diabetes reduced the predicted value f(x), thereby lowering the estimated risk. Similar contribution patterns were observed in Figures 10C,D, indicating that the model maintained consistent interpretability at the individual level.

As shown in Figure 11, the final predictive model was incorporated into a web-based application to facilitate use in clinical practice. Users can enter the actual values of the eight required features, and the system automatically computes the predicted risk of CVD mortality for initial dialysis patients. The online predictive tool is accessible at the following link: https://wangxiaoxu0817.shinyapps.io/workrun15/.