This study adhered to the Transparent Reporting of a multivariable prediction model for Individual Prognosis Or Diagnosis (TRIPOD + AI) checklist (25). The flowchart illustrates our study design ( Figure 1 ). This study encompassed three patient cohorts. The SEER database, a publicly accessible repository, is established and maintained by the National Cancer Institute. Initially, data from the SEER 17 registries research data [(2000–2020); version 8.4.3] were utilized. Given that Human epidermal growth factor receptor 2 (HER2) status and specific metastatic organ sites were incorporated from 2010 onwards, the inclusion criteria were as follows: (1) female sex; (2) diagnosis year falling within 2010 to 2020; (3) anatomical site and morphological coding aligned with the International Classification of Diseases Oncology Third Edition (ICD-O-3); (4) concurrent liver metastases at breast cancer diagnosis; and (5) survival time >0 months. Patients with multiple primary tumors were excluded. Furthermore, retrospective data were gathered from BCLM patients treated at Jiangmen Central Hospital (JCH) and Cancer Hospital of Shantou University Medical College (CHSU) between January 2010 and June 2023. This study was approved by the Ethics Committees of JCH (No. 2024312) and CHSU (No. 2024216). Our Ethics Committees, in accordance with the Guidelines for Ethical Review Committees of Clinical Research Involving Humans (2023 edition), determined that informed consent was not required and thus approved the waiver of informed consent.

The following patient characteristics were obtained: age, race, marital status, median household income (inflation adjusted), histological type, tumor location, histologic grade, molecular subtype, T stage, N stage, surgery, radiotherapy, chemotherapy, bone metastases, brain metastases, and lung metastases. Tumor-related variables comprised histological type, coded according to the ICD-O-3; primary tumor location, based on ICD-O-3 topography codes; and histological grade. Subtype was determined using available molecular or clinicopathological surrogates (hormone receptor and HER2 status). Tumor stage was assessed using the American Joint Committee on Cancer (AJCC) TNM staging system (8th edition). The primary endpoint was OS and the secondary endpoint was BCSS. The median follow-up time was 56 months (IQR: 28.0–86.0) for patients from the SEER database and 35 months (20.0–58.0) for patients from two hospitals in China. Survival rates and event counts for each prediction interval are presented in Supplementary Table S1 .

In order to reduce redundant variables and model overfitting, we employed univariate and multivariate Cox regression analyses within each training fold to select independent factors associated with prognosis in the training group. Statistically significant features were utilized for subsequent ML model construction. Five commonly used ML algorithms, including Random Forest (RF), Logistic Regression (LR), XGBoost, Decision Tree (DT), and Gradient Boosting Decision Tree (GBDT), were applied to predict 6-month, 1-year, 3-year, and 5-year Overall Survival (OS) and Breast Cancer-Specific Survival (BCSS). For each prediction horizon t ∈{0.5,1,3,5} year, we defined a binary outcome Y t = 1 if the event occurred by t  and Y t = 0 otherwise. Patients censored prior to t had unknown labels and were handled using inverse probability-of-censoring weighting (IPCW). The censoring survival function G ^ ( t ) was estimated via the Kaplan–Meier method, with right-censoring defined as alive at last follow-up. Each labeled observation received weight:

, i.e., at the time the label became known. Classifiers were then trained with these weights and evaluated using IPCW-weighted Area Under the Curve (AUC) and Brier score at each prediction horizon (26). RF is a non-parametric machine learning algorithm based on ensemble learning, which predicts by constructing multiple decision trees and integrating them through voting or averaging. Each decision tree is generated by bootstrapping the training set, with only a random subset of features considered for splitting at each node. This randomness effectively reduces model variance and enhances generalization performance. LR is a generalized linear model that estimates probabilities by mapping the output of a linear function to a sigmoid function and estimates model parameters by maximizing the likelihood function. XGBoost is an efficient gradient boosting decision tree algorithm that iteratively trains decision tree models by optimizing the negative gradient of the loss function and utilizes regularization techniques to control model complexity. DT is a classification and regression method based on tree structure, which partitions the dataset into different categories or values through a series of decision nodes. GBDT is an algorithm that enhances model performance by iteratively training decision trees and optimizing the loss function. It constructs decision trees sequentially and trains each subsequent tree using the residual of the previous tree, effectively reducing model bias.

To enhance model robustness, hyperparameters were optimized via grid search with 10-fold cross-validation. Patients from the SEER database were randomly divided into training and internal testing groups at an 8:2 ratio. Patients from two Chinese hospitals were used as an external testing group to further validate the generalizability of the optimal model. The performance of ML models was evaluated using AUC, accuracy, and F1 score of the training and testing groups. Calibration curves were used to assess the accuracy and reliability of model predictions. Decision Curve Analysis (DCA) was employed to determine the clinical utility of the models. Confusion matrices intuitively displayed the classification performance of the models. SHAPley Additive exPlanations (SHAP) values were computed using the “shap” library to explain the contribution or importance of each feature to the model.

To minimize immortal-time bias and avoid exposure misclassification around surgery, PTS was modeled as a time-dependent covariate. For patients undergoing surgery, follow-up time before surgery was classified as “non-PTS,” and time after surgery as “PTS.” Patients with a surgical history but no recorded surgery date were excluded from the main analysis, as their exposure periods could not be reliably defined. We performed multivariable analyses using a time-dependent Cox regression model based on PTS, incorporating an interaction term with log(time) to examine time-varying effects, and generated hazard ratio (HR) curves over time. For descriptive visualization, Simon–Makuch curves were additionally plotted. A segmented Cox model was applied with the follow-up time divided at 24 months, and a 2-month landmark analysis was conducted among patients who survived beyond this time point. To evaluate robustness against unmeasured confounding, the E-value was calculated using the following formula:

Categorical data were presented as counts and percentages and were compared using either the Chi-squared test or Fisher’s exact test depending on the sample size. Kaplan-Meier survival analysis was conducted using the log-rank test. Univariate and multivariate Cox regression analyses were employed to identify modeling features. The generalized variance inflation factor (GVIF) was utilized to detect multicollinearity in the Cox regression model. An adjusted GVIF value less than 2 was considered acceptable. Statistical analyses were performed using R software version 4.2.1 (r-project.org/) or Python version 3.8 (Python Software Foundation). A significance level of P< 0.05 indicated statistical significance.

Finally, we collected data on 4817 BCLM patients who met the eligibility criteria from the SEER database. As shown in Table 1 , 1545 patients (32.07%) were aged 50 years or younger, 1945 patients (40.38%) were between the ages of 51 and 65, and 1327 patients (27.55%) were 66 years or older. The majority of patients were Caucasian (72.33%). Approximately 45.84% of the patients were married, and 24.60% were single or homosexual. In terms of household income, 54.35% of the patients had incomes exceeding $70,000. About 80% of the patients were diagnosed with invasive ductal carcinoma (IDC), and 25.76% had tumors located in the upper outer quadrant. Patients classified as G2 and G3 accounted for 22.69% and 36.64%, respectively, while only 2.47% were G1. The molecular subtype HR+/HER2- was present in 44.2% of patients, with dual-positive and dual-negative subtypes accounting for 24.02% and 15.11%, respectively. The distribution of T stages from T1 to T4 was 11.21%, 34.57%, 18.21%, and 36.02%, and for N stages from N0 to N3 was 19.27%, 53.35%, 12.21%, and 15.18%, respectively. Bone metastases or lung metastases were present in 59.54% and 32.74% of patients, respectively, with only 8.95% having brain metastases. Regarding treatment, 23.92% of patients underwent primary tumor surgery, 27.40% received radiotherapy, and the majority (75.77%) received chemotherapy. Additionally, 124 BCLM patients from the JCH and CHSU cohorts were included as the external group. In contrast to the SEER cohort, all patients in the external cohorts were Asian, and only 8.06% underwent surgery.

We analyzed the correlation between variables and generated a heatmap, which indicated no multicollinearity among the variables ( Supplementary Figure S1 ). For clarity, the Cox regression results presented in Table 2 were performed on the overall training cohort to illustrate variable associations, while per-fold selections were used internally during cross-validation. Univariate Cox regression analysis revealed that age, race, marital status, median household income, histological type, tumor location, subtype, T stage, N stage, surgery, chemotherapy, bone metastases, brain metastases, and lung metastases significantly affected OS. Meanwhile, significantly affecting BCSS were age, race, marital status, median household income, histological type, subtype, tumor location, T stage, N stage, surgery, chemotherapy, bone metastases, brain metastases, and lung metastases.

Multivariate Cox regression analysis was employed to further explore independent risk factors affecting prognosis. GVIF suggested that there was no multicollinearity between the variables in the regression model ( Supplementary Table S2 ). Older age, other pathological type and tumor location, HR-/HER2-, T4 stage, and the presence of other distant metastases were associated with poorer OS. Socially, Black and divorced patients had worse OS, whereas patients with household income exceeding $70,000 exhibited better OS. Similarly, patients who underwent surgery and radiotherapy demonstrated improved OS. Moreover, factors indicating worse BCSS included older age, being black, being divorced, other pathological types, HR-/HER2-, T4 stage, and the presence of other distant metastases. In contrast, household income over $70,000, HR+/HER2+, N1 stage, and having undergone surgery and chemotherapy were associated with better BCSS.

We incorporated significant features into subsequent ML model construction to predict the 6-month, 1-year, 3-year, and 5-year OS and BCSS for BCLM patients. The performance of five ML models in the training and internal test groups was depicted in Figure 2 . The RF models demonstrated excellent performance in predicting 6-month OS (Training: AUC=0.850; Internal Test: AUC=0.765), 1-year OS (Training: AUC=0.840; Internal Test: AUC=0.774), 3-year OS (Training: AUC=0.871; Internal Test: AUC=0.774), and 5-year OS (Training: AUC=0.889; Internal Test: AUC=0.786). Additionally, the RF models also performed well in predicting BCSS at 6 months (Training: AUC=0.849; Internal Test: AUC=0.787), 1 year (Training: AUC=0.864; Internal Test: AUC=0.765), 3 years (Training: AUC=0.867; Internal Test: AUC=0.775), and 5 years (Training: AUC=0.862; Internal Test: AUC=0.774). Other ML models such as LR (6-month OS: AUC=0.779; 1-year OS: AUC = 0.748; 3-year OS: AUC = 0.735; 5-year OS: AUC = 0.768; 6-month BCSS: AUC=0.774; 1-year BCSS: AUC = 0.758; 3-year BCSS: AUC = 0.731; 5-year BCSS: AUC = 0.766), XGBoost (6-month OS: AUC=0.797; 1-year OS: AUC = 0.822; 3-year OS: AUC = 0.836; 5-year OS: AUC = 0.864; 6-month BCSS: AUC=0.786; 1-year BCSS: AUC = 0.798; 3-year BCSS: AUC = 0.822; 5-year BCSS: AUC = 0.855), DT (6-month OS: AUC=0.796; 1-year OS: AUC = 0.779; 3-year OS: AUC = 0.787; 5-year OS: AUC = 0.810; 6-month BCSS: AUC=0.793; 1-year BCSS: AUC = 0.777; 3-year BCSS: AUC = 0.826; 5-year BCSS: AUC = 0.857), and GBDT (6-month OS: AUC=0.800; 1-year OS: AUC = 0.813; 3-year OS: AUC = 0.761; 5-year OS: AUC = 0.865; 6-month BCSS: AUC=0.823; 1-year BCSS: AUC = 0.758; 3-year BCSS: AUC = 0.784; 5-year BCSS: AUC = 0.810) showed slightly inferior performance compared to RF in the training group. Moreover, the RF models’ accuracy (6-month OS: 0.776 and 0.745; 1-year OS: 0.761 and 0.714; 3-year OS: 0.764 and 0.711; 5-year OS: 0.807 and 0.732; 6-month BCSS: 0.793 and 0.759; 1-year BCSS: 0.782 and 0.708; 3-year BCSS: 0.784 and 0.714; 5-year BCSS: 0.759 and 0.786) and F1-scores (6-month OS: 0.615 and 0.498; 1-year OS: 0.686 and 0.650; 3-year OS: 0.828 and 0.768; 5-year OS: 0.863 and 0.823; 6-month BCSS: 0.600 and 0.526; 1-year BCSS: 0.703 and 605; 3-year BCSS:0.832 and 0.736; 5-year BCSS: 0.841 and 0.865) in both the training and internal testing groups were satisfactory ( Table 3 ).

DCA was used to evaluate the clinical utility of the models. The results indicated that all five models provided good net benefits in predicting survival rates at different time points, with the RF models showing a slightly higher net benefit in the training group ( Figure 3 ). Based on 1,000 bootstrap resamples, the RF models also exhibited consistent net benefits across a range of threshold probabilities ( Supplementary Table S3 ). Calibration curves further demonstrated the good agreement between the predicted survival outcomes by the RF model and actual outcomes ( Figure 4 ). The confusion matrix illustrated the performance of the RF classifier in the training and internal test groups ( Supplementary Figure S2 ). Therefore, the RF models were identified as the optimal models for predicting prognosis in BCLM patients.

To further validate the robustness and applicability of the RF models, we included a total of 124 BCLM patients from the JCH and CHSU cohorts. The RF models demonstrated excellent performance in the external test cohort, as evidenced by their high discriminatory ability ( Figures 5A–H ). Specifically, the AUC values were as follows: 6-month OS, 0.779 (95% CI: 0.688–0.871); 1-year OS, 0.803 (95% CI:0.720–0.887); 3-year OS, 0.812 (95% CI: 0.714–0.911); 5-year OS, 0.792 (95% CI: 0.536–0.995); 6-month BCSS, 0.802 (95% CI: 0.714–0.890); 1-year BCSS, 0.801 (95% CI: 0.713–0.889); 3-year BCSS, 0.784 (95% CI: 0.672–0.897); and 5-year BCSS, 0.815 (95% CI: 0.637–0.993). Moreover, calibration curves further indicated good agreement between the predicted and actual outcomes ( Figures 5I–P ). DCA revealed that the RF models also exhibited good net benefit in the external cohort ( Figure 6 ). Therefore, the RF models were considered the optimal tools for predicting the prognosis of BCLM patients.

Since 38.20% of patients had unknown grade, we excluded these cases to validate the robustness of the RF models. The RF models were retrained with the same hyperparameters as the original model ( Supplementary Table S4 ) and evaluated on the identical internal and external test groups. The retrained RF models achieved AUCs ranging from 0.858 to 0.899 in the training group, 0.741 to 0.831 in the internal test group, and 0.772 to 0.848 in the external test group ( Supplementary Figure S3 ). Given that performance remained comparable, these findings suggested that including patients with unknown grade did not materially compromise the robustness of our RF models.

The SHAP analysis provided an explanation of how the RF models predicted survival outcomes and calculated the importance of features. Figures 7A–H showed the SHAP values for each feature at different levels. As the feature values increased, the redder the color became, and vice versa, the bluer the color was. Additionally, we ranked the features of the models ( Figures 7I–P ). The higher feature ranking indicated that the feature was more important, which meant that the feature contributed more to the model. Overall, chemotherapy, age, subtype, and surgery were the most valued by the RF models for predicting 6-month, 1-year, 3-year, and 5-year OS and BCSS. Additionally, two representative samples were provided for each model to elucidate the composition of the model outputs. Each feature’s weight was depicted in blue or red, indicative of its contribution toward a positive outcome. The value of each feature denoted its proportionate influence on the final score. Supplementary Figures S4A–H presented the detailed scores of surviving patients predicted to be alive, whereas Supplementary Figures S4I–P displayed the scores of dying patients predicted to be dead.

Baseline characteristics of the PTS and non-PTS groups were shown in Supplementary Table S5 . The adjusted time-dependent Cox model demonstrated a reduced risk of OS events in the PTS group (HR=0.80, 95% CI: 0.72–0.88; P< 0.001; Table 4 ). Similarly, BCSS risk was also lower (HR=0.77, 95% CI: 0.69–0.86; P< 0.001). The Simon–Makuch curves illustrated survival differences over time between the PTS and non-PTS groups ( Figures 8A, B ). As shown in Figures 8C, D , even after accounting for time-varying effects, the instantaneous HRs for both OS and BCSS remained consistently< 1 across the majority of the follow-up period.

The piecewise Cox model further supported these findings: for OS, the HR was 0.83 (95% CI: 0.75–0.92) within the first 24 months and 0.73 (95% CI: 0.61–0.88) thereafter; for BCSS, the HR was 0.80 (95% CI: 0.72–0.89) within 0–24 months and 0.68 (95% CI: 0.56–0.83) beyond 24 months ( Supplementary Table S6 ). The 2-month landmark analysis similarly showed that patients receiving PTS after systemic therapy had improved OS (HR=0.60, 95% CI: 0.37–0.97; p = 0.038) and BCSS (HR=0.60, 95% CI: 0.36–0.98; p = 0.043), consistent with the primary model ( Supplementary Table S7 ). We further quantified the robustness of our findings using the E-value. For OS and BCSS, the E-values were 1.81 and 1.92, respectively (lower bound E-values: OS=1.53; BCSS=1.60). This suggested that to nullify the observed survival benefits of PTS, an unmeasured confounder would need to be associated with both PTS treatment and the outcomes with an HR ≥ 1.81 or 1.92. Similarly, in the 2-month landmark analysis, the E-value for OS (HR=0.60) was 2.72 (lower bound: 1.21), and for BCSS (HR=0.60) was 2.72 (lower bound: 1.16). These results indicated that only an unmeasured confounder with at least moderate strength could substantially alter our conclusions. In summary, our findings confirm the survival benefit of PTS in improving both OS and BCSS.