A time‑dependent risk prediction model for distant metastasis in early‑stage breast cancer based on explainable ensemble learning: a retrospective cohort study

Introduction

Breast cancer is the most common malignancy among women worldwide, with its incidence continuing to rise, posing a serious threat to women’s health (1,2). Although the overall survival (OS) rate for earlystage breast cancer patients has significantly improved with the advancement of screening and treatment modalities, distant metastasis remains the leading cause of treatment failure and patient death (3). Studies indicate that approximately 20–30% of earlystage breast cancer patients will eventually develop distant recurrence or metastasis (4).

Therefore, accurately identifying earlystage breast cancer patients at high risk of distant metastasis at the time of initial diagnosis is crucial for implementing intensive, personalized adjuvant treatment regimens (5). Conversely, for low‑risk patients, it may be possible to avoid unnecessary overtreatment and its associated toxicities. Currently, clinical practice relies heavily on the American Joint Committee on Cancer (AJCC) Tumor‑Node‑Metastasis (TNM) staging system and immunohistochemistry (IHC)‑based molecular subtyping [e.g., luminal A, luminal B, human epidermal growth factor receptor 2 (HER2)‑positive, triplenegative] for prognosis and treatment decisions (6). However, these traditional tools have significant limitations. The TNM staging system is primarily based on anatomical information and fails to fully incorporate the biological heterogeneity of tumors. While molecular subtyping represents a major advancement, substantial prognostic variation still exists within the same subtype (7).

With the development of high‑throughput sequencing and medical informatics, multigene expression profile assays have been used to assess recurrence risk and chemotherapy benefit. These tools have improved prediction accuracy, but their high cost and platform dependency limit global applicability (8,9). Furthermore, these genomic tools are primarily based on static tumor samples, making it difficult to integrate multi‑timepoint clinical information.

In recent years, machine learning and artificial intelligence (AI) have provided new solutions (10,11). Ensemble learning methods, such as stacking, often yield more stable and accurate predictive performance than single models (12,13). Simultaneously, there is clinical wariness towards “blackbox” models. SHapley Additive exPlanations (SHAP) can present model predictions in an intuitive, quantitative manner, clarifying each feature’s contribution to a patient’s predicted outcome (14,15). However, several challenges remain, including data heterogeneity, model interpretability, and the need for larger validation datasets (1).

Importantly, the goal of this study is not to discover novel genetic or protein biomarkers—an endeavor already accomplished by large‑scale genomic studies. Instead, our contribution lies in the methodological and translational domain: (I) we demonstrate that a stacking ensemble learning framework can extract additional prognostic signals from widely available, low‑cost clinicopathological variables beyond what traditional TNM‑ or Cox‑based models can achieve; (II) we provide, for the first time in this context, a systematic quantification of non‑linear threshold effects (e.g., the Ki‑67 escalation around 30%) and feature interactions; and (III) we translate these findings into an individualized, SHAP‑powered visual explanation tool that directly answers a clinician’s question: “Why does this specific patient have a high predicted risk?”. These features are absent in existing clinical tools (TNM, molecular subtypes, or even genomic signatures), and they address the growing demand for transparent, actionable AI in oncology. Recent studies have also highlighted the value of timedependent modeling (16) and competing risk approaches (17) in breast cancer prognostication, which we incorporate into our framework.

Based on a detailed clinicopathological database of 351 early‑stage breast cancer patients, this study aims to address the following questions: (I) can routinely available clinicopathological features be used to build a distant metastasis risk prediction model that surpasses traditional staging and subtyping? (II) Can an ensemble learning strategy further improve accuracy and robustness? (III) How can SHAP‑based interpretability make model predictions transparent and actionable? To this end, we propose and validate a time‑dependent risk prediction model based on a stacking ensemble framework and SHAP analysis. We present this article in accordance with the TRIPOD reporting checklist (available at https://gs.amegroups.com/article/view/10.21037/gs-2026-0163/rc).

Methods

Study design and patient cohort

This single-center, retrospective cohort study utilized electronic medical record data from patients diagnosed with breast cancer at Tangshan People’s Hospital between January 2016 and August 2024. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. This study was approved by the Institutional Review Board of Tangshan People’s Hospital (approval No. RMYY-LLKS-2026057). The requirement for informed consent was waived due to the retrospective nature of the study.

Inclusion criteria: (I) pathologically confirmed primary invasive breast carcinoma; (II) clinical stage I, II, or III at initial diagnosis (according to the AJCC 8th edition staging criteria); (III) underwent surgical treatment (mastectomy or breast-conserving surgery) as part of the initial treatment; (IV) had complete baseline clinicopathological data and follow‑up information.

Exclusion criteria: (I) presence of distant metastasis (stage IV) at initial diagnosis; (II) history of other malignancies; (III) missing follow-up data or loss to follow-up; (IV) received neoadjuvant therapy (as this alters the original pathological characteristics of the tumor).

Ultimately, 351 patients met the criteria and were included in the final analysis.

Data collection and variable definitions

Demographic, clinicopathological, treatment, and outcome data were extracted from electronic medical records. All pathological assessments were based on surgical specimens.

Clinicopathological variables:

• Demographics: age at diagnosis (analyzed both as a continuous variable and dichotomized at <35 years).
• Tumor characteristics: pathological TNM stage (pT, pN, per AJCC 8th edition), tumor size (maximum invasive diameter), histological type, and histological grade (Nottingham combined histologic grade, I–III). Lymphovascular invasion (LVI) was recorded as present or absent.

Biomarkers and molecular subtyping:

• ER/PR status: determined by IHC, with positivity defined as nuclear staining in ≥1% of tumor cells.
• HER2 status: defined per American Society of Clinical Oncology/College of American Pathologists (ASCO/CAP) guidelines. Positive status was assigned to cases with IHC 3+ or IHC 2+ with confirmation of gene amplification by in situ hybridization (ISH). Cases with IHC 0/1+ or IHC 2+ without ISH amplification were considered negative.
• Ki-67 proliferation index: recorded as the percentage of positively stained tumor cells. For analysis, it was treated both as a continuous variable and dichotomized using a clinically relevant cut-off of 30% (high vs. low).
• Molecular subtype: patients were categorized into intrinsic‑like subtypes based on biomarker profiles: luminal A (ER+ and/or PR+, HER2−, Ki-67 <30%), luminal B (HER2-negative) (ER+ and/or PR+, HER2−, Ki-67 ≥30%), luminal B (HER2-positive) (ER+ and/or PR+, HER2+), HER2-enriched (ER−, PR−, HER2+), and triple-negative (ER−, PR−, HER2−).
• Treatment data: information on surgery type, adjuvant chemotherapy, radiotherapy, and endocrine therapy was collected. To align with the clinical goal of pre-treatment risk prediction, treatment variables were not used as predictors in the final ensemble model. They were, however, considered in sensitivity analyses and model interpretation to contextualize findings.

Outcomes and follow-up

The primary endpoint was distant metastasis-free survival (DMFS), defined as the time from the date of surgery to the first radiologic [e.g., computed tomography (CT), bone scan, positron emission tomography (PET)-CT] or pathologic confirmation of distant metastasis. Patients without an event were censored at the date of their last disease-free assessment. Non-breast cancer deaths were treated as competing events. Therefore, we additionally performed competing risk analysis using the Fine-Gray subdistribution hazard model (18) to estimate the cumulative incidence of distant metastasis, accounting for death from other causes (17). OS was a secondary endpoint. Follow-up data were collected until November 2024, with patients typically undergoing clinical and radiologic assessments every 3 to 6 months post-surgery.

Data preprocessing, feature engineering and selection

To prepare the data for modeling, a sequential preprocessing and feature selection pipeline was implemented. Missing data were infrequent; for the Ki-67 index, missing in 12 cases (3.4%), multiple imputation was performed using a non-parametric random forest algorithm (missForest R package) to leverage the correlation structure of other variables, while records with missing categorical data (<1%) were excluded listwise. To capture additional clinical nuance, we engineered several features, including the lymph node ratio (number of positive nodes divided by total nodes examined) and pre-specified interaction terms. All continuous variables were then standardized to a mean of 0 and a standard deviation of 1 (Z-score normalization) to ensure comparability across different scales. To address the relative rarity of metastasis events, we applied inverse probability of censoring weighting (IPCW) and assigned higher misclassification costs to the event class during model training. Following this preprocessing, and exclusively on the training set to prevent data leakage, we performed feature selection using recursive feature elimination with cross-validation (RFECV). With a random survival forest (RSF) as the base estimator, this procedure iteratively eliminated the least important features and selected the optimal subset that maximized the cross-validated concordance index (C-index), thereby enhancing model parsimony and mitigating overfitting.

Prediction model construction

A two-layer stacking ensemble framework was adopted to predict DMFS. The entire cohort was randomly partitioned into a training set (70%, n=246) for model development and an independent internal validation set (30%, n=105) for performance evaluation. In the first layer, four diverse survival analysis models were trained on the training set using 5-fold cross-validation to generate robust out-of-fold predictions: (I) CoxNet, an elastic-net penalized Cox proportional hazards model; (II) RSF, a tree-based ensemble method capable of capturing non-linear relationships and interactions; (III) gradient boosting survival trees (GBST), which optimizes a Cox partial likelihood loss function; and (IV) DeepSurv, a deep neural network architecture designed to model complex, high-order patterns in survival data. Subsequently, the out-of-fold predicted risk scores from these base learners were extracted and used as input features for the second-layer meta-learner. These derived features, optionally concatenated with the originally selected feature subset, formed the training data for a LightGBM survival model, which was chosen as the meta-learner for its computational efficiency, predictive accuracy, and native support for ranking loss functions in survival settings. For comparison, each base learner was also trained and evaluated independently as a single model under the same training‑validation split. All models were implemented in Python (version 3.9) using libraries including scikit-survival, lifelines, pycox, and lightgbm, with hyperparameter tuning performed via grid search or Bayesian optimization to maximize the C-index on the training set. To prevent overfitting, we applied early stopping to the LightGBM meta-learner, set L2 regularization (lambda_l2 =0.1), and limited maximum leaves per tree to 16.

Crucially, the out-of-fold predictions for each base learner were generated using only the training set via 5-fold cross-validation; the validation set was never used in this process, ensuring no information leakage.

To benchmark the performance of our ensemble model against clinically familiar tools, we constructed two baseline models using the same training set: Model A included only TNM stage (pT + pN) and molecular subtype (standard of care); Model B was a standard Cox proportional hazards model incorporating all eight selected features without penalization. These models were evaluated on the same validation set using the same metrics.

Model evaluation and validation

Model performance was rigorously evaluated on the independent validation set. In the training set, we also computed the C-index to assess potential overfitting. All analyses were conducted in Python (version 3.9) using scikit-survival, lifelines, lightgbm, and shap, with R (version 4.2) employed for supplementary validation.

Discrimination was quantified using the C-index as a global measure of ranking accuracy, and time-dependent area under the curve (AUC) (16) at 1, 3, and 5 years to assess predictive performance at clinically relevant time points. Calibration was evaluated by plotting calibration curves, which compared the mean predicted survival probabilities against Kaplan-Meier (KM) observed estimates at the same time points, and by calculating the integrated Brier score (IBS) as an overall measure of prediction error across the entire follow-up period, with lower values indicating superior accuracy. To account for competing risks, we also plotted calibration curves for the cumulative incidence function (CIF) based on the Fine-Gray model, and reported the Fine-Gray model’s C-index as a sensitivity measure (17).

To obtain robust performance estimates, we performed 100 repeated random splits (70/30) and 500 bootstrap resamples, reporting mean C-index and AUC with 95% confidence intervals (CIs) derived from the percentile method. This addresses the concern of optimistic bias due to a single random split.

To assess clinical utility, decision curve analysis (DCA) was performed to quantify the net benefit of using the model to guide treatment decisions across a range of clinically reasonable threshold probabilities, relative to treat all or treat none strategies. Furthermore, patients were stratified into low-, intermediate-, and high-risk groups based on their final model risk scores using the K-means clustering algorithm. The number of clusters (k=3) was determined by the silhouette coefficient and elbow method, and also corresponded to clinically meaningful thresholds: low risk (predicted 5-year metastasis risk <10%), intermediate risk (10–30%), and high risk (>30%). Prognostic separation was validated by comparing DMFS across these groups via cumulative incidence curves (CIF) based on the Fine-Gray model and the log-rank test, with a two-sided P value <0.05 considered statistically significant.

Statistical analysis

Statistical analyses were performed using Python 3.9 and R 4.2. Continuous variables were standardized (Z-score). Missing Ki-67 data (3.4%) were imputed by random forest. Model performance was assessed using C-index, time-dependent AUC, IBS, and calibration curves. Internal validation used 100 repeated random splits and 500 bootstrap resamples to obtain 95% CIs. Competing risks were analyzed with the Fine-Gray model. Two-sided P<0.05 was considered significant.

Model interpretability analysis

To enhance the transparency and clinical interpretability of the final Stack-LightGBM ensemble model, we employed the SHAP framework, implemented via the shap Python library (14,15). For global interpretability, a SHAP summary plot was generated to visualize feature importance, displaying the mean absolute contribution of each variable to the model output across the entire cohort, ranked in descending order. SHAP dependence plots were constructed to illustrate the relationship between individual feature values and their corresponding SHAP values, thereby revealing potential non-linear effects. These plots were further enhanced by coloring points according to a second feature (e.g., molecular subtype) to uncover interaction effects between variables. For local interpretability, SHAP force plots and waterfall plots were generated for individual patients to decompose their predicted risk scores. These visualizations quantitatively demonstrated how each clinical characteristic contributed—either positively or negatively—to deviate from the baseline prediction, thereby providing a personalized explanation for why a specific patient was classified as high-risk. This individualized risk profiling offers clinically actionable insights at the patient level and directly answers the clinician’s question: “Why does this specific patient have a high predicted risk?” (19).

Results

Patient baseline characteristics

A total of 351 patients with early‑stage breast cancer were included in this study. The median age was 52 years (range, 24–85 years). The median follow‑up time was 42 months (range, 6–96 months), during which 89 distant metastasis events occurred. The 5‑year DMFS rate was 74.6%. The main baseline clinicopathological characteristics of the patients are detailed in Table 1. Luminal B (HER2‑negative) was the most prevalent subtype (40.2%), followed by luminal B (HER2‑positive, 21.1%) and triple-negative breast cancer (TNBC) (18.5%). Most patients (65.8%) were node‑negative (pN0), but 34.2% had axillary lymph node involvement.

Feature selection results

Through recursive feature elimination with cross‑validation, eight features with the highest predictive value for DMFS were selected from an initial pool of 20 candidate features, forming the core variable set for the final model: pN stage (number of positive lymph nodes); Ki‑67 index group (≥30% vs. <30%); molecular subtype (with luminal A as reference); pT stage; pathological grade (III vs. I/II); age group (<35 years vs. ≥35 years); tumor size (continuous); and LVI status. The effective events per variable (EPV) was 89/8≈11.1, which is borderline but within acceptable range; regularization and feature selection were applied to mitigate overfitting.

Comparison of predictive model performance

The predictive performance of individual models and the Stack‑LightGBM ensemble model was compared on the independent internal validation set (n=105) (Table 2, Figure 1). The Stack‑LightGBM ensemble model demonstrated the highest discrimination at all time points. In the training set, the Stack‑LightGBM achieved a C‑index of 0.85 (95% CI: 0.81–0.89). In the validation set, the Stack‑LightGBM model achieved a global C‑index of 0.82 (95% CI: 0.77–0.87). The difference of 0.03 between training and validation C‑indices indicates minimal overfitting. The time‑dependent AUC for predicting 1‑, 3‑, and 5‑year DMFS were 0.85 (95% CI: 0.80–0.90), 0.82 (95% CI: 0.77–0.87), and 0.79 (95% CI: 0.74–0.84), respectively, significantly outperforming the best single model (GBST, C‑index =0.78, 95% CI: 0.73–0.83). Calibration curves showed good agreement between the predicted survival probabilities by the Stack‑LightGBM model and the actual observed survival probabilities via KM estimation at the 1‑, 3‑, and 5‑year time points (Figure 2).

Figure 1 Comparison of time-dependent AUC for each model on the validation set. AUC, area under the curve; GBST, gradient boosting survival trees; RSF, random survival forest.

Figure 2 Calibration curves of the Stack-LightGBM ensemble model on the validation set. The diagonal line represents perfect calibration. The model’s predicted probabilities show good agreement with the actual KM-estimated probabilities at 1-year (A), 3-year (B), and 5-year (C) time points. KM, Kaplan-Meier.

The benchmark models showed inferior performance: Model A (TNM + molecular subtype) achieved a C‑index of 0.71 (95% CI: 0.66–0.76); Model B (Cox‑8f) achieved a C‑index of 0.76 (95% CI: 0.71–0.81). The Stack‑LightGBM significantly outperformed both (P<0.01, DeLong test). The detailed performance of the benchmark models is shown in Table 3.

As a sensitivity analysis to account for competing risks, the FineGray subdistribution hazard model (17,18) yielded a Cindex of 0.81 (95% CI: 0.76–0.86), similar to the main model, and its calibration curves for the CIF showed good agreement (Figure S1).

Model interpretability analysis (SHAP)

The SHAP summary plot (Figure 3A) revealed that the number of positive lymph nodes (pN) was the most important feature influencing DMFS risk, with the largest range of SHAP values, indicating its most significant and variable impact on the risk score. This was followed by high Ki‑67 expression and molecular subtype (particularly triple‑negative and HER2‑enriched subtypes). High pN stage, Ki‑67 ≥30%, and triple‑negative or HER2‑enriched subtypes were all associated with higher SHAP values (i.e., increased metastasis risk). Age <35 years also showed an association with higher risk.

Figure 3 Results of SHAP interpretability analysis. (A) SHAP summary plot showing global feature importance and the direction of each feature’s effect on the model output; (B) SHAP dependence plot showing the nonlinear association between the Ki-67 index and model output, with points colored by molecular subtype; (C) SHAP contribution plot for a representative high-risk patient (ID #247), showing feature-specific contributions to the model output. N, node; SHAP, SHapley Additive exPlanations; T, tumor; TNBC, triple-negative breast cancer.

SHAP dependence plots revealed non‑linear effects of Ki‑67: when Ki‑67 values were below 20%, their impact on risk was minimal and variable; however, when exceeding 30%, SHAP values increased sharply, demonstrating a clear pro‑risk effect (Figure 3B). Coloring points by molecular subtype revealed interactions: for example, within the triple‑negative subtype, even moderate Ki‑67 levels were associated with higher SHAP values, suggesting that the “harm” of proliferative activity is amplified in this aggressive subtype.

Figure 3C presents a SHAP contribution plot for an actual high-risk patient (ID #247, who developed liver metastasis 18 months post-surgery). The largest positive contributors to the model output were pN3 (11 positive lymph nodes; SHAP value +3.5), Ki-67 of 80% (+2.8), and the triple-negative subtype (+2.2), whereas the absence of lymphovascular invasion had a small negative contribution (−0.2). The baseline and final model outputs were 0.30 and 13.10, respectively. This visualization enables clinicians to understand the feature-specific composition of the patient’s high predicted risk.

Risk stratification and clinical validation

Based on the risk scores calculated by the Stack‑LightGBM model for all patients, K‑means clustering divided them into low‑, intermediate‑, and high‑risk groups (Figure 4). The number of clusters (k=3) was determined by the silhouette coefficient and elbow method, and also corresponded to clinically meaningful thresholds: low risk (predicted 5‑year metastasis risk <10%), intermediate risk (10–30%), and high risk (>30%).

Figure 4 Distribution of risk scores and K-means clustering (low, intermediate, high).

KM analysis showed 5-year DMFS rates of 95.2%, 78.5%, and 51.3% for the low-, intermediate-, and high-risk groups, respectively (log-rank P<0.001). To account for competing risks (non-breast cancer death), Fine-Gray subdistribution hazard analysis yielded 5-year cumulative incidence of distant metastasis of 4.8%, 21.5%, and 48.7%, respectively (Figure 5). These findings were consistent with the KM-based risk stratification and confirmed significant prognostic separation among the three groups. The high-risk group concentrated most patients with adverse features such as pN2–3, triple-negative/HER2‑positive, and high Ki-67.

Figure 5 Cumulative incidence curves of distant metastasis for the low-, intermediate-, and high-risk groups based on the Fine-Gray competing-risk analysis.

Clinical decision support and tool development

DCA (Figure 6) indicated that across a wide range of threshold probabilities from 5% to 50%, using the Stack‑LightGBM model to guide decisions on treatment intensification provided a net benefit higher than the extreme strategies of treat all and treat none. For example, at a decision point that might correspond clinically to intensify treatment for patients with a predicted 3‑year metastasis risk >20%, the model provided a clear net clinical benefit.

Figure 6 Decision curve analysis.

Discussion

This study successfully developed and validated an ensemble learning model for predicting distant metastasis risk in early‑stage breast cancer based on routine clinicopathological data. The core strengths of this model lie in the organic integration of high predictive performance, SHAP‑enhanced interpretability, and clinical utility.

Regarding predictive performance, our Stack‑LightGBM ensemble model (validation Cindex =0.82) significantly outperformed the traditional Cox model and advanced single machine learning models. This result aligns with the recent trend that ensemble learning can effectively integrate the strengths of different algorithms, enhancing robustness and accuracy in complex medical prediction tasks (12,13,20). We acknowledge that a direct performance comparison between our model and commercial genomic assays (e.g., Oncotype DX) is not methodologically appropriate, as the latter are based on entirely different data types (gene expression) and have been validated in large prospective trials (8). Therefore, we do not claim superiority over these well‑established tools. Instead, our model should be viewed as a complementary, costeffective alternative for settings where genomic testing is not accessible or affordable. Its novelty resides not in discovering new biological drivers, but in repurposing routine clinical variables with advanced ensemble learning and explainable AI to achieve a level of predictive granularity and transparency that approaches—but does not replace—genomic tools, while being immediately deployable worldwide (1).

What actionable knowledge does our model provide beyond existing clinical tools? Although the individual predictors (pN, Ki‑67, molecular subtype) are well‑established, traditional tools such as TNM staging or molecular subtyping do not offer:

• Quantified non‑linear thresholds—our SHAP analysis identifies a clear risk escalation point for Ki‑67 at approximately 30%, suggesting that patients exceeding this threshold may require closer surveillance even within the same molecular subtype.
• Interaction‑aware risk modification—for example, young age (<35 years) is a much stronger risk factor in triple‑negative breast cancer than in luminal subtypes, a nuance that standard subtyping alone misses.
• Individualized risk decomposition—the SHAP contribution plot (Figure 3C) breaks down a patient’s model output into feature-specific contributions (e.g., “pN3 contributed +3.5, Ki-67 contributed +2.8, and the triple-negative subtype contributed +2.2”), enabling targeted clinical discussions. None of these actionable insights are available from conventional staging or molecular subtyping, nor from genomic recurrence scores that output a single numeric risk estimate without featurewise explanation (21). Therefore, the novelty of our study lies not in discovering new biological markers, but in demonstrating how existing routine data can be re‑engineered with ensemble learning and SHAP to generate transparent, individualized, and actionable prognostic information—a contribution that is methodologically innovative and clinically relevant, especially in resource‑constrained healthcare systems.

However, SHAP is a post‑hoc attribution method and does not fully explain the internal decision‑making of complex ensembles. Therefore, we describe our approach as “SHAPenhanced interpretability” rather than a full explainable AI framework (15,19).

The risk stratification based on the model demonstrated excellent discriminatory ability. The high‑risk group (16.8% of patients) had a 5‑year cumulative incidence of distant metastasis of 48.7%, while the low‑risk group (40.5%) had a cumulative incidence of 4.8% (Fine‑Gray model). This polarized risk distribution is advantageous for clinical decision‑making. For the low‑risk group, our results provide strong support for considering treatment de‑escalation (e.g., shortening endocrine therapy duration, omitting certain chemotherapy regimens), potentially reducing treatmentrelated toxicity and healthcare costs (22). For the high‑risk group, it reinforces the necessity of intensive adjuvant therapy and close follow‑up. DCA further confirmed, from a clinical utility perspective, that applying this model at reasonable risk thresholds can yield positive net benefit.

There are several limitations in this study. First, as a single‑center retrospective study, potential selection and information biases exist. Although rigorous internal validation was performed, the model still requires external validation in multi‑center, prospective cohorts with different demographic characteristics and treatment protocols to demonstrate its generalizability. We are currently conducting a multicenter external validation (n=220, ethics approved). Second, our data were primarily based on baseline characteristics at surgery. Integrating dynamic information during and after treatment (e.g., pathological response after neoadjuvant therapy, adherence to adjuvant endocrine therapy, circulating tumor DNA changes) could enable the construction of more refined dynamic prediction models (23). Third, this study did not incorporate radiomics or digital pathology features. The fusion of such multimodal data with clinical information is a direction for future improvement in prediction accuracy (24,25). Fourth, the effective EPV was approximately 11, which is borderline; however, we applied regularization and feature selection to mitigate overfitting. Although an EPV of 10–15 is considered acceptable for logistic regression, survival models may require more events; our use of penalized methods partially addresses this concern. Fifth, our use of KM for DMFS may overestimate survival; we have now repeated analyses using FineGray competing risk models (17,18). which yielded similar risk stratification (see “Supplementary Results” in Appendix 1). Sixth, we did not identify any new biological markers; all predictors used are already well‑known in clinical practice. This is by design—our aim is to maximize the utility of existing data, not to discover novel drivers. Readers seeking biological discovery should refer to dedicated genomic studies. Finally, while SHAP greatly enhances interpretability, the translation of model predictions into final clinical decisions still requires clinicians to integrate the patient’s overall condition, values, and preferences.

Conclusions

This study developed and preliminatively validated an SHAP‑enhanced interpretable ensemble learning model based on multi‑dimensional routine clinicopathological features for predicting the risk of distant metastasis in early‑stage breast cancer patients. The model (Stack‑LightGBM) surpassed traditional methods in predictive accuracy and achieved transparency in prediction logic through SHAP analysis, enabling the identification of key risk factors, quantification of their non‑linear threshold effects and feature interactions, and providing individualized risk contribution interpretations. Risk stratification based on this model effectively discriminated patient subgroups with markedly different prognoses. This study provides a proof‑of‑concept for integrating high‑performance, interpretable machine learning tools into the clinical decision‑making workflow for breast cancer, particularly in settings where genomic testing is not available or affordable. Future work will focus on multi‑center external validation, prospective evaluation of clinical effectiveness, and integration with other modal data, with the ultimate goal of achieving truly personalized and precise breast cancer management.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://gs.amegroups.com/article/view/10.21037/gs-2026-0163/rc

Data Sharing Statement: Available at https://gs.amegroups.com/article/view/10.21037/gs-2026-0163/dss

Peer Review File: Available at https://gs.amegroups.com/article/view/10.21037/gs-2026-0163/prf

Funding: This work was supported by a research project entitled Application of SPECT/CT Image Fusion Technology in the Localization of Sentinel Lymph Nodes in Breast Cancer (project No. 12130229b).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://gs.amegroups.com/article/view/10.21037/gs-2026-0163/coif). H.W. and X.X. report support from a research project entitled Application of SPECT/CT Image Fusion Technology in the Localization of Sentinel Lymph Nodes in Breast Cancer. The other authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. This study was approved by the Institutional Review Board of Tangshan People’s Hospital (approval No. RMYY-LLKS-2026057). The requirement for informed consent was waived due to the retrospective nature of the study.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Abbas GH, Khouri ER, Thaher O, et al. Predictive modeling for metastasis in oncology: current methods and future directions. Ann Med Surg (Lond) 2025;87:3489-508. [Crossref] [PubMed]
2. Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin 2021;71:209-49. [Crossref] [PubMed]
3. Harbeck N, Penault-Llorca F, Cortes J, et al. Breast cancer. Nat Rev Dis Primers 2019;5:66. [Crossref] [PubMed]
4. Noman SM, Fadel YM, Henedak MT, et al. Leveraging survival analysis and machine learning for accurate prediction of breast cancer recurrence and metastasis. Sci Rep 2025;15:3728. [Crossref] [PubMed]
5. Cardoso F, Kyriakides S, Ohno S, et al. Early breast cancer: ESMO Clinical Practice Guidelines for diagnosis, treatment and follow-up†. Ann Oncol 2019;30:1194-220. [Crossref] [PubMed]
6. Giuliano AE, Edge SB, Hortobagyi GN. Eighth Edition of the AJCC Cancer Staging Manual: Breast Cancer. Ann Surg Oncol 2018;25:1783-5.
7. Zarean Shahraki S, Azizmohammad Looha M, Mohammadi Kazaj P, et al. Time-related survival prediction in molecular subtypes of breast cancer using time-to-event deep-learning-based models. Front Oncol 2023;13:1147604. [Crossref] [PubMed]
8. Pilgram L, Yang K, Beltran-Bless AA, et al. Transfer Learning and Machine Learning for Training Five-Year Survival Prognostic Models in Early Breast Cancer: Development and Validation Study. J Med Internet Res 2026;28:e88665. [Crossref] [PubMed]
9. Clift AK, Dodwell D, Lord S, et al. Development and internal-external validation of statistical and machine learning models for breast cancer prognostication: cohort study. BMJ 2023;381:e073800. [Crossref] [PubMed]
10. Esteva A, Robicquet A, Ramsundar B, et al. A guide to deep learning in healthcare. Nat Med 2019;25:24-9. [Crossref] [PubMed]
11. Deo RC. Machine Learning in Medicine. Circulation 2015;132:1920-30. [Crossref] [PubMed]
12. Gurcan F. Enhancing breast cancer prediction through stacking ensemble and deep learning integration. PeerJ Comput Sci 2025;11:e2461. [Crossref] [PubMed]
13. Li X, Gao M, Zhang C, et al. A robust stacked neural network approach for early and accurate breast cancer diagnosis. Front Med (Lausanne) 2025;12:1644857. [Crossref] [PubMed]
14. Lundberg SM, Lee SI. A Unified Approach to Interpreting Model Predictions. In: Advances in Neural Information Processing Systems; 2017:4765-74.
15. Zaheer Sajid M, Fareed Hamid M, Qureshi I. Explainable and uncertainty-aware ensemble framework with causal analysis for breast cancer detection. Front Oncol 2025;15:1751090. [Crossref] [PubMed]
16. Park E, Park JH, Lee HK, et al. Pattern of Metastasis as a Risk Factor for Brain Metastasis in Patients With Breast Cancer: A Time-Dependent Cox Regression. Asia Pac J Clin Oncol 2026; Epub ahead of print. [Crossref]
17. Mariotto AB, Botta L, Bernasconi A, et al. Prediction of Risk of Metastatic Recurrence for Female Breast Cancer Patients in the Presence of Competing Causes of Death. Cancer Epidemiol Biomarkers Prev 2023;32:1683-9. [Crossref] [PubMed]
18. Fine JP, Gray RJ. A proportional hazards model for the subdistribution of a competing risk. J Am Stat Assoc 1999;94:496-509.
19. Adnan N, Zand M, Huang THM, et al. Construction and Evaluation of Robust Interpretation Models for Breast Cancer Metastasis Prediction. IEEE/ACM Trans Comput Biol Bioinform 2022;19:1344-53. [Crossref] [PubMed]
20. Buyrukoğlu G. Survival analysis in breast cancer: evaluating ensemble learning techniques for prediction. PeerJ Comput Sci 2024;10:e2147. [Crossref] [PubMed]
21. Guan Z, Huang T, McCarthy AM, et al. Combining Breast Cancer Risk Prediction Models. Cancers (Basel) 2023;15:1090. [Crossref] [PubMed]
22. Andre F, Ismaila N, Allison KH, et al. Biomarkers for Adjuvant Endocrine and Chemotherapy in Early-Stage Breast Cancer: ASCO Guideline Update. J Clin Oncol 2022;40:1816-37. [Crossref] [PubMed]
23. Pascual J, Attard G, Bidard FC, et al. ESMO recommendations on the use of circulating tumour DNA assays for patients with cancer: a report from the ESMO Precision Medicine Working Group. Ann Oncol 2022;33:750-68. [Crossref] [PubMed]
24. Huang Y, Wang X, Cao Y, et al. Multiparametric MRI model to predict molecular subtypes of breast cancer using Shapley additive explanations interpretability analysis. Diagn Interv Imaging 2024;105:191-205. [Crossref] [PubMed]
25. Chen Y, Chen S, Tang W, et al. Multiparametric MRI Radiomics With Machine Learning for Differentiating HER2-Zero, -Low, and -Positive Breast Cancer: Model Development, Testing, and Interpretability Analysis. AJR Am J Roentgenol 2025;224:e2431717. [Crossref] [PubMed]