This retrospective cohort study used data from hospitalized patients at Hainan General Hospital between January 2020 and December 2023. The Hainan General Ethics Committee reviewed and approved the study (Approval No. EC-YLY-2026-83-01). The study population comprised hospitalized patients who met the diagnostic criteria for liver cancer during the study period and contemporaneous hospitalized non-liver cancer controls. The control group included patients admitted for liver-related and non-liver-related conditions who had no liver cancer diagnosis in the electronic medical record. To improve the clinical interpretability of model specificity, we mapped the disease spectrum of non-liver cancer controls from electronic diagnosis codes and summarized it in Supplementary Table S2. This disease-spectrum analysis also clarified whether controls represented healthy individuals or non-liver cancer patients with liver disease or liver cancer-related risk backgrounds.

The diagnosis of liver cancer was required to conform to the relevant Chinese guidelines for the diagnosis and treatment of primary liver cancer (Zhou et al., 2020) and was confirmed by contrast-enhanced CT/MRI and/or pathological findings. In addition, the discharge diagnosis was required to include the corresponding code from the International Statistical Classification of Diseases and Related Health Problems, 10th Revision (ICD-10). The inclusion criteria were as follows: (1) age ≥18 years; (2) completion of relevant laboratory tests during hospitalization; and (3) complete clinical data. The exclusion criteria were as follows (Wang et al., 2025; Feng et al., 2026): (1) incomplete clinical or laboratory data; (2) other active malignancies; and (3) concomitant autoimmune diseases requiring long-term glucocorticoid or immunosuppressive therapy.

Patient data were extracted from electronic medical records from the hospital. The outcome measure was a confirmed diagnosis of liver cancer. The collected data included demographic characteristics (age and sex), laboratory indicators (complete blood count, liver function, renal function, coagulation function, inflammatory markers, and tumor markers), and clinical features (liver cirrhosis, ascites, and hepatic encephalopathy). To ensure data quality, all variables were extracted based on the original records from the electronic medical system and independently cross-checked by two researchers.

This study included 74 multidimensional features in three major categories: demographics, laboratory test results, and clinical characteristics. Outliers in continuous variables were identified using the 1.5 × interquartile range (IQR) rule. Specifically, values below Q1 − 1.5 × IQR or above Q3 + 1.5 × IQR were treated as outliers and set to missing before imputation. Features with a missing rate exceeding 40% were excluded from the analysis. For the remaining variables with missing data, continuous variables were imputed using multiple imputation by chained equations (MICE) with predictive mean matching (PMM), generating five imputed datasets with a maximum of 10 iterations. Missing categorical values were treated as a separate “missing” category where applicable. Missing-data handling was performed after dataset splitting and independently within the training and internal validation sets to reduce the risk of information leakage. Categorical variables were converted to factors, continuous variables were converted to numeric types, and variable names were standardized using the make.names function. Variables included in model development were required to have nonzero variance in the training set.

The dataset was split into training and internal validation sets at a 7:3 ratio using outcome-stratified random sampling with a fixed random seed. This split preserved similar liver cancer proportions in the two cohorts. Feature selection was performed exclusively in the training set using LASSO regression with 5-fold cross-validation (Xi et al., 2023). The minimum-error criterion (lambda.min; λ = 8.49 × 10⁻⁴), defined as the penalty parameter that minimized the cross-validated binomial deviance, was used rather than lambda.1se. Model terms with nonzero coefficients at lambda.min were retained for subsequent model construction. To evaluate collinearity, the correlation structure of continuous variables was examined using Spearman correlation (Supplementary Figure S1), and variance inflation factors (VIFs) were calculated for the retained terms. Feature-selection robustness was further assessed using bootstrap resampling of the LASSO procedure (Supplementary Figure S2).

The retained feature list, LASSO coefficients, VIF diagnostics, bootstrap sensitivity results, correlation heatmap, LASSO stability plots, and model-importance rankings are reported in the Section 3, Supplementary Table S3, Supplementary Figures S1, S2, and Figures 1–3 to improve transparency and reproducibility.

Six classification models were developed using the retained terms: logistic regression, random forest, support vector machine (SVM), k-nearest neighbors (KNN), extreme gradient boosting (XGBoost), and Elastic Net. Model tuning was performed within the training set using 5-fold cross-validation where applicable. Logistic regression was fitted using a binomial family. Random forest was configured with 500 trees. SVM used a radial basis function kernel with probability estimates enabled. For KNN, the optimal k value was selected by cross-validation. XGBoost was trained using the binary:logistic objective and AUC as the evaluation metric, with max_depth = 3, eta = 0.05, subsample = 0.7, colsample_bytree = 0.7, and gamma = 0.1. Elastic Net was configured with alpha = 0.5, and the optimal lambda was selected by cross-validation. All feature-selection and model-tuning procedures were conducted within the training set; the internal validation set was reserved for model evaluation.

The primary model-performance metric was the area under the receiver operating characteristic (ROC) curve. The 95% confidence intervals (CIs) for AUCs were calculated using the DeLong method (Zou et al., 2024). Pairwise DeLong tests were performed in the internal validation cohort to compare XGBoost with each of the other five models, and the resulting P values were adjusted using the Benjamini-Hochberg procedure. Additional performance metrics, including accuracy, sensitivity, specificity, and F1 score, were calculated at a prespecified probability threshold of 0.5. Positive predictive value (PPV) and negative predictive value (NPV) were also evaluated for XGBoost under the prespecified threshold and exploratory alternative threshold strategies. Because threshold selection influences clinical utility, decision curve analysis (DCA) was used to evaluate net benefit across threshold probabilities (Netto Flores Cruz and Korthauer, 2024). All metrics were calculated separately for the training and validation sets. SHapley Additive exPlanations (SHAP) was applied to the optimal model for global interpretation. SHAP values quantify feature contributions to individual predictions and were interpreted as model-based associations rather than evidence of biological causality.

All statistical analyses were performed using R 4.0.1 (https://www.r-project.org). Continuous variables were summarized as medians with interquartile ranges (IQRs). Normality was assessed using the Shapiro-Wilk test. Normally distributed continuous variables were compared using Student’s t-test, whereas non-normally distributed continuous variables were compared using the Mann-Whitney U test. Categorical variables were expressed as frequencies and percentages and compared using the chi-square test or Fisher’s exact test, as appropriate. Statistical significance was defined as a two-tailed P < 0.05.

A total of 9,284 hospitalized patients were enrolled in this study. The training set comprised 6,426 patients, including 3,130 (48.7%) with liver cancer, and the internal validation set comprised 2,858 patients, including 1,394 (48.8%) with liver cancer. This near-balanced case proportion facilitated model development and internal evaluation but differs from the prevalence expected in population-based screening settings. Diagnosis-code mapping showed that the non-liver-cancer control records were dominated by cirrhosis/chronic liver disease (2,167 records, 34.2%), viral hepatitis with cirrhosis (1,672 records, 26.4%), viral hepatitis without recorded cirrhosis (417 records, 6.6%), and endocrine, metabolic, or nutritional diseases (250 records, 3.9%) (Supplementary Table S2). Because this disease-spectrum summary was based on diagnosis records, the category counts should not be interpreted as unique-patient counts. Overall, the controls were not healthy population controls but hospitalized patients without liver cancer, many of whom had underlying liver disease or liver-cancer-related risk backgrounds.

Within the training set, patients in the liver cancer group had a median age of 57.0 years, which was significantly higher than that in the non-liver cancer group (55.0 years). The proportion of male patients was significantly higher in the liver cancer group than in the non-liver cancer group (P < 0.01). Regarding laboratory indicators, the liver cancer group exhibited significantly higher neutrophil counts, platelet counts, ALT, AST, GGT, ALP, and AFP and significantly lower lymphocyte percentage, albumin, albumin/globulin ratio, prothrombin time, and INR than the non-liver cancer group (all P < 0.01). The demographic and clinical characteristics of the patients with and without liver cancer were presented in Table 1.

Feature selection was performed using LASSO regression, with the minimum-error criterion (lambda.min; λ = 8.49 × 10⁻⁴), rather than lambda.1se, selected through 5-fold cross-validation (Figure 1; Supplementary Table S3; Supplementary Figure S2). At lambda.min, 57 nonzero model terms were retained. These included serum sialic acid (SSA), basophil percentage (BASO%), monocyte-to-lymphocyte ratio (MLR), alkaline phosphatase (ALP), retinol-binding protein (RBP), the aspartate aminotransferase-to-alanine aminotransferase ratio (AST/ALT), LDH, GGT, INR, PDW, FIB, ALT, PTA, IBIL, UA, BUN/Cr, P, NEUT%, CK-MB, TBA, and the additional terms listed in Supplementary Table S3. The maximum VIF was 6.77, and no retained term exceeded the conventional VIF >10 threshold (Supplementary Table S3; Supplementary Figure S1). In the bootstrap sensitivity analysis, 49 retained terms had a selection frequency of at least 0.80; the principal predictors highlighted in the interpretation analysis, including SSA, BASO%, ALP, RBP, and AST/ALT, were consistently retained. The six classification models were then constructed using the retained terms.

Model performance in the training and internal validation sets is summarized in Figures 4, 5. In the training set, XGBoost demonstrated the highest AUC, achieving 0.909 (95% CI, 0.903–0.915), followed by SVM and KNN. In the internal validation set, XGBoost also had the highest AUC (0.802; 95% CI, 0.788–0.817), while SVM ranked second with an AUC of 0.794. The validation-set AUCs for random forest, logistic regression, KNN, and Elastic Net were 0.718, 0.724, 0.747, and 0.725, respectively. At the prespecified probability threshold of 0.5, XGBoost achieved validation-set accuracy, sensitivity, specificity, and F1 score values of 0.730, 0.758, 0.698, and 0.753, respectively. Exploratory threshold analyses, including Youden-index and sensitivity-prioritized strategies, are presented in Supplementary Table S4; these candidate thresholds require prospective validation before clinical implementation. Pairwise DeLong testing in the validation cohort showed that XGBoost had a higher AUC than logistic regression, random forest, KNN, and Elastic Net (all P < 0.001), and a modestly higher AUC than SVM (P = 0.037; Benjamini-Hochberg-adjusted P = 0.046) (Supplementary Table S1).

DCA was used to evaluate the net clinical benefit of each model across risk thresholds (Figure 6). Across a wide range of threshold probabilities, XGBoost showed a higher net benefit than the other models, indicating a potentially more favorable trade-off between false-positive and false-negative classifications in this hospitalized cohort. SVM showed favorable performance in the high-risk threshold range (>0.6), whereas logistic regression, Elastic Net, and KNN yielded relatively low net benefits across most thresholds. Exploratory threshold analysis showed the expected trade-off: lowering the threshold increased sensitivity and NPV but reduced specificity and PPV. The optimal threshold should therefore be selected according to the intended clinical setting and validated prospectively.

The SHAP method and the built-in gain metric of XGBoost were used to evaluate feature importance for the optimal model. Features were ranked by mean absolute SHAP value (Figure 2). The five leading features were SSA, the AST/ALT ratio, ALP, BASO%, and MLR. The SHAP summary plot displays both the magnitude and direction of feature contributions at the individual-sample level. Higher values of SSA, the AST/ALT ratio, and ALP were mainly located on the positive side of the SHAP-value axis, indicating contributions toward a higher model-predicted probability of liver cancer. These patterns represent model-derived associations and should not be interpreted as causal biological effects. The gain metric, which quantifies the average improvement in tree-split performance attributed to a feature, ranked SSA, BASO%, MLR, ALP, and RBP as the five leading features (Figure 3); SSA ranked first by both importance measures.

Taken together, SSA, the AST/ALT ratio, ALP, BASO%, and MLR emerged as influential model predictors spanning liver function, inflammatory status, and metabolic or nutritional domains. Their importance supports their contribution to model discrimination within this hospitalized cohort but does not establish independent causal mechanisms.