Inflammation-lipid metabolism biomarker-based machine learning models for predicting postoperative pulmonary metastasis in colorectal cancer

Abstract

BACKGROUND

Colorectal cancer (CRC) ranks third in global incidence and second in cancer-related mortality. The lung is the second most common distant metastatic site, and pulmonary metastasis substantially worsens long-term survival. Although tumor-associated inflammation and lipid metabolic reprogramming synergistically drive the metastatic cascade, systematic integration of these biomarkers into predictive models remains limited. Whether preoperative routine hematological indices alone can achieve preliminary risk stratification - particularly before pathological results become available or in resource-limited primary care settings - has not been systematically investigated.

AIM

To develop and validate an interpretable machine learning model integrating inflammation-lipid metabolism biomarkers for CRC pulmonary metastasis prediction.

METHODS

A total of 1246 CRC patients from two centers were enrolled. The development cohort was split 7:3 into training (n = 621) and internal validation (n = 267) sets; an external cohort (n = 358) independently assessed generalizability. Sixteen candidate predictors spanning demographics (5), systemic inflammatory indices (8), and lipid metabolism indices (3) were incorporated via a four-layer progressive feature combination strategy. Fourteen machine learning algorithms and three ensemble strategies (Soft Voting, Weighted Voting, Stacking) were systematically evaluated, with SHapley Additive exPlanations (SHAP) analysis applied for multilevel interpretability.

RESULTS

Progressive modeling demonstrated incrementally improved performance with sequential introduction of inflammatory and lipid metabolism features. LightGBM achieved the best individual performance [area under the curve (AUC) = 0.7507, 95% confidence interval: 0.6872-0.8157]. Ensemble methods further enhanced performance; Weighted Voting attained the highest internal validation AUC (0.7513) and maintained robust external generalization (AUC = 0.7112, Brier = 0.1850). SHAP global analysis identified lymphocyte-to-C-reactive protein ratio (LCR), body mass index (BMI), and apolipoprotein B-to-apolipoprotein A1 ratio as the top three predictors, indicating that inflammatory-immune balance, metabolic status, and lipoprotein dysregulation constitute core dimensions driving risk prediction. Dependence plots revealed key thresholds (LCR < 0.5, BMI > 25) informing potential clinical risk stratification.

CONCLUSION

This multimodal inflammation-lipid metabolism model effectively predicts CRC pulmonary metastasis risk, with cross-center transportability confirmed through external validation. SHAP analysis supports the model's clinical credibility and biological plausibility. All predictors are obtainable through routine clinical testing without specialized equipment, enabling deployment in primary care settings. The model has been implemented as an interactive web application for individualized risk assessment and visual decision support.

Key Words: Colorectal cancer; Pulmonary metastasis; Machine learning; Inflammatory biomarkers; Lipid metabolism; Prediction model

Core Tip: No prior model predicts colorectal cancer pulmonary metastasis using only preoperative routine blood-derived inflammation and lipid metabolism biomarkers. In this 1246-patient, dual-center study, we systematically compared 14 machine learning algorithms and three ensemble strategies - the most extensive benchmark reported for this task. The optimal Weighted Voting ensemble achieved validated discrimination [internal area under the curve (AUC) = 0.7513; external AUC = 0.7112] using 16 readily obtainable predictors. SHapley Additive exPlanations interpretability analysis revealed the lymphocyte-to-C-reactive protein ratio and body mass index as dominant predictors, offering a practical risk-stratification tool before pathological staging becomes available.

INTRODUCTION

Colorectal cancer (CRC) is the third most common malignancy worldwide and ranks second in cancer-related mortality, representing a major public health burden[1]. In countries and regions undergoing socioeconomic transition, the incidence of CRC continues to rise, driven by shifts in dietary patterns and the increasing prevalence of obesity[2]. Distant metastasis is the leading cause of death in CRC patients, with the liver and lung being the most common metastatic sites[3]. Epidemiological studies have shown that pulmonary metastasis occurs in approximately 10%-15% of CRC patients, manifesting either synchronously at initial diagnosis or metachronously following curative resection[4,5]. A subset of patients with resectable pulmonary metastases may achieve meaningful long-term survival following local treatment, although high-quality evidence remains limited and the criteria for identifying optimal candidates remain controversial[4]. Therefore, early identification and accurate risk stratification of pulmonary metastasis in CRC patients are indispensable for the development of individualized surveillance strategies and optimization of therapeutic decision-making[6].

Tumor-associated inflammation is recognized as one of the hallmarks of cancer[7], and the systemic inflammatory response not only reflects the interplay between the tumor and host immune system but also facilitates the colonization of circulating tumor cells in distant organs through remodeling of the pre-metastatic niche[8,9]. Peripheral blood-based composite inflammatory indices, including the lymphocyte-to-C-reactive protein ratio (LCR), systemic inflammation response index (SIRI), platelet-to-lymphocyte ratio (PLR), and modified Glasgow Prognostic Score (mGPS), have been demonstrated to be closely associated with prognosis and metastatic risk across multiple solid tumor types[10,11]. Concurrently, metabolic reprogramming is a key mechanism by which tumor cells adapt to the metastatic microenvironment, with aberrant lipid metabolism being particularly prominent[12]; lipids not only provide tumor cells with energy and membrane structural components but also participate in regulating immune evasion and angiogenesis[13]. Metabolic indices reflecting insulin resistance, such as the triglyceride-glucose-body mass index (TyG-BMI)[14] and the apolipoprotein B-to-apolipoprotein A1 ratio (ApoB/ApoA1)[15], have also been associated with malignant tumor progression. Extensive biological crosstalk exists between inflammation and lipid metabolism: Pro-inflammatory cytokines secreted by adipose tissue under conditions of obesity can exacerbate systemic inflammation, while chronic inflammation can in turn disrupt lipid metabolic homeostasis[16,17]. However, studies systematically integrating both dimensions for CRC pulmonary metastasis risk prediction remain limited. Given that inflammation and lipid metabolism constitute a bidirectional positive feedback loop - wherein pro-inflammatory adipokines exacerbate systemic inflammation while chronic inflammation disrupts lipid homeostasis - their joint integration is biologically imperative for comprehensive risk characterization, as capturing only a single dimension would omit the synergistic contribution of this crosstalk to the metastatic cascade.

Despite the demonstrated potential predictive value of these biomarkers, methodological limitations persist in prior research. Previous CRC pulmonary metastasis risk prediction studies have predominantly relied on traditional statistical methods, presenting results in nomogram format[18], with inherent limitations in capturing nonlinear effects and higher-order interactions[19]. Machine learning approaches can effectively capture complex nonlinear feature relationships[20], yet existing models have largely focused on comparing a small number of classical algorithms and leave room for improvement in the systematic capture of feature interactions. Notably, oncological variables such as tumor-node-metastasis (TNM) staging, pathological grading, and lymphovascular/perineural invasion are recognized prognostic factors for CRC metastasis and have been incorporated into comprehensive prediction models in multiple studies[18-21]. However, the acquisition of these variables depends on surgical pathology or biopsy results. In current clinical practice, postoperative surveillance intensity and imaging follow-up frequency are primarily determined by pathological staging, meaning that risk stratification is deferred until pathological results become available. Two practically important clinical scenarios remain underserved by pathology-dependent models: (1) The preoperative or initial-diagnosis window - typically spanning days to weeks - during which surgical pathology results are pending and no quantitative risk assessment tool is available to guide early clinical decision-making; and (2) Resource-limited primary care or community hospital settings where advanced pathological evaluation may be delayed, incomplete, or inaccessible due to infrastructure constraints. Therefore, exploring a set of predictor variables obtainable solely through preoperative routine physical examination and standard laboratory testing, and evaluating their risk stratification capacity in the absence of pathological information, carries practical significance for expanding the clinical coverage and temporal reach of CRC pulmonary metastasis risk assessment. This study does not aim to replace comprehensive prediction models incorporating full pathological information but rather to develop a complementary tool capable of providing preliminary risk signals prior to or in the absence of pathological data - thereby bridging the temporal and infrastructural gaps that limit the applicability of pathology-dependent models.

Against this background, the present study aimed to develop and validate a machine learning model based on multimodal inflammation-lipid metabolism biomarkers using multicenter real-world retrospective cohort data, for predicting the risk of pulmonary metastasis in CRC patients. Specifically, we designed a four-layer progressive feature combination strategy to quantify the incremental predictive value of each biomarker category based on 16 prespecified candidate variables spanning three domains (demographic characteristics, systemic inflammatory indices, and lipid metabolism indices); employed 14 machine learning algorithms spanning diverse learning paradigms for model development and validation, and further explored multiple ensemble strategies to optimize predictive performance; assessed cross-center generalizability through an independent external validation cohort; conducted multilevel model interpretability analysis using the SHapley Additive exPlanations (SHAP) framework; and deployed the final model as an interactive web application to facilitate clinical translation (Figure 1).

Figure 1 Flow chart of the study design. CRC: Colorectal cancer; KNN: K-Nearest Neighbors; SVM: Support Vector Machine; MLP: Multi-Layer Perceptron; GBDT: Gradient Boosting Decision Tree; AUC: Area under the curve; AUROC: Area under the receiver operating characteristic; PPV: Positive predictive value; NPV: Negative predictive value; SHAP: SHapley Additive exPlanations; LCR: Lymphocyte-to-C-reactive protein ratio; BMI: Body mass index; ApoB/ApoA1: Apolipoprotein B-to-apolipoprotein A1 ratio.

MATERIALS AND METHODS

Study design and patient population

This was a multicenter retrospective cohort study designed to develop and validate a machine learning model based on multimodal inflammation-lipid metabolism biomarkers for predicting the risk of pulmonary metastasis in patients with CRC. The development cohort (Guang’anmen Hospital, China Academy of Chinese Medical Sciences; Center A) enrolled 888 CRC patients treated between April 2018 and February 2023; the external validation cohort (Gansu Provincial Cancer Hospital; Center B) enrolled 358 CRC patients treated during the same period.

Inclusion criteria: (1) Histopathologically confirmed colorectal adenocarcinoma; (2) CRC as the sole diagnosed primary malignancy; (3) Patients underwent curative resection of the primary tumor (R0 resection); (4) Staged as M0 (no distant metastasis) at initial diagnosis, with no documented evidence of recurrence or metastasis prior to diagnosis; (5) Age ≥ 18 years; and (6) Availability of at least 3 years of standardized follow-up records, or occurrence of a confirmed pulmonary metastasis event during the follow-up period.

Exclusion criteria: (1) Development of pulmonary metastasis within 6 months after curative surgical resection; (2) Missing proportion of baseline clinical data or laboratory parameters exceeding 20%; follow-up duration less than 3 years without an observed pulmonary metastasis event; or death from non-pulmonary metastasis causes during the follow-up period; (3) Receipt of neoadjuvant chemotherapy or radiotherapy prior to surgery; (4) Concurrent diagnosis of another primary malignancy; and (5) Presence of severe underlying pulmonary disease (e.g., diffuse pulmonary interstitial fibrosis, active pulmonary tuberculosis), precluding reliable differentiation between metastatic lesions and pre-existing pulmonary abnormalities on imaging.

This study was conducted in accordance with the principles of the Declaration of Helsinki (2024 revision)[22] and was approved by the Ethics Committee of Guang’anmen Hospital, China Academy of Chinese Medical Sciences (approval number: 2025-238-KY), which included approval for the sub-center at Gansu Provincial Cancer Hospital. Study design, analysis, and reporting adhered to the recommendations of the Prediction model Risk of Bias Assessment Tool for Artificial Intelligence (PROBAST+AI)[23] and the Transparent Reporting of a Multivariable Prediction Model for Individual Prognosis or Diagnosis for Artificial Intelligence (TRIPOD+AI)[24] to ensure methodological rigor, bias control, and reporting transparency.

Outcome definition

The primary outcome was the occurrence of metachronous pulmonary metastasis in CRC patients, defined as a pulmonary metastatic lesion first identified ≥ 6 months after curative treatment (typically R0 resection) of the primary CRC, confirmed by contrast-enhanced chest CT or PET-CT[25]. According to the Expert Consensus on Multidisciplinary Treatment of Pulmonary Metastasis from Colorectal Cancer (2018 edition)[26], radiological diagnostic criteria for pulmonary metastasis included: Solid or ground-glass nodules in the peripheral zones and lower lobes of bilateral lungs, > 5 mm in diameter, with well-defined margins, lobulation, or short spiculation. For indeterminate pulmonary nodules, independent review was performed by ≥ 2 radiologists, integrating clinical risk factors, dynamic follow-up changes, and pathological examination when necessary. The inclusion criteria ensured that all patients had a definitive outcome status (pulmonary metastasis within 3 years or completion of 3-year follow-up without an outcome event); therefore, no censored data existed and survival analysis was not required.

Candidate predictor variables and selection rationale

A total of 16 candidate predictor variables were included in this study, covering three domains: Demographic characteristics, systemic inflammatory indices, and lipid metabolism parameters (Supplementary Table 1). Variable selection was guided by a comprehensive consideration of clinical evidence and biological plausibility (all candidate variables had been reported in prior literature as being associated with CRC prognosis or metastatic risk), clinical accessibility (all obtainable through routine physical examination and standard laboratory testing without requiring specialized molecular assay equipment, suitable for deployment in primary care settings), and informational complementarity (the three categories of indices characterized the tumor-host interaction from three distinct dimensions - host baseline status, immune-inflammatory microenvironment, and metabolic reprogramming - thereby minimizing information redundancy). All predictor variables were collected within a 14-day window prior to curative surgery.

Both centers’ clinical laboratories had passed the External Quality Assessment administered by the National Center for Clinical Laboratories of the National Health Commission, with satisfactory performance across assessed analytes during the evaluation cycles. As the two centers utilized different analytical platforms and reagent systems, systematic bias between raw assay values may exist. This study did not perform cross-center standardization or calibration of raw values; instead, external validation was employed to evaluate the model’s transportability across different analytical platforms.

Demographic and basic clinical characteristics (5 variables): This study included five demographic variables: Age, sex, body mass index (BMI), smoking history, and alcohol consumption history. Advanced age is associated with immunosenescence and reduced treatment tolerance and is an established risk factor for adverse CRC outcomes[27]. Epidemiological studies have shown that male CRC patients carry a higher risk of distant metastasis than female patients, potentially attributable to the anti-inflammatory protective effects of estrogen[28]. BMI reflects not only nutritional status but also the capacity of obesity to promote tumor metastasis through mechanisms including chronic low-grade inflammation, insulin resistance, and adipokine dysregulation[29]. Tobacco exposure can induce DNA damage and promote epithelial-mesenchymal transition and angiogenesis, thereby accelerating tumor progression[30]. Chronic alcohol consumption has also been associated with increased CRC risk and elevated systemic inflammation[31].

Systemic inflammatory indices (8 variables): Tumor-associated inflammation is widely recognized as one of the hallmarks of cancer. Systemic inflammatory responses not only reflect the interaction between the tumor and host immune system but also promote metastatic colonization of distant organs by remodeling the pre-metastatic niche[9]. Accordingly, this study incorporated eight composite inflammatory indices.

Among these, the SIRI (SIRI = neutrophils × monocytes/Lymphocytes) comprehensively reflects the balance between innate and adaptive immunity[32]; the PLR captures platelet-mediated shielding of circulating tumor cells[33]; and the lymphocyte-to-monocyte ratio (LMR) and prognostic nutritional index (PNI) (PNI = albumin + 5 × lymphocyte count) reflect host anti-tumor capacity from the perspectives of immune-inflammatory balance and immune-nutritional status, respectively[34]. The C-reactive protein-to-albumin ratio (CAR) and LCR both integrate information from CRP, an acute-phase protein: The former reflects imbalance along the inflammation-nutrition axis[35], while the latter comprehensively reflects the dynamic relationship between immune function and inflammatory burden. Low LCR indicates augmented inflammation coupled with immunosuppression and is closely associated with poor CRC prognosis[10]. The mGPS is a well-established inflammation-nutrition scoring system based on CRP and albumin levels, validated as an independent prognostic factor for solid tumors in multiple studies[11]. The fibrinogen-to-albumin ratio (FAR) bridges the coagulation and inflammation dimensions, as fibrinogen promotes angiogenesis and metastasis within the tumor microenvironment[36].

Lipid metabolism-related indices (3 variables): Lipid metabolic reprogramming is a key biological process by which tumor cells adapt to the metastatic microenvironment[12]. Lipids not only provide rapidly proliferating tumor cells with energy and membrane structural components but also function as signaling molecules involved in regulating inflammation, immune evasion, and angiogenesis[13]. Accordingly, this study incorporated three composite lipid metabolism indices.

The TyG-BMI [TyG-BMI = ln(triglycerides × fasting blood glucose/2) × BMI] serves as a surrogate marker for insulin resistance, and obesity combined with insulin resistance can promote CRC proliferation and metastasis through activation of the PI3K/Akt/mTOR pathway[14]. The ApoB/ApoA1 reflects the balance between pro-atherogenic lipoproteins [low-density lipoprotein (LDL), very low-density lipoprotein (VLDL)] and anti-atherogenic lipoproteins [high-density lipoprotein (HDL)]; an elevated ApoB/ApoA1 ratio is independently associated with increased CRC risk[15] and may indirectly promote tumor metastasis by fostering chronic low-grade inflammation and oxidative stress. The fibrinogen-to-high-density lipoprotein cholesterol ratio (FHR) is a composite biomarker integrating fibrinogen and HDL-C, reflecting the combined burden of coagulation/inflammation and lipid dysregulation[37].

Regarding collinearity, BMI appears simultaneously as a demographic characteristic and as a component variable of TyG-BMI, resulting in high correlation between the two. The decision to retain both was based on the following rationale: (1) BMI reflects overall adiposity as a systemic measure, whereas TyG-BMI integrates information from three components - blood glucose, triglycerides, and BMI, characterizing metabolic status from distinct biological perspectives; (2) The tree-based algorithms employed in this study possess inherent robustness to collinearity, automatically mitigating its effects through the stochasticity of feature splitting; and (3) The SHAP method can reasonably allocate contribution shares among correlated features when computing individual feature contributions.

Hierarchical feature modeling strategy

To systematically evaluate the incremental predictive value of different categories of biomarkers for CRC pulmonary metastasis, a four-layer progressive feature combination strategy was designed:

Layer 1 (baseline demographic layer): Included only 5 demographic features (age, sex, BMI, smoking, alcohol consumption) as a baseline reference model.

Layer 2 (demographic + inflammatory layer): Added 8 systemic inflammatory indices to layer 1, totaling 13 features, to evaluate the incremental predictive contribution of inflammatory biomarkers.

Layer 3 (demographic + lipid metabolism layer): Added 3 lipid metabolism indices to layer 1, totaling 8 features, to evaluate the incremental predictive contribution of lipid metabolism biomarkers.

Layer 4 (full-feature layer): Integrated all 16 features to evaluate the combined predictive efficacy of inflammatory and lipid metabolism biomarkers.

By comparing performance differences among the best models at each layer and DeLong test results, the incremental contribution of inflammatory and lipid metabolism biomarkers to the baseline model could be quantitatively assessed, along with whether synergistic effects emerged from their combination.

Sample size adequacy assessment

Sample size adequacy was evaluated using the minimum sample size criteria for multivariable prediction models proposed by Riley et al[38]. For the full-feature model (16 variables), the development cohort comprised 888 patients with 285 events, exceeding the minimum required thresholds of 567 patients and 177 events, thereby confirming adequate sample size.

Data partitioning

Patients from Guang’anmen Hospital were randomly split in a 7:3 ratio (StratifiedShuffleSplit, random seed = 42) into a training set (n = 621; 199 events; event rate 32.0%) and an internal validation set (n = 267; 86 events; event rate 32.2%), ensuring balanced event rates between the two subsets. All 358 patients from Gansu Provincial Cancer Hospital (103 events; event rate 28.8%) served as the independent external validation set. For data with less than 20% missingness, multiple imputation was performed at the raw variable level using Multiple Imputation by Chained Equations, employing predictive mean matching for continuous variables and logistic regression (LR) for binary variables.

Feature collinearity diagnostics

Prior to modeling, multicollinearity diagnostics were conducted on all 16 candidate variables. A Pearson correlation coefficient matrix was used to assess pairwise linear associations among variables, with |r| > 0.9 set as the high collinearity warning threshold. Variance inflation factors (VIF) were also calculated, with VIF > 10 designated as the criterion for severe multicollinearity[39].

Machine learning algorithms

This study systematically compared 14 machine learning algorithms spanning multiple learning paradigms, including LR, Generalized Additive Model (GAM), K-Nearest Neighbors (KNN), Support Vector Machine (SVM), Gaussian Naive Bayes (NB), Decision Tree (DT), Random Forest (RF), Extra Trees, Gradient Boosting DT (GBDT), XGBoost, LightGBM, AdaBoost, Rotation Forest (RotF), and Multi-Layer Perceptron (MLP). For algorithms requiring feature scale normalization (LR, GAM, KNN, SVM, MLP), StandardScaler standardization was applied prior to training. Standardization parameters were fitted exclusively on the training set and subsequently applied to the validation and external validation sets, strictly preventing data leakage[40]. Systematic comparison of multiple algorithms is a recommended practice in contemporary prediction modeling research[41], facilitating the identification of the learning strategy best suited to a given data structure.

Systematic performance evaluation and comparison of all candidate models were conducted on the internal validation set. As multiple model evaluations on the same validation set may introduce optimism bias due to multiple comparisons, 5 repetitions of 5-fold stratified cross-validation were performed on the training set, and the resulting algorithm rankings were compared for consistency with those obtained on the internal validation set.

Model evaluation metrics

Each model underwent multidimensional performance evaluation on the internal validation set, encompassing discrimination [area under the curve (AUC) with Bootstrap 95% confidence interval (95%CI), area under the precision-recall curve], calibration (Brier score, calibration curve), clinical utility (decision curve analysis), and classification performance (accuracy, F1, sensitivity, specificity, positive predictive value (PPV), and negative predictive value (NPV) based on the Youden index-derived threshold). Additionally, a composite score, ABscore = 0.5 × AUROC_norm + 0.5 × (1 - Brier_norm), was introduced as an integrated metric[42].

Cross-layer comparison and deLong test

The DeLong test[43] was employed to quantitatively assess the statistical significance of AUC differences across feature layers. For each algorithm, AUC differences were compared between layer 1 and layer 2, layer 1 and layer 3, and layer 1 and layer 4, with P < 0.05 considered statistically significant.

Ensemble model strategies

Ensemble models were constructed from the top 4 individual models ranked by AUC in the best-performing layer[44,45]: (1) Soft Voting: Unweighted arithmetic mean of the predicted probabilities from the 4 base models; (2) Weighted Voting: Weighted average of predicted probabilities using each base model’s internal validation AUC as the weight; and (3) Stacking ensemble: A two-layer architecture in which the first layer generated out-of-fold predictions via 5-fold cross-validation, and the second layer employed LR as the meta-learner. This two-stage approach - first systematically comparing individual algorithms to identify the best-performing base learners, then combining them through ensemble strategies - follows the principle of bias-variance decomposition: Whereas individual models may overfit to idiosyncratic patterns in the training data, ensembles of diverse models can attenuate such variance, thereby improving prediction stability and cross-center generalizability.

External validation

The final models were retrained on the complete dataset from Guang’anmen Hospital (Center A), and external generalization performance was evaluated on the Gansu Provincial Cancer Hospital (Center B) dataset. In external validation, all classification performance metrics (accuracy, sensitivity, specificity, F1, PPV, NPV) were calculated using the optimal classification threshold determined by the Youden index on the internal validation set, without re-optimization on the external validation set, thereby simulating the prospective application scenario of real-world clinical deployment.

SHAP interpretability analysis

SHAP is a unified model interpretation framework based on Shapley values from cooperative game theory[46]. In this study, KernelExplainer-based SHAP analysis was applied to the best-performing ensemble model in the optimal layer, including: Global feature importance ranking, beeswarm plots, dependence plots, local explanations (waterfall plots and force plots), SHAP interaction value heatmaps, partial dependence plots, and feature correlation network plots.

Clinical application tool

To facilitate clinical translation and practical utility, the final prediction model was deployed as an interactive web application built on the Streamlit framework. The application utilized the top 10 features ranked by SHAP importance as input variables from the optimal model, balancing predictive accuracy with clinical operational simplicity.

Statistical analysis

Continuous variables were expressed as M (Q1, Q3), with between-group comparisons performed using the Mann-Whitney U test or Kruskal-Wallis test. Categorical variables were expressed as n (%), with between-group comparisons performed using the χ² test. All machine learning modeling and analyses were completed in a Python 3.10 environment. P < 0.05 was considered statistically significant.

RESULTS

Baseline characteristics

The 888 CRC patients from Guang’anmen Hospital were split in a 7:3 ratio into a training set (n = 621) and an internal validation set (n = 267), while patients from Gansu Provincial Cancer Hospital (n = 358) served as the external validation set. Comparative results of baseline characteristics across the three groups are presented in Supplementary Table 2. The event rate was 32.0% (199/621) in the training set, 32.2% (86/267) in the internal validation set, and 28.8% (103/358) in the external validation set, indicating that stratified random splitting effectively ensured balanced event rates within the development cohort.

Feature collinearity assessment

Pearson correlation analysis revealed a very high positive correlation between BMI and TyG-BMI (r = 0.933) and a moderately high positive correlation between FAR and FHR (r = 0.842) (Figure 2A). VIF analysis showed that TyG-BMI (VIF = 9.39) and BMI (VIF = 9.06) approached but did not exceed the threshold of 10.0, while VIF values for all remaining variables were below 5.0 (Figure 2B).

Figure 2 Feature collinearity assessment. A: Pearson correlation coefficient heatmap; B: Variance inflation factor. BMI: Body mass index; SIRI: Systemic inflammation response index; LMR: Lymphocyte-to-monocyte ratio; PLR: Platelet-to-lymphocyte ratio; CAR: C-reactive protein-to-albumin ratio; LCR: Lymphocyte-to-C-reactive protein ratio; mGPS: Modified Glasgow Prognostic Score; PNI: Prognostic nutritional index; FAR: Fibrinogen-to-albumin ratio; TyG-BMI: Triglyceride-glucose-body mass index; FHR: Fibrinogen-to-high-density lipoprotein cholesterol ratio; ApoB/ApoA1: Apolipoprotein B-to-apolipoprotein A1 ratio.

Hierarchical machine learning modeling results

Summary of layer 1-3 results: Layer 1 (demographic features only): The best-performing model was KNN (AUC = 0.7038, ABscore = 0.8361), followed by MLP (AUC = 0.6957, ABscore = 0.8754). Overall predictive performance in this layer was limited, with AUC values ranging from 0.5408 to 0.7038 across the 14 algorithms (Supplementary Table 3), indicating that demographic features alone were insufficient to adequately capture pulmonary metastasis risk.

Layer 2 (demographic features + inflammatory indices): The best-performing model was LightGBM (AUC = 0.7448, ABscore = 0.9953), representing an AUC improvement of 0.0410 over the best layer 1 model. The introduction of inflammatory indices was associated with improved predictive performance for most algorithms (Supplementary Table 4), suggesting potential incremental value of systemic inflammatory biomarkers in CRC pulmonary metastasis prediction.

Layer 3 (demographic features + lipid metabolism indices): The best-performing model was MLP (AUC = 0.7459, ABscore = 0.8301); with only 3 lipid metabolism indices, this layer approached the performance of the 13-feature model in layer 2, suggesting that lipid metabolism biomarkers possess high predictive information density (Supplementary Table 5).

Layer 4: Full-feature model: After integrating all 16 features, LightGBM ranked first with AUC = 0.7507, followed closely by RotF (AUC = 0.7454) and GBDT (AUC = 0.7367). LightGBM simultaneously achieved the lowest Brier score (0.1856) and the highest ABscore (1.0000). It should be noted that ABscore = 1.0000 indicates that this model achieved the best normalized AUC and Brier values among all 14 candidate models within layer 4, i.e., LightGBM ranked first in both discrimination and calibration dimensions, which is statistically plausible (Table 1).

Table 1 Layer 4 (full-feature model): Internal validation performance comparison across 14 algorithms.

Algorithm	AUC	AUC_95CI	AUPRC	Brier	ABscore	Accuracy	Sensitivity	Specificity	F1	PPV	NPV	Opt_Threshold
LightGBM	0.7507	0.6872-0.8157	0.5592	0.1856	1.0	0.7191	0.7907	0.6851	0.6445	0.544	0.8732	0.3144
RotF	0.7454	0.6783-0.8091	0.5333	0.1875	0.9583	0.7266	0.7326	0.7238	0.6332	0.5575	0.8506	0.3459
GBDT	0.7367	0.6723-0.7967	0.5439	0.1876	0.9113	0.6704	0.7791	0.6188	0.6036	0.4926	0.855	0.2922
RF	0.7252	0.6620-0.7862	0.5171	0.1935	0.8075	0.6779	0.8023	0.6188	0.6161	0.5	0.8682	0.3333
SVM	0.7222	0.6581-0.7879	0.5327	0.1938	0.7892	0.6779	0.7093	0.663	0.5865	0.5	0.8276	0.294
AdaBoost	0.7172	0.6472-0.7794	0.509	0.1976	0.7349	0.6479	0.8023	0.5746	0.5948	0.4726	0.8595	0.3209
ExtraTrees	0.7169	0.6519-0.7812	0.5125	0.1988	0.7248	0.6966	0.686	0.7017	0.593	0.5221	0.8247	0.3487
KNN	0.7162	0.6553-0.7778	0.4889	0.1977	0.7287	0.6517	0.7326	0.6133	0.5753	0.4737	0.8284	0.4286
MLP	0.7132	0.6436-0.7778	0.522	0.2539	0.3016	0.6667	0.7442	0.6298	0.5899	0.4885	0.8382	0.0523
DT	0.7096	0.6463-0.7688	0.4685	0.1951	0.7131	0.6667	0.7442	0.6298	0.5899	0.4885	0.8382	0.4217
XGBoost	0.6956	0.6276-0.7603	0.4896	0.1971	0.6243	0.633	0.7791	0.5635	0.5776	0.4589	0.843	0.2972
GAM	0.6848	0.6158-0.7573	0.4956	0.1997	0.5487	0.6554	0.6395	0.663	0.5446	0.4741	0.7947	0.3507
LR	0.6689	0.5988-0.7375	0.4832	0.2022	0.446	0.6667	0.5581	0.7182	0.5189	0.4848	0.7738	0.3677
NB	0.656	0.5881-0.7276	0.4897	0.2376	0.1198	0.6592	0.593	0.6906	0.5285	0.4766	0.7812	0.3976

AUC: Area under the curve; AUPRC: Area under the precision-recall curve; PPV: Positive predictive value; NPV: Negative predictive value; RotF: Rotation Forest; GBDT: Gradient Boosting Decision Tree; RF: Random Forest; SVM: Support Vector Machine; NB: Naive Bayes; KNN: K-Nearest Neighbors; MLP: Multi-Layer Perceptron; DT: Decision Tree; LR: Logistic regression; GAM: Generalized Additive Model.

The ROC curves, precision-recall curves (PRC), decision curve analysis (DCA), calibration curves, and Brier score bar charts for the four-layer models are summarized in Figures 3, 4, 5, 6, and 7. As shown in Figure 3, the ROC curves shifted progressively toward the upper-left corner with advancing feature layers, with AUC values increasing incrementally. The PRC curves in Figure 4 exhibited a similar trend, indicating that the introduction of multimodal features also improved the precision-recall tradeoff in the class-imbalanced setting. The DCA in Figure 5 demonstrated that the best models in layers 2-4 provided positive net benefit exceeding the “treat all” and “treat none” strategies across a clinically relevant threshold probability range of approximately 10%-40%. The calibration curves in Figure 6 indicated good calibration for the optimal models in layers 2-4, with predicted probabilities broadly concordant with observed frequencies. The Brier score bar chart in Figure 7 visually displayed the combined calibration-discrimination performance of each algorithm across layers, with LightGBM in layer 4 achieving the lowest Brier score (0.1856).

Figure 3 Receiver operating characteristic curves for the four-layer models. A: Layer 1 receiver operating characteristic (ROC) curves (demographics); B: Layer 2 ROC curves (demographics + inflammation); C: Layer 3 ROC curves (demographics + lipid); D: Layer 4 ROC curves (all features). DT: Decision Tree; RotF: Rotation Forest; GAM: Generalized Additive Model; RF: Random Forest; LR: Logistic regression; NB: Naive Bayes; KNN: K-Nearest Neighbors; SVM: Support Vector Machine; MLP: Multi-Layer Perceptron; GBDT: Gradient Boosting Decision Tree; TPR: True positive rate; FPR: False positive rate.

Figure 4 Precision-recall curves for the four-layer models. A: Layer 1 precision-recall curves (demographics); B: Layer 2 precision-recall curves (demographics + inflammation); C: Layer 3 precision-recall curves (demographics + lipid); D: Layer 4 precision-recall curves (all features). PPV: Positive predictive value; RotF: Rotation Forest; GBDT: Gradient Boosting Decision Tree; MLP: Multi-Layer Perceptron; RF: Random Forest; KNN: K-Nearest Neighbors; GAM: Generalized Additive Model; SVM: Support Vector Machine; DT: Decision Tree; LR: Logistic regression; NB: Naive Bayes.

Figure 5 Decision curve analysis for the four-layer models. A: Layer 1 decision curve analysis (demographics); B: Layer 2 decision curve analysis (demographics + inflammation); C: Layer 3 decision curve analysis (demographics + lipid); D: Layer 4 decision curve analysis (all features). KNN: K-Nearest Neighbors; MLP: Multi-Layer Perceptron; GBDT: Gradient Boosting Decision Tree; RotF: Rotation Forest; SVM: Support Vector Machine; RF: Random Forest.

Figure 6 Calibration curves for the four-layer models. A: Layer 1 calibration curves (demographics); B: Layer 2 calibration curves (demographics + inflammation); C: Layer 3 calibration curves (demographics + lipid); D: Layer 4 calibration curves (all features). KNN: K-Nearest Neighbors; MLP: Multi-Layer Perceptron; GBDT: Gradient Boosting Decision Tree; RotF: Rotation Forest; RF: Random Forest; DT: Decision Tree; SVM: Support Vector Machine.

Figure 7 Brier scores for the four-layer models. A: Layer 1 Brier scores (demographics); B: Layer 2 Brier scores (demographics + inflammation); C: Layer 3 Brier scores (demographics + lipid); D: Layer 4 Brier scores (all features). KNN: K-Nearest Neighbors; MLP: Multi-Layer Perceptron; GBDT: Gradient Boosting Decision Tree; RotF: Rotation Forest; RF: Random Forest; DT: Decision Tree; SVM: Support Vector Machine; LR: Logistic regression; NB: Naive Bayes; GAM: Generalized Additive Model.

Cross-layer incremental analysis: DeLong test results showed that the ΔAUC between layer 1 and layer 4 was positive for most algorithms but did not reach statistical significance (P > 0.05), which may be attributable to limited statistical power of the DeLong test’s asymptotic framework given the internal validation set sample size. Consequently, the incremental predictive value of inflammatory and lipid metabolism indices should be regarded as exploratory findings warranting confirmation in larger studies. Figure 8A presents the AUC comparison among the best models at each layer, and the DeLong test heatmap in Figure 8B provides a detailed P-value matrix for cross-layer comparisons across all algorithms.

Figure 8 Cross-layer incremental analysis. A: Area under the curve (AUC) comparison of best models across layers; B: P-value matrix of cross-layer AUC comparisons for each algorithm. AUC: Area under the curve; KNN: K-Nearest Neighbors; MLP: Multi-Layer Perceptron; LR: Logistic regression; GAM: Generalized Additive Model; GBDT: Gradient Boosting Decision Tree; RotF: Rotation Forest; RF: Random Forest; DT: Decision Tree; SVM: Support Vector Machine; NB: Naive Bayes.

Ensemble model results

The ensemble strategies Weighted Voting (AUC = 0.7513) and Soft Voting (AUC = 0.7513) outperformed the best individual model LightGBM (AUC = 0.7507) in internal validation. RotF (AUC = 0.7454), GBDT (AUC = 0.7367), Stacking ensemble (AUC = 0.7321), and RF (AUC = 0.7252) also demonstrated strong performance. Regarding Brier scores, LightGBM (Brier = 0.1856) achieved the lowest Brier score within the ensemble framework, while the ensemble models held a slight advantage in discrimination (Table 2 and Figure 9).

Figure 9 Internal performance of ensemble models. A: Receiver operating characteristic curves; B: Precision-recall curves; C: Calibration curves; D: Decision curve analysis; E: Brier scores. PPV: Positive predictive value; ROC: Receiver operating characteristic; GBDT: Gradient Boosting Decision Tree; RotF: Rotation Forest; RF: Random Forest; TPR: True positive rate; FPR: False positive rate.

Table 2 Internal validation performance comparison of ensemble models.

Algorithm	Type	AUC	AUC_95CI	AUPRC	Brier	ABscore	Accuracy	Sensitivity	Specificity	F1	PPV	NPV	Opt_Threshold
WeightedVoting	Ensemble	0.7513	0.6884-0.8149	0.5508	0.1862	0.9632	0.7303	0.7093	0.7403	0.6289	0.5648	0.8428	0.3507
SoftVoting	Ensemble	0.7513	0.6882-0.8149	0.5506	0.1862	0.9600	0.7303	0.7093	0.7403	0.6289	0.5648	0.8428	0.3512
LightGBM	Single	0.7507	0.6872-0.8157	0.5592	0.1856	0.9889	0.7191	0.7907	0.6851	0.6445	0.544	0.8732	0.3144
RotF	Single	0.7454	0.6783-0.8091	0.5333	0.1875	0.769	0.7266	0.7326	0.7238	0.6332	0.5575	0.8506	0.3459
GBDT	Single	0.7367	0.6723-0.7967	0.5439	0.1876	0.5944	0.6704	0.7791	0.6188	0.6036	0.4926	0.855	0.2922
Stacking	Ensemble	0.7321	0.6701-0.7964	0.5317	0.1897	0.3698	0.6929	0.7093	0.6851	0.598	0.5169	0.8322	0.3221
RF	Single	0.7252	0.6620-0.7862	0.5171	0.1935	0	0.6779	0.8023	0.6188	0.6161	0.5	0.8682	0.3333

AUC: Area under the curve; AUPRC: Area under the precision-recall curve; PPV: Positive predictive value; NPV: Negative predictive value; RotF: Rotation Forest; GBDT: Gradient Boosting Decision Tree; RF: Random Forest.

External validation results

The Stacking ensemble achieved the best performance in external validation (AUC = 0.7127, ABscore = 0.8906), followed by Weighted Voting (AUC = 0.7112, ABscore = 0.9731) and Soft Voting (AUC = 0.7111, ABscore = 0.9713). All three ensemble strategies outperformed all individual models in external validation (Table 3 and Figure 10).

Figure 10  External validation performance. A: Receiver operating characteristic curves; B: Precision-recall curves; C: Calibration curves; D: Decision curve analysis; E: Brier scores. PPV: Positive predictive value; ROC: Receiver operating characteristic; GBDT: Gradient Boosting Decision Tree; RotF: Rotation Forest; RF: Random Forest; TPR: True positive rate; FPR: False positive rate.

Table 3 External validation (Gansu Provincial Cancer Hospital, Center B): Performance comparison across models.

Algorithm	Type	Ext_AUC	Ext_AUC_95CI	Int_AUC	AUPRC	Brier	ABscore	Accuracy	Sensitivity	Specificity	F1	PPV	NPV	Threshold (from internal validation)
Stacking	Ensemble	0.7127	0.6533-0.7722	0.732	0.4476	0.1857	0.8881	0.6955	0.6796	0.702	0.5622	0.4795	0.8443	0.3221
Weighted Voting	Ensemble	0.7112	0.6498-0.7702	0.7498	0.439	0.185	0.9734	0.7039	0.699	0.7059	0.576	0.4898	0.8531	0.3507
Soft Voting	Ensemble	0.7111	0.6495-0.7699	0.7499	0.4389	0.185	0.9712	0.7067	0.699	0.7098	0.5783	0.4932	0.8538	0.3512
LightGBM	Single	0.7087	0.6492-0.7664	0.7507	0.4331	0.1863	0.725	0.676	0.7087	0.6627	0.5573	0.4591	0.8492	0.3144
GBDT	Single	0.7075	0.6473-0.7699	0.7367	0.4452	0.1874	0.528	0.648	0.7476	0.6078	0.55	0.435	0.8564	0.2922
RF	Single	0.6961	0.6352-0.7541	0.7252	0.4168	0.1882	0.204	0.6369	0.7476	0.5922	0.5423	0.4254	0.8531	0.3333
RotF	Single	0.6848	0.6228-0.7436	0.7454	0.4384	0.188	0.037	0.6732	0.6893	0.6667	0.5483	0.4551	0.8416	0.3459

AUC: Area under the curve; AUPRC: Area under the precision-recall curve; PPV: Positive predictive value; NPV: Negative predictive value; RotF: Rotation Forest; GBDT: Gradient Boosting Decision Tree; RF: Random Forest.

Regarding Brier scores, Weighted Voting (Brier = 0.1850) and Soft Voting (Brier = 0.1850) were jointly optimal, with Stacking (Brier = 0.1857) closely following. All candidate models achieved Brier scores below 0.19 on the external validation set (Figure 10E).

The selection of the final optimal model followed a systematic multi-criteria decision framework prioritizing: (1) Balanced internal and external discrimination; (2) Calibration accuracy; and (3) Model parsimony. In the external validation cohort, Stacking achieved marginally higher AUC (0.7127 vs 0.7112, a difference of 0.0015), whereas Weighted Voting achieved the best Brier score (0.1850 vs 0.1857), the calibration slope closest to the ideal value, and the smallest calibration intercept deviation, indicating that its predicted probabilities more closely approximated actual event rates across the full risk spectrum. The AUC difference of 0.0015 was clinically negligible and fell within the overlap of the two models’ 95%CIs. We prioritized calibration over this marginal AUC gain for two reasons: First, well-calibrated predicted probabilities are essential for patient communication and shared clinical decision-making in real-world deployment; second, Weighted Voting’s parameter-free architecture (using pre-determined AUC-based weights without additional model fitting) is inherently more robust to distribution shifts in new clinical settings than Stacking’s two-layer meta-learning architecture. On this basis, the Weighted Voting ensemble model was identified as the optimal prediction model (Table 3 and Figure 10).

SHAP interpretability analysis

Global feature importance: Based on the SHAP analysis of the Weighted Voting model, the top 10 features ranked by mean absolute SHAP value (mean |SHAP value|) were, in descending order: LCR, BMI, ApoB/ApoA1, CAR, TyG-BMI, FHR, FAR, LMR, PLR, and age (Figure 11 and Supplementary Table 6).

Figure 11  SHapley Additive exPlanations global analysis. A: Feature importance bar chart and Nightingale rose chart. The bar chart ranks features in descending order by mean absolute SHapley Additive exPlanations value; bar length represents global importance. The Nightingale rose chart displays the relative contribution proportions of each feature; B: Beeswarm plot showing the directional impact and distribution of each feature on model predictions. LCR: Lymphocyte-to-C-reactive protein ratio; BMI: Body mass index; CAR: C-reactive protein-to-albumin ratio; FHR: Fibrinogen-to-high-density lipoprotein cholesterol ratio; ApoB/ApoA1: Apolipoprotein B-to-apolipoprotein A1 ratio; SHAP: SHapley Additive exPlanations; LMR: Lymphocyte-to-monocyte ratio; FAR: Fibrinogen-to-albumin ratio; PLR: Platelet-to-lymphocyte ratio; TyG-BMI: Triglyceride-glucose-body mass index.

LCR ranked first by a substantial margin, underscoring the central role of inflammatory-immune balance in CRC pulmonary metastasis prediction. This finding is consistent with recent reports on the independent predictive value of LCR in solid tumor prognosis[10]. BMI ranked second, and combined with TyG-BMI ranking fifth, this suggests that metabolic status (obesity/insulin resistance) represents the second most important predictive dimension after inflammation. The lipid metabolism index ApoB/ApoA1 ranked third, indicating that pro-atherogenic lipoprotein profile dysregulation is a significant predictor independent of traditional inflammatory indices.

Global interpretability of feature effects: Dependence plots showed the relationship between SHAP values and original feature values for the top 6 features, revealing nonlinear associations (Figure 12).

Figure 12  SHapley Additive exPlanations dependence plots. A: Lymphocyte-to-C-reactive protein ratio; B: Body mass index; C: Apolipoprotein B-to-apolipoprotein A1 ratio; D: C-reactive protein-to-albumin ratio; E: Triglyceride-glucose-body mass index; F: Fibrinogen-to-high-density lipoprotein cholesterol ratio. LCR: Lymphocyte-to-C-reactive protein ratio; BMI: Body mass index; ApoB/ApoA1: Apolipoprotein B-to-apolipoprotein A1 ratio; SHAP: SHapley Additive exPlanations; CAR: C-reactive protein-to-albumin ratio; FHR: Fibrinogen-to-high-density lipoprotein cholesterol ratio; TyG-BMI: Triglyceride-glucose-body mass index.

The LCR dependence plot (Figure 12A) exhibited a prominent nonlinear effect. At lower LCR levels, SHAP values were predominantly positive and of large magnitude, suggesting that low LCR was associated with high model-predicted metastasis risk; as LCR increased, SHAP values gradually shifted to negative territory.

The BMI dependence plot (Figure 12B) revealed a non-monotonic relationship between BMI and SHAP values. When BMI ranged between 20 and 25, SHAP values were less than 0, whereas SHAP values rose markedly when BMI exceeded 25 - an observation directionally consistent with existing biological evidence linking obesity to tumor metastasis. The ApoB/ApoA1 dependence plot (Figure 12C) displayed a relatively clear positive correlation trend, with SHAP values rising notably when ApoB/ApoA1 exceeded 0.9, suggesting a largely linear association between severe pro-atherogenic lipoprotein imbalance and metastatic risk.

In the CAR dependence plot (Figure 12D), SHAP values remained relatively flat when CAR was below 0.5 but shifted sharply positive when CAR exceeded 1.0, indicating that CAR becomes a strong predictor of metastatic risk only at higher inflammatory burdens. The TyG-BMI dependence plot (Figure 12E) showed that SHAP values began to shift positively when TyG-BMI exceeded 200, consistent with the pro-metastatic effects of insulin resistance. The FHR dependence plot (Figure 12F) showed that SHAP values increased when FHR exceeded 3, suggesting that imbalance along the fibrinogen-HDL axis beyond this threshold contributes substantially to metastatic risk.

It must be emphasized that the visual inflection points observed in the above dependence plots reflect nonlinear trends in SHAP values across the feature space, with their numerical values jointly influenced by the data distribution of the development cohort, model architecture, and the SHAP approximation algorithm. These inflection points may serve as exploratory leads for generating clinical hypotheses but should not be directly applied to clinical risk stratification decisions without validation in independent cohorts. Future studies may employ subgroup comparisons in external cohorts or recursive partitioning methods to evaluate the clinical reliability of these thresholds.

Individual-level local explanations: To demonstrate the decomposition of model predictions at the individual patient level, four representative cases were selected (Figure 13). For each case, waterfall plots and force plots were used to display the positive (red; driving metastasis risk upward) and negative (blue; driving metastasis risk downward) SHAP contribution values of each feature. In the waterfall plots, E[f(X)] at the bottom represents the model’s baseline predicted value, with features sequentially stacked to yield the patient’s final predicted value f(x). Force plots present all feature contributions along a single horizontal axis: Red arrows represent features pushing the prediction toward metastasis, while blue arrows represent features pushing the prediction away from metastasis; the width of each arrow reflects the absolute magnitude of that feature’s SHAP contribution.

Figure 13  SHapley Additive exPlanations local model explanations. A: Case 1; B: Case 2; C: Case 3; D: Case 4. BMI: Body mass index; ApoB/ApoA1: Apolipoprotein B-to-apolipoprotein A1 ratio; LCR: Lymphocyte-to-C-reactive protein ratio; CAR: C-reactive protein-to-albumin ratio; FHR: Fibrinogen-to-high-density lipoprotein cholesterol ratio; TyG-BMI: Triglyceride-glucose-body mass index; LMR: Lymphocyte-to-monocyte ratio; FAR: Fibrinogen-to-albumin ratio; PLR: Platelet-to-lymphocyte ratio.

Case 1 (Figure 13A): A high-risk true-positive patient (predicted probability P = 0.748; actual metastasis = 1). The final model output was f(x) = 1.082. The waterfall plot indicated that BMI and ApoB/ApoA1 were the two primary drivers of the elevated risk prediction, followed by TyG-BMI and LCR.

Case 2 (Figure 13B): A high-risk true-positive patient (P = 0.712; actual metastasis = 1). The largest SHAP contribution came from ApoB/ApoA1 rather than BMI, with FHR and PLR also contributing positive risk, suggesting that the coagulation-inflammation axis and platelet-mediated immune evasion played roles in this patient’s metastatic risk.

Case 3 (Figure 13C): A low-risk true-negative patient (P = 0.095; actual metastasis = 0), with model output f(x) = -2.194. In stark contrast to the preceding cases, SHAP values for nearly all features were negative. LCR contributed the largest protective effect, followed by CAR and FAR. This patient’s LCR was 2.8, indicating robust lymphocyte immune function and very low systemic inflammation; both BMI and TyG-BMI were at low levels, suggesting the absence of obesity and insulin resistance.

Case 4 (Figure 13D): This patient presented with concurrent risk-increasing factors (FHR = 3.66, LMR = 1.86, CAR = 0.08) and risk-decreasing factors (BMI = 20.76, PLR = 271.0, ApoB/ApoA1 = 0.48, age = 51).

Feature interaction effects: The SHAP interaction value heatmap (Figure 14A) revealed the most prominent second-order feature interaction effects. The strongest interaction pairs included: BMI × LCR, BMI × age, TyG-BMI × LCR, and ApoB/ApoA1 × age.

Figure 14  Feature interaction effects. A: SHapley Additive exPlanations feature interaction value heatmap (Weighted Voting); B: Three-dimensional partial dependence plots. BMI: Body mass index; SIRI: Systemic inflammation response index; PLR: Platelet-to-lymphocyte ratio; LCR: Lymphocyte-to-C-reactive protein ratio; CAR: C-reactive protein-to-albumin ratio; LMR: Lymphocyte-to-monocyte ratio; FAR: Fibrinogen-to-albumin ratio; mGPS: Modified Glasgow Prognostic Score; PNI: Prognostic nutritional index; ApoB/ApoA1: Apolipoprotein B-to-apolipoprotein A1 ratio; FHR: Fibrinogen-to-high-density lipoprotein cholesterol ratio; TyG-BMI: Triglyceride-glucose-body mass index.

The strong interaction effect between BMI and LCR carries important biological significance. Obesity (high BMI) exacerbates systemic inflammation and suppresses lymphocyte function through pro-inflammatory adipokines secreted by adipose tissue, thereby forming a positive feedback loop with decreased LCR that synergistically promotes pre-metastatic niche formation[2]. The interaction effect between TyG-BMI and LCR further supports the hypothesis of insulin resistance-inflammation crosstalk playing a synergistic role in CRC metastasis. Three-dimensional partial dependence plots visualized their joint effects (Figure 14B).

Clinical application tool

An interactive web application was constructed based on the Weighted Voting model and all 16 features. Upon entering the actual values of 10 key clinical indicators for a given patient, the application automatically calculates and displays the patient’s predicted probability of pulmonary metastasis in a semicircular gauge format (Figure 15). The application is accessible online at https://cccrrrccclllmmm.streamlit.app/.

Figure 15  Clinical application tool. CRC: Colorectal cancer; BMI: Body mass index; LCR: Lymphocyte-to-C-reactive protein ratio; CAR: C-reactive protein-to-albumin ratio; ApoB/ApoA1: Apolipoprotein B-to-apolipoprotein A1 ratio; FAR: Fibrinogen-to-albumin ratio; PLR: Platelet-to-lymphocyte ratio; FHR: Fibrinogen-to-high-density lipoprotein cholesterol ratio; LMR: Lymphocyte-to-monocyte ratio; TyG-BMI: Triglyceride-glucose-body mass index.

DISCUSSION

This study systematically developed and validated a machine learning prediction model based on multimodal inflammation-lipid metabolism biomarkers using real-world multicenter retrospective cohort data, for assessing the risk of pulmonary metastasis in postoperative patients with CRC. Pulmonary metastasis represents one of the primary patterns of distant dissemination in CRC; once established, it typically signifies disease progression into a more complex phase of systemic dissemination, necessitating comprehensive adjustments to surveillance strategies, selection of locoregional and systemic therapeutic approaches, and prognostic assessment. Consequently, the early identification of individuals at high risk of pulmonary metastasis through interpretable and generalizable prediction models carries substantial clinical significance for achieving more refined risk stratification and more efficient allocation of healthcare resources.

From a study design perspective, the present study offers several notable strengths. First, the multicenter design more faithfully simulates the distribution shift encountered in real-world deployment settings compared with single-center studies, thereby providing a more rigorous assessment of cross-institutional generalizability. Second, we proposed a four-tier progressive feature combination strategy that systematically quantified the incremental predictive contributions of demographic characteristics, systemic inflammatory biomarkers, and lipid metabolism biomarkers to CRC pulmonary metastasis prediction - an approach that is relatively uncommon in prior studies of this nature. Furthermore, the dataset encompassed 16 candidate variables spanning three dimensions: Demographic characteristics, systemic inflammatory indices, and lipid metabolism indices. All variables are obtainable through routine clinical examination without the need for specialized molecular testing platforms, thereby laying the groundwork for model deployment in primary care settings.

With respect to modeling strategy, this study systematically compared 14 machine learning algorithms spanning diverse learning paradigms and further explored three ensemble strategies (Soft Voting, Weighted Voting, and Stacking). This comprehensive algorithmic comparison framework facilitates identification of the learning strategy most suited to the current data structure. Among the 14 algorithms evaluated, LightGBM demonstrated superior performance in internal validation under Layer 4 (full-feature model), achieving an AUC of 0.7507 (95%CI: 0.6872-0.8157), the lowest Brier score (0.1856), and the highest ABscore (1.0000), reflecting the most balanced overall ranking across discrimination and calibration metrics. The superior performance of LightGBM can be attributed to its histogram-based approximate split algorithm and leaf-wise growth strategy[47], which enable effective capture of nonlinear feature interactions in moderately sized datasets containing multiple correlated features, while controlling overfitting through built-in L1/L2 regularization mechanisms.

Regarding ensemble modeling, Weighted Voting achieved the highest performance in internal validation (AUC = 0.7513, 95%CI: 0.6884-0.8149) and also demonstrated the best overall performance in external validation. The AUC degradation of all candidate models in external validation was less than 0.07, falling within the acceptable range recommended by Debray et al[48] (ΔAUC < 0.10), indicating that the models developed in this study possess reliable cross-institutional generalizability.

An important consideration for interpreting the external validation results is the heterogeneity between the two participating centers. Guang’anmen Hospital (Center A) is a national tertiary hospital in Beijing, whereas Gansu Provincial Cancer Hospital (Center B) is a regional cancer center in northwestern China, serving a geographically, socioeconomically, and potentially demographically distinct population. The two centers utilized different laboratory analytical platforms and reagent systems, and center-specific variations in surgical techniques, adjuvant chemotherapy regimens, and follow-up imaging protocols may exist. Notably, because all predictor variables were collected preoperatively (within 14 days before curative surgery), post-surgical treatment protocol differences did not affect predictor values, although they could influence outcome ascertainment through differences in surveillance intensity. Despite these sources of heterogeneity, the model demonstrated consistent external performance with an AUC degradation of only 0.04 and a Brier score of 0.1850, suggesting that the learned predictor-outcome associations are robust to the inter-center variation present in this study. Nevertheless, the two-center design cannot fully characterize the spectrum of clinical heterogeneity encountered in broader practice, and validation across a greater number of geographically and institutionally diverse centers remains necessary.

The robust performance of Weighted Voting can be attributed to three mechanisms. First, performance-proportional weighting concentrates ensemble predictions around the most discriminative base learner (LightGBM) while retaining complementary information from architecturally diverse models (RotF, GBDT, RF), whose partially decorrelated prediction errors reduce overall variance upon averaging. Second, unlike Stacking - which introduces a second-layer meta-learner with additional model-fitting steps that risk overfitting to development cohort-specific patterns - Weighted Voting employs pre-determined AUC-based weights without additional parameter estimation, thereby preserving calibration accuracy during cross-center transfer. This explains the observed pattern whereby Stacking achieved marginally higher external AUC (0.7127 vs 0.7112) but inferior calibration (Brier = 0.1857 vs 0.1850). Third, the simplicity and transparency of the weighting scheme facilitate clinical deployment and regulatory review.

SHAP interpretability analysis revealed the central role of LCR in predicting CRC pulmonary metastasis. LCR integrates both host immune surveillance function and systemic inflammatory burden. SHAP analysis demonstrated a strong association between low LCR values and increased pulmonary metastasis risk, with marked interindividual heterogeneity. This finding is consistent with the report by Okugawa et al[10] in a CRC cohort, which identified LCR as an independent predictor of surgical and oncological outcomes. From a biological mechanism perspective, lymphocytopenia reflects a tumor-induced immunosuppressive microenvironment characterized by T-cell exhaustion, regulatory T-cell (Treg) expansion, and upregulation of immune checkpoint pathways[49]; conversely, elevated CRP signifies activation of a systemic inflammatory response driven by tumor cells and tumor-associated macrophages, which promote pre-metastatic niche formation via the IL-6/STAT3 signaling pathway[50]. The bidirectional imbalance - diminished immune surveillance coupled with exaggerated inflammatory response - collectively creates favorable conditions for the colonization and survival of circulating tumor cells within the pulmonary parenchyma.

From a clinical application perspective, the SHAP-derived inflection points - particularly LCR < 0.5, BMI > 25, and ApoB/ApoA1 > 0.9 - offer preliminary signals that may inform postoperative surveillance intensity. For instance, a patient presenting with concurrent low LCR, elevated BMI, and high ApoB/ApoA1 ratio could be flagged for more frequent chest CT follow-up. However, these thresholds are exploratory and require prospective validation before clinical adoption. The interactive web application we developed enables individualized risk assessment by integrating all features simultaneously, transcending single-feature threshold-based approaches.

BMI ranked second in SHAP feature importance, with TyG-BMI ranking fifth, suggesting that metabolic status constitutes the second most important predictive dimension after inflammation. The mechanisms by which obesity promotes CRC metastasis involve multiple interconnected biological pathways: (1) Pro-inflammatory adipokines secreted by adipose tissue promote tumor cell proliferation and invasion via the AMPK/mTOR signaling pathway; (2) Obesity-related insulin resistance activates the PI3K/Akt/mTOR pathway through hyperinsulinemia and insulin-like growth factor-1 (IGF-1), enhancing tumor cell proliferation, survival, and migration; and (3) The chronic low-grade inflammatory state in obese patients acts synergistically with tumor microenvironment inflammation, further promoting pre-metastatic niche formation[2,51]. SHAP interaction value analysis confirmed that BMI × LCR represents the strongest second-order interaction pair, providing data-driven support for the hypothesis that obesity and inflammation act synergistically to facilitate metastasis.

ApoB/ApoA1 ranked third in SHAP feature importance, indicating that an atherogenic lipoprotein profile imbalance constitutes an important predictor independent of traditional inflammatory markers. Elevated ApoB/ApoA1 reflects the excess of LDL/VLDL relative to HDL, suggesting a shift in lipid metabolism toward a pro-atherogenic and pro-tumorigenic phenotype. Prior research has demonstrated that ApoB promotes oxidized LDL generation and endoplasmic reticulum stress responses, indirectly facilitating tumor metastasis through activation of NF-κB inflammatory signaling[15]; meanwhile, ApoA1 has been shown to possess anti-inflammatory, antioxidant, and antitumor properties[52]. Accordingly, an elevated ApoB/ApoA1 ratio may simultaneously reflect both enhanced pro-tumorigenic lipid signaling and attenuated antitumor protective mechanisms. FHR ranked sixth, bridging the coagulation-inflammation axis with the lipid metabolism axis, highlighting the role of concurrent dyslipidemia and hypercoagulability in CRC pulmonary metastasis.

Compared with prior studies on CRC pulmonary metastasis prediction, the present study offers several notable contributions. First, regarding predictor selection, to our knowledge, this study is among the first to systematically integrate systemic inflammatory indices and lipid metabolism indices as dual predictive dimensions for constructing a CRC pulmonary metastasis prediction model. Previous studies have predominantly focused on traditional clinicopathological variables[4] or single-dimension biomarkers, whereas the stratified modeling approach employed here suggests that inflammatory and lipid metabolism biomarkers may provide incremental predictive contributions, offering a new biomarker dimension for CRC pulmonary metastasis risk assessment. Second, regarding modeling methodology, this study systematically compared 14 machine learning algorithms and 3 ensemble strategies, a substantially more comprehensive evaluation than prior studies that typically assessed only 2-5 algorithms. Third, regarding model interpretability, SHAP analysis was employed to elucidate the internal logic underlying model decisions at multiple levels, overcoming the limitations of the black-box nature of conventional machine learning. Fourth, the final model was deployed as an interactive web application, enabling translation from a research tool to a clinical instrument, allowing frontline clinicians to conveniently access individualized risk assessments and interpretable outputs.

In contextual benchmarking, prior CRC pulmonary metastasis prediction models - predominantly nomograms incorporating pathological variables such as TNM staging, histological grade, and lymph node status - have reported C-indices ranging from 0.70 to 0.78 in their respective validation cohorts[18,21]. Our model’s AUC of 0.7507 (internal) and 0.7112 (external) falls within this performance range despite intentionally excluding all pathological variables. This comparison underscores that routinely obtainable inflammation-lipid metabolism biomarkers carry substantial predictive information for CRC pulmonary metastasis. Moreover, unlike many prior models that lacked independent external validation, our model demonstrated consistent cross-center generalizability with an AUC degradation of only 0.04, within the acceptable criterion (ΔAUC < 0.10) recommended by Debray et al[48]. Nevertheless, the model is positioned as a complementary - not replacement - tool for use before pathological results become available or in settings where pathological evaluation is delayed.

An important consideration is whether the observed AUC of approximately 0.75 is sufficient to support clinical decision-making. For a screening or triage tool - where the objective is to identify high-risk patients for intensified surveillance rather than to make definitive therapeutic decisions - moderate discrimination can provide meaningful clinical utility, provided that the model demonstrates positive net benefit. In our study, decision curve analysis confirmed positive net benefit across a clinically relevant threshold probability range of approximately 10%-40%, indicating that model-guided surveillance decisions would yield greater benefit than either universal intensive follow-up or uniform standard follow-up. In practical terms, a patient flagged as high-risk by the model could receive more frequent chest CT surveillance (e.g., every 6 months rather than annually), a low-cost, low-risk intervention with potential for earlier detection of pulmonary metastasis. Furthermore, achieving this level of discrimination using only routinely obtainable preoperative biomarkers - in the absence of pathological staging information - is notable, as it provides risk stratification in a clinical window where no alternative quantitative tool currently exists. Nevertheless, the model should be regarded as a preliminary screening tool rather than a definitive diagnostic instrument, and integration with pathological variables, radiomics features, or molecular biomarkers in future studies may further enhance predictive accuracy.

Although the prediction models developed in this study demonstrated satisfactory discrimination and calibration in both internal and external validation, several limitations warrant acknowledgment. First, the retrospective cohort study design, despite the use of multicenter external validation to assess model generalizability, cannot fully eliminate the inherent selection bias and information bias associated with retrospective data collection. Future prospective, multicenter cohort studies - enrolling consecutive patients at geographically diverse centers with prespecified performance endpoints (AUC, calibration slope, net benefit) and a fixed classification threshold derived from the current study - are warranted to further evaluate model predictive performance and temporal stability. Such studies could additionally assess whether model-guided surveillance strategies translate into earlier detection of pulmonary metastasis and improved patient outcomes through randomized or cluster-randomized comparisons. Relatedly, although the two participating centers exhibited meaningful differences in geographic location, patient demographics, laboratory platforms, and likely treatment protocols, the two-center design provides only a limited assessment of the model’s generalizability across the full spectrum of clinical practice heterogeneity. Differences in follow-up imaging frequency between centers may have introduced ascertainment bias for pulmonary metastasis detection. Future multi-center validation should include centers with diverse follow-up protocols and treatment paradigms to more rigorously characterize the model’s transportability. Second, the optimal model achieved an internal validation AUC of approximately 0.75, reflecting moderate predictive performance. Although all models demonstrated positive net benefit in decision curve analysis, predictive accuracy remains below that of certain studies incorporating radiomics features or molecular biomarkers, leaving room for improvement. Future investigations may consider integrating radiomics features and molecular biomarkers to further enhance predictive performance. Third, the study population was drawn from two medical centers in China, and its demographic and clinical characteristics - including BMI distribution, lipid metabolism profiles, and baseline inflammatory states - may differ from patients in other geographic regions or ethnic backgrounds, potentially limiting the external generalizability of the model to diverse populations. Future cross-population validation should encompass multi-ethnic cohorts from geographically diverse regions, with subgroup transportability analyses stratified by ethnicity and metabolic phenotype. Where discrimination is preserved but calibration degrades, established recalibration techniques may be applied to adapt the model without complete re-derivation. Federated learning approaches could further enable privacy-preserving international multicenter validation. Fourth, although SHAP analysis revealed statistical associations and effect directions between individual features and pulmonary metastasis prediction, these associations should not be interpreted as causal relationships.

CONCLUSION

This study developed and externally validated a machine learning framework integrating inflammation and lipid metabolism biomarkers for predicting postoperative pulmonary metastasis in CRC. Through a hierarchical four-layer modeling strategy encompassing 14 algorithms and three ensemble approaches, the Weighted Voting ensemble incorporating all 16 candidate features demonstrated the best predictive performance, with an AUC of 0.7513 in internal validation and 0.7112 in external validation, confirming acceptable cross-center generalizability. SHAP-based interpretability analysis identified LCR, BMI, and ApoB/ApoA1 ratio as the most influential predictors, with clinically actionable thresholds (LCR < 0.5, BMI > 25, ApoB/ApoA1 > 0.9) that align with established biological mechanisms linking systemic inflammation and lipid metabolic dysregulation to metastatic progression. Notably, all predictor variables are derived from routine preoperative laboratory tests, facilitating seamless clinical integration without additional cost. An interactive web-based application was further developed to enable individualized risk stratification at the point of care. Future prospective, multicenter studies are warranted to validate these findings across diverse populations and to evaluate the clinical impact of model-guided surveillance strategies on patient outcomes.

ACKNOWLEDGEMENTS

We would like to express our gratitude to each patient who participated in this study.