Clinicopathological data of patients with site recode ICD-O-3/WHO 2008 “pancreas” and AYA site recode 2020 Revision “9.3.9.2 Pancreas – adenocarcinoma” between 2000 and 2021 was retrieved from the SEER database. Additionally, clinicopathological information of pancreatic The First People’s Hospital of Lianyungang (2015–2024) was retrospectively collected through electronic medical record system. The study was conducted according to the guidelines of the Declaration of Helsinki and was approved by the Ethics Committee of The First People’s Hospital of Lianyungang (protocol code: KY-20210910004, approved on 2021-09-10). Informed consent was obtained from all subjects involved in this study. Inclusion criteria comprised individuals with ICD-O-3/WHO 2008 “pancreas” and AYA site recode 2020 Revision “9.3.9.2 Pancreas – adenocarcinoma” which are older than 18 years old. Exclusion criteria comprised patients lacking follow-up information of survival months and death cause, not diagnosed with positive histology, no surgery information, not first malignant tumor, without TNM stage or grade details. In SEER database, metastasis is characterized by spreading to distant organs during the initial cancer diagnosis. And we define the outcome of predictive model as three-year survival, indicating that patients are still alive at the timepoint of 36 months follow-up. The positive outcome was death of patient in three-year follow-up.

Extracted data were gathered on demographic data (age, gender, race, marital status household location and income), cancer characteristics (pathological grade, summary stage, TNM stage, tumor size, tumor primary location, pathology, metastasis information), therapeutic information (surgery, lymph node surgery, positive lymph node, radiotherapy, chemotherapy) and follow-up information (overall and cancer-specific survival status, survival months). Two continuous variables, age and tumor size, were divided into categorical variables. The age was split into five groups: “<50”, “50-59”, “60-69”, “70-79” and “>=80”. The tumor size was split into “<2cm”, “2-3.9cm”, “4-5.9cm”, “6-7.9cm”, “>8cm” and “Unknown”. “Metastasis” was defined as “yes” with metastasis either in brain, bone, liver, lung, and distant lymph nodes, as well as tumor categorized as M1 stage. The missing rate for each categorical variable is calculated and reported. For those classified data that is unknown, we classify its missing value into the “unknown” category. This processing ensures data integrity and avoids information loss due to missing data. We determined the minimum sample size needed for an external validation cohort by formula of Riley et al. (12).

In the preliminary analysis, variables with P < 0.05 in the univariate and multivariate logistic analysis in the training set were included for the feature selection process. Subsequently, we employed Recursive Feature Elimination (RFE) method based on 6 ML algorithms, involving categorical boosting (CatBoost), random forest (RF), support vector machine (SVM), extreme gradient boosting (XGB), decision tree (DT) and gradient boosting machine (GBM), combined with 5-fold cross-validation, to sift through the clinical features. RFE works by building a model and identifying the most significant features in feature selection phase. This selection process is then iteratively repeated on the subset of remaining features until all features have been evaluated and ranked (13). Then Robust rank aggregation (RRA) algorithm was utilized to integrate the rank of variable importance from six ML algorithms utilized in RFE method to obtain a comprehensive ranking of all variables (14). We set random seed as “2024” in our analysis. In model development phase, we applied 13 ML algorithms, including CatBoost, RF, SVM, XGB, DT, GBM, k-nearest neighbor (KNN), logistic regression (LR), naive bayes classifier (NBC), linear discriminant analysis (LDA), quadratic discriminant analysis (QDA), neural network (NNET) and generalized linear model (GLM) to predict three-year survival via “mlr3” R package (15). This approach allows us to compare the performance of various models and select the best predictive model. To tackle the issue of class imbalance, which could significantly skew performance metrics, we implemented the Synthetic Minority Over-sampling Technique (SMOTE) for training our model (16). We further refined our approach by employing nested resampling, which involved a two-tiered k-fold cross-validation process: one for optimizing model hyperparameters and another nested within it for model selection. Meanwhile, we utilized a 1000-evaluation random search across a 5-fold cross-validation framework, repeated five times for each model. Subsequently, area under the curve (AUC), area under precision-recall curve (PRAUC), accuracy, sensitivity, specificity, precision, cross-entropy, Brier scores and Balanced Accuracy (bacc) and F Beta Score (fbeta) were calculated to select the best ML model. Internal validation was carried out through 5-fold cross-validation. Precision-recall curve (PRC) was employed to evaluate the performance of classification models in handling imbalanced datasets. Calibration curve was utilized to appraise model’s discriminative ability, and decision curve analysis (DCA) was applied to verify the clinic benefit of ML model via “runway” R package (https://github.com/ML4LHS/runway/). We set the selection criteria of our best model: highest AUC, highest PRAUC, and lowest Brier score, while also ensuring a good calibration curve, as well as outperforming balanced accuracy and F Beta Score. To quantify the impact of each variable, we calculated its mean contribution to the AUC as a percentage relative to the full model via “DALEX” R package (17). SHapley Additive exPlanations (SHAP) value were used to explain the best model predictions and to interpret the black-box ML model via “shapviz” R package (https://github.com/ModelOriented/shapviz) (18).

Univariate and multivariate cox analysis were employed to define clinical variables with significant prognosis value in overall survival (OS). We integrated 10 ML algorithms involving random survival forest (RSF), elastic network (Enet), Lasso, Ridge, stepwise Cox, CoxBoost, partial least squares regression for Cox (plsRcox), supervised principal components (SuperPC), GBM and survival support vector machine (survival-SVM) to predict prognosis (in terms of OS) of pancreatic cancer patients. Altogether 101 prognostic ML algorithm combinations were trained in the training cohort, to develop the prognostic ML model according to the leave-one-out cross-validation (LOOCV) framework. Models with <3 clinical variable were removed. Subsequently, the concordance index (C-index) of every ML combination in training, testing and external validation cohorts was obtained (19). The top five ML combinations yielding the highest average C-index across three cohorts were selected for model evaluation via k-fold cross-validation, to mitigate overfitting and ensure the robustness and generalizability of model. Logarithmic loss, recall and decision calibration were utilized to select the best prognostic ML combination via “mlr3proba” R package (20). We incorporated variables from various feature selection patterns to compute risk scores using a linear combination function for each prognostic ML combination. The median risk score from the training cohort was chosen as the threshold to categorize patients in training, testing and external validation cohorts into high or low-risk groups. We utilized the Kaplan-Meier (KM) survival analysis and the log-rank test on these groups, using the “survival” and “survminer” R packages. AUC, time-dependent receiver operating characteristic (ROC) curves, calibration curves and DCA were employed to evaluate the precision, discrimination and clinical benefit of the model.

In the predictive model for three-year survival, a total of 20064 pancreatic cancer patients from SEER database and 103 patients from The First People’s Hospital of Lianyungang were included. We divided patients from SEER database randomly into training and internal validation sets in a 7:3 ratio, respectively. And pancreatic cancer patients from The First People’s Hospital of Lianyungang were assigned as the external validation set. In the trainset from SEER database, 2579 cases (18.3%) were alive at three-year follow-up, while 11548 cases (81.7%) did not. Detailed clinical information regarding the training and validation sets to predict three-year survival can be found in Table 1 . For the outcomes (in terms of OS) of prognostic model, 13157 cases (93.1%) were dead at the time of follow-up, while 970 cases (6.87%) were alive ( Table 1 ). In the training, internal validation and external validation sets, the median follow-up time was 12.0 [5.00;26.0], 12.0 [5.00;27.0] and 16.0 [6.00;30.5] ( Table 1 ). The specific selection process of patients from SEER database is shown in Figure 1 .

We utilized “autoplot” function in “mlr3” R package to visualized the correlation coefficients of the baseline characteristics with three-year survival, which revealed that “AJCC stage” had the most significant correlation with three-year survival ( Figure 2A ). Based on our clinical experiences, we selected 24 variables for the logistic regression analysis ( Table 2 ), while the variable with a correlation coefficient > 0.6 was removed. Subsequently, we performed univariate and multivariate logistic regression analysis in the training cohort to find the effective variables to predict three-year survival, which revealed that “Age” (OR 1.67(1.35-2.06)), “Marital_Status” (OR 1.27(1.14-1.42)), “Household_Income” (OR 0.75(0.68-0.83)), “Histology” (OR 0.5(0.43-0.59)), “Grade” (OR 2.4(1.59-3.72)), “Summary_Stage” (OR 3.12(2.27-4.31)), “Tumor_Size” (OR 2.63(2.08-3.34)), “AJCC_Stage” (OR 1.59(1.09-2.35)), “Surgery_Type” (OR 0.15(0.11-0.2)), “Radiotherapy” (OR 0.79(0.72-0.87)), “Chemotherapy” (OR 0.57(0.52-0.64)), “Lung_Metastasis” (OR 0.2(0.06-0.99)), “M_Stage” (OR 1.26(1.07-1.48)) were significantly powerful to predict three-year survival (P < 0.05, Table 2 ). The correlation analysis between the variables and three-year survival showed that “AJCC_stage” is the most influential factor ( Figure 2B ). Due to high correlation between “AJCC_stage” and “Summary_Stage”, we only choose “AJCC_stage” in the following analysis. Afterwards, we utilized Recursive Feature Elimination (RFE) method based on six ML algorithms (GBM, SVM, RF, DT, XGB and CatBoost), combined with 5-fold cross-validation, to sift through the clinic features ( Figures 2C–H ). Feature selection based on RFE found that the optimal selection was according to GBM algorithm, remaining 12 variables, with the highest AUC (0.819, Figure 2C ). We utilized RRA algorithm to obtain the comprehensive ranking of the clinic variables in six ML algorithms, with the “AJCC_stage” considered most important ( Supplementary Table 1 ). We finally select eight variables with frequencies more than 4, which indicates that these variables are important in most of the ML selection process, into the following procedures of model development ( Supplementary Table 1 ).

To establish a precise model to predict three-year survival, we utilized the eight variables (“AJCC_Stage”, “Chemotherapy”, “Age”, “Grade”, “Lung_Metastasis”, “M_Stage”, “Surgery_Type”, “Tumor_Size”) selected by RFE and RRA. A total of 13 ML models, comprising CatBoost, RF, SVM, XGB, DT, GBM, KNN, LR, NBC, LDA, QDA, NNET and GLM algorithms, were developed by incorporating the above selection of eight variables in the training set. Hyperparameters were fine-tuned by performing 5-cross validation and random searches. Then we evaluated the 13 ML models in the internal validation and external validation cohorts, respectively. Finally, ROC curves analysis found that CatBoost model had the highest AUC in the training (0.932 [0.924, 0.939]), internal validation (0.899 [0.873, 0.934]) and external validation (0.826 [0.735, 0.919]) cohorts ( Figures 3A , 4A , 5A ). CatBoost model has the accuracy of 0.839 [0.831, 0.847], sensitivity of 0.872 [0.858, 0.887], specificity of 0.803 [0.784, 0.825] and precision of 0.832 [0.821, 0.853]. After grid search in hyperparameter tuning, the best hyperparameter metric of CatBoost was depth, 5; learning_rate, 0.01678325; iterations, 548; 12_leaf_reg, 7.409126. The precision-recall curves (PRC) revealed that CatBoost model was powerful in handling imbalanced datasets ( Figures 3B , 4B , 5B ). Calibration plots showed that CatBoost algorithm had the best fitting ability and could accurately predict three-year survival ( Figures 3C , 4C , 5C ). This indicates that the model’s probability estimates are reliable and well-calibrated, as it ensures that the risk estimates provided by the model can be trusted to reflect the true likelihood of patient outcomes. DCA curves suggested that CatBoost algorithm had the best clinical application value and could effectively help predict three-year survival ( Figures 3D , 4D , 5D ). This implies that using the CatBoost model to guide clinical decision-making would result in more effective identification of patients who are likely to benefit from certain interventions, such as more aggressive treatment or intensive monitoring. The accuracy, sensitivity, specificity, precision, cross-entropy, Brier scores, Balanced Accuracy (bacc) and F Beta Score (fbeta) of the 13 ML models were calculated to comprehensively evaluate the model performance, which revealed that CatBoost model was robust in predicting three-year survival ( Figures 3E , 4E , 5E ). The results of tenfold cross-validation indicated that CatBoost exhibited the best performance ( Figure 3F ). Confusion matrix displayed the outstanding predictive ability of CatBoost in the internal validation and external validation cohorts ( Figures 4F , 5F ). Therefore, CatBoost was chosen as the best model for the next step.

We calculated the feature importance rankings of each ML models and illustrated eight of them, including CatBoost, GBM, GLM, NB, KNN, RF, NNET and SVM ( Figure 6A ). The importance scores were determined by leveraging the inherent attributes of various ML algorithms, which revealed that the risk factors most associated with three-year survival were “Surgery Type”, “AJCC Stage” and “M Stage”. Subsequently, we utilized SHAP framework to interpret CatBoost model. We illustrated all of the risk factors evaluated by the mean absolute SHAP value, which revealed that “Surgery Type” was the most impactful variable ( Figure 6B ). Besides, beeswarm plot elucidated the influence of various risk factors on three-year survival ( Figure 6C ). The y-axis denotes the magnitude of the risk factors, while the x-axis represents their impact on the model’s output, specifically three-year survival, as measured by the SHAP value. It was observed that no surgery, higher grade, older age, have lung metastasis, no chemotherapy, higher AJCC stage and M1 stage are associated with an increased likelihood of death in three-year follow-up. To illustrate the model’s interpretability, we highlighted two representative cases. SHAP values were used to understand the impact of each feature on the model’s prediction. In our study, lower SHAP values indicate a higher likelihood of three-year survival, while higher SHAP values suggest a higher probability of death within the three-year follow-up. We chose median value (0.0962) as the cut-off point for predicting the low or high probability of three-year survival. For instance, the first patient with three-year survival had a lower SHAP value and a prediction score of 0.0276, indicating a higher likelihood of three-year survival ( Figure 6D ). In contrast, the second patient without three-year survival showed a higher SHAP value and a prediction score of 0.187, suggesting a higher probability of death in three-year follow-up ( Figure 6E ).

To explore the prognostic values of multiple variables, we performed univariate and multivariate Cox analysis to found that “Sex” (HR 0.934(0.907-0.963)), “Race” (HR 1.071(1.022-1.123)), “Age” (HR 1.19(1.119-1.264)), “Marital_Status” (HR 1.137(1.086-1.191)), “Household_Income” (HR 0.866(0.837-0.896)), “Household_Location” (HR 0.945(0.899-0.993)), “Tumor_Primary_Site” (HR 0.946(0.91-0.983)), “Histology” (HR 0.706(0.66-0.755)), “Grade” (HR 1.334(1.271-1.401)), “Tumor_Size” (HR 1.325(1.238-1.417)), “AJCC_Stage” (HR 0.712(0.643-0.789)), “T_Stage” (HR 1.127(1.024-1.242)), “Surgery_Type” (HR 0.491(0.45-0.535)), “Lymph_Nodes_Surgery” (HR 0.846(0.761-0.94)), “Regional_Lymph_Nodes” (HR 0.709(0.65-0.774)), “Radiotherapy” (HR 0.905(0.873-0.938)), “Chemotherapy” (HR 0.582(0.562-0.602)), “Bone_Metastasis” (HR 1.196(1.006-1.422)), “Liver_Metastasis” (HR 1.338(1.224-1.463)), “Lung_Metastasis” (HR 1.338(1.224-1.463)) and “Metastasis” (HR 0.781(0.709-0.861)) were independent prognosis variables for predicting OS in pancreatic cancer patients (P < 0.05, Table 3 ). Incorporating these clinical variables, 101 prognostic ML algorithm combinations were constructed via LOOCV framework. The C-index of each ML combination was calculated in training, internal validation and external validation datasets ( Figure 7A ). Among top five ML combinations with highest C-index across three cohorts, logarithmic loss, recall and decision calibration were calculated to assess the model performances, discovering the well calibration and precision of “RSF+GBM” model ( Supplementary Figure 1 ). The best ML model combination was “RSF+GBM”, which was established based on RSF algorithm in feature selection ( Figure 7B ), and GBM algorithm in model construction ( Figure 7C ), with the highest average C-index (0.723) across three datasets ( Figure 7A ). Finally, a 20-variable “RSF+GBM” prognostic ML model was accordingly established to predict OS of pancreatic patients, with “Surgery Type” being the most significant variable both in the feature importance visualization of RSF and GBM model ( Figures 7B, C ). ROC curves of 1-, 3- and 5-year OS showed well specificity of “RSF+GBM” model ( Figure 7D ). Time dependent ROC curves indicated that the curve of “RSF+GBM” model was upper than other curves at most of the time points, indicating that “RSF+GBM” model remarkably outperformed conventional clinical variables in capability of discrimination and prediction ( Figure 7E ). Calibration curves ( Figure 7F ) and DCA curves ( Figure 7G ) showed that “RSF+GBM” model is well-behaved in accuracy and clinical benefit. Based on risk scores calculated by GBM algorithm, we utilized the median risk score to divide patients in the training, internal validation and external validation cohorts into low-risk and high-risk groups, respectively. Obliviously, the low-risk group owned a relatively longer OS than the high-risk group in the training, internal validation and external validation cohorts, respectively ( Figure 7H ). The K-M curves validated the capability of risk stratification of “RSF+GBM” model. All these metrics collectively indicated that “RSF+GBM” model demonstrated stability and robustness in model performances. In conclusion, we have successfully developed a “RSF+GBM” model to predict OS in pancreatic cancer patients, which outperforming other models and was well behaved in model performances.