Predicting Lymph Node Metastasis in T1 Colorectal Cancer Patients Using Interpretable Machine Learning Models: A Multicenter Retrospective Study

Globally, CRC is recognized as the third most frequently diagnosed malignancy and ranks second in terms of cancer-related mortality [1]. With widespread CRC screening, diagnosed cases of T1 CRC have increased [2, 3]. In China, CRC has become the second most common malignant tumor (29.51 per 100,000) and the fourth primary cause of cancer-associated fatalities (14.14 per 100,000) [4]. LNM, defined as cancer cells found in lymph nodes surrounding the primary tumor, indicates tumor spread [5]. Clinically, the presence of LNM dictates treatment: patients without LNM may achieve curative results through endoscopic resection, whereas those with LNM require radical surgery [6]. Currently, LNM evaluation mainly relies on MRI, CT, and biopsy [7-9]. However, these methods fail to detect LNM in 8-12% of patients due to imaging limitations [10], and their high cost and long duration limit widespread use. Hence, there is an urgent need for more efficient, imaging-independent detection methods. Compared to traditional logistic regression (LR), recent advancements in statistical theory and machine learning have led to predictive models with superior performance, enhancing prediction of LNM in cancers such as gastric, thyroid, and breast [11-13].

The history of machine learning dates to 1957, when Rosenblatt [14] proposed the "Perceptron", recognized as the earliest machine learning model. Simultaneously, LR, developed between the late 19th and early 20th centuries, gained prominence for binary classification tasks [15]. LR remains a common baseline method, though it is often outperformed by more complex models like Support Vector Machines (SVMs) and Neural Networks (NNs). Linear models such as logistic and Least Absolute Shrinkage and Selection Operator (LASSO) regression have limitations. Hartwig et al. [16] employed a LASSO model using data from 35,812 Danish patients, achieving a relatively low AUC (0.64), indicating poor model selection. Wang et al. [17] used LR on 825 T1-stage CRC patients, achieving an AUC of 0.793, yet faced overfitting concerns due to limited data. Fujino et al. [18] assessed LNM risks with LR on 934 Japanese patients, yielding an AUC of 0.786 (training) and 0.721 (validation), but lacked robustness due to the absence of cross-validation. Similarly, Niu and Cao [19] achieved an AUC of only 0.708 from a large SEER database without cross-validation, limiting its generalizability.

SVM is a machine learning method using optimization to avoid the "curse of dimensionality" and "overfitting" by finding an optimal hyperplane that maximally separates classes. Ichimasa et al. [20] analyzed data from 690 T1-stage CRC patients (590 training, 100 validation) between 2001 and 2016. Their model incorporated 45 clinicopathological factors—such as age, Body Mass Index (BMI), tumor size, location, LVI, tumor markers (CEA, CA19-9), and biochemical indicators—to predict LNM risk. ROC analysis showed good accuracy (AUC=0.821). However, small sample size, low LNM prevalence (55 positive cases in training, 9 in validation), data imbalance, and inclusion of clinically ambiguous features limited model reliability and practical applicability.

The training mechanism of NNs comprises two essential phases: forward propagation of signals and back-propagation of errors. During the forward pass, input data traverse through multiple hidden layers before reaching the output layer. If actual outputs differ from expected, error signals propagate backward for weight adjustments. NNs evolved into models like TabNet, an NN combining attention mechanisms and tree-based feature selection. TabNet’s main steps include (a) sparse feature selection; (b) sequential multi-step structures for incremental decision contributions; (c) nonlinear processing to enhance learning; and (d) ensemble processing through multiple steps. Song et al. [21] developed a deep learning model leveraging attention mechanisms to analyze whole-slide images (WSIs) for predicting LNM in T1 CRC, attaining AUC values ranging from 0.781 to 0.824. Despite promising results, limitations include single-center WSI data, small samples, potential preparation biases, complexity in interpreting high-dimensional data, and challenging feature associations, complicating model interpretability.

Decision Trees (DTs) are flowchart-like structures with nodes representing attribute tests, branches representing outcomes, and leaves representing class distributions. AdaBoost sequentially trains a series of weak learners and aggregates them into a composite strong model. Random Forest (RF) integrates multiple DTs via bagging, while Gradient Boosted Decision Trees (GBDTs) iteratively aggregate weak learners for improved predictions. Takamatsu et al. [22] proposed a hybrid approach, combining Convolutional Neural Networks (CNNs) with RF classifiers to predict LNM from histopathological WSIs in T1 CRC cases. Among 783 cases divided into training (548) and validation (235), AUC values were 0.971 (training) and 0.760 (validation). However, training the CNN on roughly 500 cases led to an AUC drop of 21%, indicating significant overfitting and limited reliability.

XGBoost, proposed by Chen [23], is a powerful tree-based ensemble learning algorithm widely adopted for classification and regression problems. Unlike RF and GBDT, XGBoost applies first- and second-order derivatives via a second-order Taylor expansion of the loss function, introducing L1 and L2 regularization and column sampling to reduce overfitting. XGBoost represents trees using leaf weights, directly expressing the loss function's extremum through derivatives [24]. The "gain" metric, reflecting the relative contribution of features, determines optimal split points. Ahn et al. [25] applied XGBoost to data from 26,733 T1-stage CRC patients (SEER database), using eight prognostic variables and achieving a modest AUC of 0.659 (sensitivity 0.242, accuracy 0.604, F1 score 0.352). Despite the large sample and relevant variables, the low performance suggests inadequate parameter tuning, limiting the model’s reliability.

The LightGBM model, a highly optimized implementation of the GBDT algorithm, is designed to handle large-scale data with superior computational efficiency. It supports parallel learning, operates at high speed, consumes less memory, and delivers improved predictive accuracy. In a study conducted by Piao et al. [26], a LightGBM-driven model was constructed using data from 651 patients diagnosed with T1-stage CRC who underwent radical resection across six medical institutions in Ningbo between 2016 and 2022. The model demonstrated excellent discriminative ability, with an AUC of 0.960. Despite this high performance, the small sample size raises concerns about overfitting. Additionally, reliance on pathologist-dependent features, such as submucosal invasion area, limits the model’s standardization and broader applicability.

CatBoost, a gradient boosting algorithm released as an open-source project by the Russian tech company Yandex in 2017, represents another advancement in the GBDT family. While XGBoost is widely used and LightGBM improves computational efficiency, CatBoost claims superior accuracy [27]. It employs symmetric DTs (oblivious trees) as base learners, featuring fewer parameters and built-in support for categorical variables. This enables efficient handling of categorical data—reflected in its name, combining "Categorical" and "Boosting." Furthermore, CatBoost mitigates issues such as gradient bias and prediction shift, thereby reducing the likelihood of overfitting and enhancing both accuracy and model generalizability.

Machine learning methods are advantageous in capturing higher-order nonlinear interactions among predictive factors, often leading to more stable and robust predictions. However, their inherent "black box" nature makes it difficult to explain the rationale behind specific predictions, limiting their acceptance in clinical decision-making. SHAP addresses this challenge by precisely quantifying the contribution and influence of each input feature on the model's output [28]. While numerous machine learning models have shown high performance, a recurring shortcoming in previous studies is the absence of transparent and intuitive feature attribution. Without clear attribution, clinicians may struggle to identify key features, limiting the model’s practical utility in guiding treatment decisions.

Accurate prediction of LNM remains a clinical challenge for patients with T1 CRC. Although the incidence of LNM is relatively low (approximately 10.9%), many patients still undergo radical surgery, which may be unnecessary in the absence of LNM [29]. Endoscopic resection has become an alternative treatment option for T1 CRC patients without LNM [30]. However, patients with T1 CRC who have concomitant LNM generally have a poorer prognosis. Therefore, establishing an accurate, reliable, and reproducible method for predicting LNM risk is crucial for optimizing preoperative treatment strategies and reducing overtreatment [31-33]. In this study, we combined interpretable machine learning techniques with common demographic and pathological data to develop a model that predicts the LNM risk in patients with T1 CRC, providing an innovative approach for personalized treatment decisions.

T1 CRC refers to tumors confined to the mucosal or submucosal layers without reaching or penetrating the muscularis propria [7, 34]. Currently, the following five high-risk factors are generally used to guide clinical practice: tumor invasion depth, LVI, tumor differentiation grade, tumor budding, and incomplete or positive resection margins [7]. Tumor invasion depth is one of the most important high-risk factors: the deeper the invasion into the submucosa, the higher the likelihood of cancer cells invading lymph nodes. It is generally considered that tumor invasion exceeding 1,000 micrometers (1ămm) is a high-risk marker for LNM [31]. LVI refers to cancer cells invading lymphatic or blood vessels, indicating that cancer cells have already spread through these vessels and may metastasize to distant sites [35]. Poorly differentiated or undifferentiated cancer cells (including mucinous adenocarcinoma or signet ring cell carcinoma) often have higher aggressiveness and are more likely to metastasize to lymph nodes [10]. Tumor budding refers to isolated or small clusters of cancer cells formed at the tumor margin, which is a marker of tumor aggressiveness and metastatic behavior; high levels of tumor budding are significantly associated with an increased risk of LNM [36]. If pathological examination shows incomplete tumor resection margins (R1 resection) or positive margins, it implies that residual tumor cells may be present, increasing the risk of LNM [37]. If any of the above high-risk factors are found in a patient’s specimen, salvage radical surgery is recommended [7]. Even so, a considerable proportion of patients will undergo unnecessary radical surgery. Researchers have reported that T1 cancer patients with high-risk factors account for 70-80% of all T1 CRC cases. Therefore, considering that only 8-12% of patients have true LNM, more than 60% of patients still receive overtreatment [10, 38]. In the guidelines for CRC treatment by the Japanese Society for Cancer of the Colon and Rectum (JSCCR), when using the predictive model in the guidelines to assess LNM in T1 CRC patients, the AUC value is about 0.588 [7]. This indicates that in practical applications, the model has certain reliability in predicting LNM in T1 CRC patients but also has significant room for improvement. Therefore, we need to develop a more accurate and feasible algorithmic model to assist clinicians in formulating surgical plans.

Our study shows that the CatBoost model achieves the highest AUC, F1 Score, Recall, Precision, and Accuracy, demonstrating the best predictive performance among the 29 models compared. Traditional machine learning models, such as artificial Neural Networks (ANNs), SVMs, and RFs, face a fundamental issue known as the "black box problem" [39]. While these models can calculate metastasis probabilities based on patient input data, their internal decision-making processes are difficult to interpret. This makes it challenging for clinicians to clearly understand which demographic characteristics or pathological report details play key roles in the prediction, and they cannot effectively utilize this information to decide on additional treatments. Therefore, the clinical application of these models is somewhat limited, making it difficult to fully meet the medical field’s demand for transparency and interpretability in decision-making.

The CatBoost model effectively addresses this problem. It relies on changes in the internal loss function and pre-diction shifts within the model to evaluate feature importance, offering high computational efficiency and being particularly suitable for tree-based models. However, it primarily provides global interpretations and cannot deeply explain the feature contributions for each individual prediction [27]. In contrast, SHAP, based on Shapley values from cooperative game theory, provides consistency and local interpretability by considering the marginal contributions of features. It is applicable to various model types but has higher computational complexity when handling a large number of features [28]. Therefore, CatBoost is more suitable for scenarios requiring rapid global feature importance, while SHAP offers more comprehensive global and local explanations [40]. To reflect the contributions of all factors in LNM, we used machine learning combined with the SHAP method to evaluate the risk of LNM in T1 CRC patients. This method displays a list of important features, ranked from most to least important (from top to bottom). Consequently, the important features identified by CatBoost were further plotted. The results showed that the model’s predictions are significantly influenced by features such as LVI, age, total number of lymph nodes, tumor size, and depth of tumor invasion, while features like perineural invasion and specific tumor location contributed less to the model. We provided an example to illustrate the model’s interpretability, presenting the prediction results in an intuitive manner that allows clinicians to clearly observe the weights of the included features in the model’s predictions.

In this study, the positive SHAP values indicate that the presence of LVI is associated with an increased likelihood of LNM, highlighting LVI as an important predictor of metastasis risk. The SHAP values for LVI are mainly positive, concentrated on the right side of zero, emphasizing its role in increasing the risk of LNM. Conversely, the lack of negative value distribution further confirms the risk-increasing effect of LVI in predicting LNM. Basic experimental research shows that LVI plays a crucial role in tumor cell invasion. By analyzing lymphatic vessel markers (such as D2-40 and VEGF-C) in tumor tissues, it can be determined that LVI-positive tumors are more prone to LNM [41].

In our study, younger patients were more likely to exhibit LNM. A previous study shows that younger patients have a higher likelihood of LNM due to higher tumor biological activity [42]. This result is consistent with earlier research.

This study demonstrates that tumor size significantly influences the model’s prediction of LNM. Larger tumors correspond to positive SHAP values, indicating an increased risk of metastasis, while smaller tumors correspond to negative SHAP values, suggesting a higher probability of no metastasis. Previous studies have shown that larger tumor sizes are generally associated with a higher risk of LNM, whereas smaller tumors have a lower risk. This may be because larger tumors often have deeper invasion depths and higher aggressiveness, increasing the likelihood of cancer cells metastasizing to lymph nodes [43]. Earlier research has indicated that the deeper the invasion depth of cancer cells, the higher the probability of lymph node invasion [31]. The findings of this study align with existing studies, further confirming that tumor invasion depth is positively correlated with LNM risk. Invasion depth is not only an independent predictor of LNM but has also become a key indicator in CRC staging [44].

Additionally, our study reveals a direct relationship between the total number of lymph nodes and patient LNM. The more lymph nodes are removed during radical surgery, the greater the risk of detecting LNM in postoperative pathology. For radical surgery in T1 CRC, previous studies have shown that the number of lymph nodes obtained during systematic lymphadenectomy can be used to assess the comprehensiveness of the surgeon’s operation. Removing at least 12 lymph nodes allows for more accurate cancer staging, improved prognostic evaluation, and increased long-term survival rates [45].