The investigation focused on HCC patients who received curative hepatectomy at Guangxi Medical University Cancer Hospital from January 2014 through December 2021. Curative hepatectomy was defined as an R0 resection, with no microscopic tumor cells at the surgical margin, according to the Diagnosis and Treatment Guidelines for Primary Liver Cancer (24); The study’s criteria for participant selection were defined by specific inclusion and exclusion parameters. Inclusion criteria included: (1) patients were underwent R0 resection and all enrolled patients had adequate liver function reserve, as defined by an indocyanine green 15-minute retention rate (ICG R15) ≤30%; (2) Child–Pugh score of 5-7; (3) Eastern Cooperative Oncology Group performance status (ECOG PS) of 0 or 1; and (4) BCLC stage A or B. Exclusion criteria comprised: (1) history of other malignancies, (2) any preoperative anticancer treatment such as adjuvant chemotherapy, targeted therapy, immunotherapy, interventional therapy, or radiotherapy; (3) postoperative therapy, including aforementioned treatments; and (4) incomplete clinical data and follow-up duration of less than 2 months.

The Guangxi Medical University Ethics Committee (KY2025413) approved this study, which adhered to the Declaration of Helsinki principles.

Clinicopathological information of patients with HCC were collected, including (1) demographic information: gender, age, height, weight, etc. (2) laboratory parameters: total bilirubin, alpha-fetoprotein (AFP), albumin, platelets, etc. (3) liver disease-related information: Hepatitis B virus (HBV) infection status, HBV DNA level, etc. (4) tumor-related information: tumor number, tumor size, postoperative pathology, etc. The HCC stage was evaluated according to the BCLC staging classification system (25).

The main purpose was to assess overall survival (OS), tracked from the date of curative hepatectomy to either death from any cause or the last follow-up. The secondary endpoint was recurrence-free survival (RFS), defined as the time from curative hepatectomy until to the first occurrence of either disease recurrence or death from any cause.

Postoperative follow-ups were conducted at intervals of 1–2 months for the first year, followed by every 3 months thereafter until recurrence occurred. The follow-up programs included regular evaluations of liver function, AFP levels, and at least one contrast-enhanced imaging. HCC recurrence was diagnosed through a thorough evaluation of clinical history, AFP tests, and imaging results. The follow-up continued until 26 January 2025.

Continuous variables were expressed as means along with standard deviations (SD) and analyzed via Student’s t-test. Alternatively, they were presented as medians together with interquartile ranges (IQR) and analyzed using the Mann–Whitney U test. Categorical variables were presented as n (%) and compared with the Chi-square test. The Kaplan-Meier method was employed to generate the OS curves, and the log-rank test was utilized for their analysis. A multivariable Cox regression analysis model was developed to estimate the likelihood of hepatectomy risk, incorporating predictive factors identified through univariate analysis.

The dataset used for this analysis was complete, with no missing values for the analyzed variables. To develop and validate the predictive models, the entire dataset was randomly split into a training group (70%) and an internal validation group (30%) using the createDataPartition function from the caret package in R. This function employs a stratified random sampling strategy based on the OS status to ensure an equal distribution of deaths between the training and validation sets, thereby improving the robustness of the model evaluation. The random seed was set to 123 to ensure the complete reproducibility of the data partitioning.

Univariate Cox regression analyses were first performed to identify potential prognostic factors. To control the false discovery rate resulting from multiple testing, the p-values from the univariate analysis were further adjusted using the False Discovery Rate (FDR) correction via the Benjamini-Hochberg method. Variables with an FDR-adjusted p-value (P_FDR) < 0.05 were considered statistically significant and selected for inclusion in the subsequent multivariate Cox regression analysis. The proportional hazards assumption for the final multivariate model was verified using Schoenfeld residual tests, and no significant violations were found (global test p >0.05). Multicollinearity among the covariates in the multivariate Cox regression was assessed using the variance inflation factor (VIF). All VIF values were well below the threshold of 5 (BCLC: 1.23, MVI: 1.15, Size: 1.08), indicating no severe multicollinearity that would adversely affect the model estimates.

The independent risk factors identified from the multivariate Cox regression (BCLC stage, MVI, and tumor size) were used as input features for constructing eight ML models. These models include least absolute shrinkage and selection operator regression (Lasso_Cox), gradient boosting machine (GBM), random survival forests (RSF), boosting for Cox’s proportional hazards model (Coxboost), survival support vector machine (Survivalsvm), extreme gradient boosting (xgboost), super-predictor Cox model (superpc), and partial least squares with Cox’s proportional hazards model (plsRcox). The details of the algorithms and their hyperparameters are summarized in Supplementary Table S1 and Supplementary Figure S1 .

Hyperparameter tuning is a critical step to optimize model performance and prevent overfitting. For all models, hyperparameter tuning was conducted exclusively on the training group using resampling methods to avoid any information leakage. The specific tuning strategy, search spaces for key hyperparameters, and the criterion for evaluating model performance are described in detail for each algorithm in Supplementary Table S1 and Supplementary Figure S1 . We employed a systematic approach based on K-fold cross-validation (with K = 10) for hyperparameter exploration. The internal validation for hyperparameter tuning was performed using stratified 10-fold cross-validation (CV) on the training set. The stratification was based on the OS status to maintain the proportion of events (deaths) consistent across all folds. The performance of each hyperparameter combination was evaluated using the concordance index (C-index). The model configurations identified through this CV process were then evaluated on the independent internal validation set. The final optimal hyperparameter configuration for each algorithm was selected based on the highest average C-index across the 10 stratified CV folds. Critically, the independent internal validation set (30% of the data, held out from all tuning processes) was used only for post-selection evaluation of generalizability to unseen data—this strict separation ensures no information leakage into model selection. The mean and standard deviation of the C-index from both the 10-fold CV process and the independent internal validation for each final model are reported in Supplementary Table S2 .

The final model for each algorithm, with its hyperparameters fixed to the optimized values, was then refit on the entire training group and subsequently applied to the held-out validation group for an unbiased assessment of its performance. The performance of the ML models was comprehensively assessed using multiple metrics: the C-index, the area under the receiver operating characteristic curve (AUC), calibration curves, decision curve analysis (DCA), the Integrated Brier Score (IBS), and the Net Reclassification Index (NRI). In this study, a Cox proportional hazards model containing no predictor variables (the Null Model) was selected as the reference for NRI calculation. This model represents the average risk of the entire study cohort. This setup allows us to evaluate the absolute incremental value of all ML models over a “no-information” baseline. The NRI was calculated on the training set at 1-, 3-, and 5-year post-hepatectomy. The risk stratification threshold for all models was set at the median of their predicted risk probabilities.

The GBM model was also further compared with previously reported predictive models, including the BCLC staging system, metroticket Cox regression, Tumor-burden score (TBS), ERASL-pre score, and ERASL-post score, using similar evaluation metrics. These comparator models were applied based on their original published algorithms and were not re-trained on our dataset. Patients diagnosed with LHCC were divided into high-risk and low-risk groups based on the median risk scores obtained from the ML models.

The shapviz package was used to visualize contributions in the GBM model. SHAP summary plots showed feature impacts on predictions. Higher SHAP values in the GBM model meant a higher death likelihood. The SHAP feature importance plot orders features based on their average absolute SHAP values. SHAP force plots used colors (orange for positive, dark red for negative) to denote feature contributions.

For two-tailed tests, statistical significance was defined as a p-value of less than 0.05. All the statistical analyses were carried out using R version 4.4.2 (http://www.r-project.org/).

The LHCC group’s OS was much shorter than that of the SHCC group. (Hazard Ratio, HR: 1.81, 95% Confidence Interval, CI: 1.585-2.068, Figure 1A ). A total of 1,017 people died in the LHCC group, resulting in an overall mortality rate of 66.6%. In the LHCC group, the survival rates were 81.5% at the 1-year mark, 60.5% at the 3-year mark, and 52.2% at the 5-year mark. Similarly, the RFS of the LHCC group was significantly shorter compared to the SHCC group (HR: 1.526, 95% CI: 1.376-1.693, Figure 1B ). A total of 1,035 people experienced recurrence in the LHCC group, leading to an overall recurrence rate of 67.8%. Among the patients in the LHCC group, the rates of RFS were 47.6% after 1 year, 32.3% after 3 years, and 30.0% after 5 years.

Between 2014 and 2021, 2,565 HCC patients were enrolled, including 1,526 patients with LHCC and 1,039 patients with SHCC. For additional analysis, LHCC patients were randomly split into the training group (n = 1,069) and the validation group (n = 457) in a 7:3 ratio ( Figure 2 ). The baseline characteristics for both groups were presented ( Table 1 ). Notably, 41.3% of patients in the training group were classified as BCLC stage B, whereas 43.5% of the validation group fell into the same category. MVI was present in 52.5% of the training group and 51.0% of the validation group. In the training group, the average tumor size was 9.20 cm, while in the validation group, it was 9.05 cm. No statistically significant differences were observed between the two groups (p>0.05). The median OS follow-up times were 42.2 months for the training and 41.5 months for the validation group.

Univariate Cox regression analysis identified significant risk factors for OS in HCC patients, including the BCLC stage (p_FDR<0.001), number of tumors (p_FDR <0.001), MVI (p_FDR <0.001), and tumor size (p_FDR <0.001). After FDR adjustment, variables including BCLC stage, number of tumors, MVI, and tumor size remained significant (P_FDR <0.05). The multivariate Cox regression analysis, which satisfied the proportional hazards assumption (Schoenfeld global test p=0.439, Supplementary Table S3 ), identified BCLC stage (HR: 1.87, 95% CI: 1.48-2.36, p<0.001), MVI (HR: 1.55, 95% CI: 1.29-1.88, p < 0.001), and tumor size (HR: 1.04, 95% CI: 1.01-1.07, p = 0.005) were associated with an independent risk factors for OS in patients with LHCC ( Supplementary Table S4 ). No significant multicollinearity was detected among these variables in the multivariate model (all VIFs <5).

Among the eight ML models, the GBM model achieved AUC of 0.738 in the training group and 0.750 in the validation group, with corresponding C-index values of 0.715 and 0.737 ( Figures 3A-D ). The GBM model attained an IBS of 0.27 on the test set, indicating low overall prediction error. Compared to the null model, the GBM model demonstrated better NRI at 1-, 3-, and 5-year post-surgery, with NRI values of 32.84%, 32.74%, and 33.91%, respectively ( Supplementary Table S5 ). The GBM model attained AUC values of 0.738 (95% CI: 0.696-0.780), 0.708 (95% CI: 0.675-0.740), and 0.700 (95% CI: 0.668-0.732) for 1-, 3-, and 5-year OS, respectively ( Supplementary Figures S2A-C ). The validation groups showed AUC values for 1-, 3-, and 5-year OS of 0.742 (95% CI: 0.690-0.794), 0.744 (95% CI: 0.697-0.791), and 0.750 (95% CI: 0.706-0.795), respectively ( Supplementary Figures S2D-F ). Additionally, the calibration curves of the GBM model showed improved alignment between actual observations and model predictions for 1-, 3-, and 5-year OS in both the training ( Supplementary Figure S3 ) and validation groups ( Supplementary Figure S4 ). The DCA curves for 1-, 3-, and 5-year OS in both the training ( Supplementary Figures S5A-C ) and validation groups ( Supplementary Figures S5D-F ) suggested potential clinical utility of the GBM model, indicating a certain degree of positive benefit under the study’s evaluation framework. Patients with low GBM scores exhibited significantly better OS outcomes than those with high GBM scores ( Figures 3E, F ). Compared to other ML models, this was evident in both the training and validation groups ( Supplementary Figures S6 and S7 , respectively).

The SHAP summary plot for the GBM model illustrated the influence of individual features on the predictive outcomes ( Figure 4A ). The most influential features, ranked in descending order, were BCLC stage, tumor size, and MVI ( Figure 4B ). Furthermore, SHAP force plots ( Figures 4C, D ) were applied to explain the individual predictions. For example, for patient A, with MVI, BCLC stage B, and a 6 cm tumor, the GBM model predicted a risk score of “-0.911”, so he was classified into the high-risk group. For patient B, with MVI, BCLC stage A and a 6 cm tumor, the GBM model predicted a risk score of “-1.6”, so he was classified into the low-risk group.

A user-friendly website ( Figure 4E , https://ytx000.shinyapps.io/GBM-Shinyapp/) has been developed to facilitate the application of the GBM model in clinical practice. Practitioners can easily calculate individualized predicted risk scores for patients with LHCC by entering each patient’s clinical data into an online web server. Clinicians can input three key clinical parameters for LHCC patients—BCLC stage, MVI status, and tumor size—to generate individualized predictions of 1-, 3-, and 5-year OS rates. To illustrate its functionality, we present a representative case: a patient with BCLC stage A and absence of MVI, but with a large tumor size of 9.1 cm, received a GBM risk score of -1.33. As this score is higher than the mean risk score threshold of -1.34 used in our study, the model classified this patient into the high-risk group. The corresponding predicted 1-, 3-, and 5-year OS rates for this individual were 76.70%, 45.12%, and 26.54%, respectively. This example underscores the model’s ability to identify high-risk patients even among those with otherwise favorable clinical features (early BCLC stage and no MVI), highlighting the critical prognostic weight of tumor size captured by the GBM algorithm.

We compared the performance of the GBM model with previous postoperative predictive models ( Supplementary Table S6 ). The GBM model achieved AUC values of 0.714 (95% CI: 0.679-0.749), 0.708 (95% CI: 0.679-0.738), and 0.705 (95% CI: 0.674-0.736) for 1-, 3-, and 5-year OS, respectively ( Figures 5A-D ). The GBM model reached C-index values of 0.680, 0.663, and 0.656 that correspond to 1-, 3-, and 5-year OS, respectively ( Supplementary Figure S8A ). The DCA curves ( Supplementary Figures S8B-D ) for 1-, 3-, and 5-year OS showed the GBM model’s strong clinical utility and superior net benefit. The calibration curves ( Supplementary Figure S9 ) for the GBM model showed improved alignment between predicted probabilities and observed outcomes for 1-, 3-, and 5-year OS.

Notably, the GBM model demonstrated differential predictive performance between patients with high (NLR ≥5, n=156) and low (NLR <5, n=1,370) neutrophil-lymphocyte ratios, using a cutoff value supported by prior literature (26–28). The GBM model achieved a C-index of 0.819 in the high NLR group and 0.718 in the low NLR group, with corresponding AUC values of 0.814 and 0.732 ( Figures 6A-D ). The GBM model attained AUC values of 0.814 (95% CI: 0.729-0.899), 0.762 (95% CI: 0.685-0.840), and 0.774 (95% CI: 0.699-0.849) for 1-, 3-, and 5-year OS in high NLR groups, respectively ( Figure 6E ).The low NLR groups showed AUC values for 1-, 3-, and 5-year OS of 0.732 (95% CI: 0.697-0.768), 0.712 (95% CI: 0.683-0.740), and 0.708 (95% CI: 0.680-0.736), respectively ( Figure 6F ). Patients with low GBM scores exhibited significantly better OS outcomes than those with high GBM scores in both the high NLR group (log-rank p <0.01; Figure 6G ) and the low NLR group (log-rank p <0.01; Figure 6H ).