This retrospective cohort study was approved by the Ethics Committee of the First Affiliated Hospital of Chongqing Medical University (Ethics Approval Number: 2024-087-01). Informed consent was waived due to the retrospective design. We reviewed 336 patients who underwent PD between October 2018 and June 2023. Inclusion criteria were: (1) complete clinical and pathological data, (2) preoperative contrast-enhanced CT within 1 month before surgery. Exclusion criteria included: (1) non-curative resection, (2) prior neoadjuvant therapy, (3) poor CT image quality. After screening (Figure 1), 241 patients were included and stratified into CR-POPF (N = 55, 22.8%) and non-CR-POPF (N = 186, 77.2%) groups based on ISGPS 2016 criteria (Bassi et al., 2017).

Demographic and clinical comparisons between groups are summarized in Table 1. Age, gender, diabetes, hypertension, cardiovascular/pulmonary diseases, smoking history, and prior abdominal surgery showed no significant differences (P > 0.05). However, CR-POPF patients exhibited higher BMI (23.2 vs. 22.2 kg/m², P = 0.035), increased alcohol consumption (43.6% vs. 28.0%, P = 0.042), and smaller MPD diameter (2.77 mm vs. 4.25 mm, P < 0.001). Preoperative laboratory tests revealed elevated platelet-to-albumin ratio (PAR: 6.55 vs. 5.53, P = 0.012) and bilirubin levels (129 vs. 78.2 μmol/L, P = 0.004) in the CR-POPF group. Pancreatic head lesions were more frequent in CR-POPF patients (67.3% vs. 32.7%, P = 0.006). No differences were observed in preoperative biliary drainage, ASA classification, or surgical approach (P > 0.05).

The cohort was randomly split into training (n = 193, 80%) and test sets (n = 48, 20%) using an 8:2 ratio. Reporting followed TRIPOD guidelines (Collins et al., 2015).

Contrast-enhanced abdominal CT scans were performed using Siemens SOMATOM Force, GE Discovery CT750 HD, or GE LightSpeed VCT. Scanning parameters: 120 kV, 200 mA, 5 mm slice thickness. All images were reconstructed using a standard reconstruction kernel with the following parameters: pitch of 1, rotation time of 0.5 s, field of view of 350 mm × 350 mm, matrix size of 512 × 512, slice thickness of 5 mm, interval of 5 mm, and reconstruction slice thickness of 1 mm. Patients were required to fast and avoid drinking for at least 3 h prior to the examination. A non-ionic iodinated contrast agent (300–400 mgI/ml) was administered intravenously at a dose of 1–1.5 mL/kg with an injection rate of 3 mL/s. Arterial phase scanning was delayed by 15–18 s. Portal venous and delayed phase scans were performed with delays of 33–36 s and 180 s, respectively. Enhanced CT images were exported from the Picture Archiving and Communication System (PACS) in DICOM format for further analysis.

Image preprocessing included artifact removal, grayscale normalization (0–255), and enhancement via contrast adjustment, sharpening, and noise reduction.

Two radiologists (>5 years of experience) manually delineated pancreatic parenchyma (body and tail) as regions of interest (ROIs) on portal venous phase images using ITK-SNAP (v3.6.0). The portal vein served as the anatomical landmark to differentiate the pancreatic head from the body. Segmentation masks were saved in Nifti format. A senior radiologist (>10 years of experience) validated 50 randomly selected samples. Intraclass and interclass correlation coefficients (ICCs) were calculated, with ICC >0.8 indicating satisfactory reproducibility.

Radiomics feature extraction was performed using the PyRadiomics library (v3.0.1) in Python, based on original CT images and their preprocessed variants, including those filtered with Laplacian of Gaussian (LoG) and wavelet transforms. The extracted features encompassed first-order statistics (e.g., mean, variance, skewness), shape features (e.g., volume, sphericity, maximum diameter), and texture features derived from matrices such as the gray-level co-occurrence matrix (GLCM), gray-level run-length matrix (GLRLM), gray-level size zone matrix (GLSZM), and neighborhood gray-tone difference matrix (NGTDM).

To optimize feature selection and reduce dimensionality, the Lasso regression combined with cross-validation (Lasso-CV) was applied. Regularization parameters were optimized using grid search, and features with non-zero coefficients were retained. Additionally, the Random Forest (RF) algorithm was employed to rank the importance of clinical variables and identify the most relevant features for model development.

In the training cohort, five machine learning predictors—XGBoost, Random Forest (RF), Extra Trees (ET), Gradient Boosting (GB), and AdaBoost—were employed to develop 15 CR-POPF prediction models (the model parameters are presented in Table 2). These models were trained on three datasets: radiomics-only, clinical-only, and a combined radiomics-clinical dataset. Feature selection and model training were performed exclusively on the training set. The test set remained entirely independent and was only used for final model evaluation to prevent data leakage. Hyperparameters, including learning rate, maximum tree depth, subsampling rate, and regularization terms, were optimized via 5-fold cross-validated grid search to balance model complexity and generalizability. To mitigate overfitting, early stopping mechanisms and maximum iteration limits (10,000 iterations) were enforced during training. Model performance was rigorously evaluated using accuracy, AUC, precision, recall, and F1 score.

Pairwise comparisons of AUC values between models were conducted using DeLong’s test, with results visualized as a heatmap (Supplementary Material S1) to highlight statistically significant differences (P < 0.05). Calibration curves quantified the agreement between predicted probabilities and observed outcomes, while decision curve analysis (DCA) assessed clinical utility by quantifying net benefits across threshold probabilities (Supplementary Material S2). Model interpretability was enhanced via SHAP analysis, elucidating feature contributions globally and locally. The workflow is summarized in Figure 2.

To comprehensively evaluate the performance of the predictive models, five standard metrics were employed: accuracy, precision, recall, F1-score, and AUC. The definitions and corresponding formulas are as follows:

• Accuracy: The proportion of correctly classified instances relative to the total instances (Sokolova and Lapalme, 2009).

The statistical evaluations were executed employing Python software (version 3.7; https://www.python.org/). Quantitative data, conforming to a normal distribution, are articulated as the mean ± standard deviation (SD), while quantitative data that do not follow a normal distribution are represented as the median, along with the interquartile range. Categorical data are denoted as numbers and percentages (N, %). To assess the efficacy of the constructed models, several widely utilized metrics were chosen, encompassing accuracy, precision, recall, F1 score, and the area under the Receiver Operating Characteristic (ROC) curve (AUC). Pairwise comparisons of AUC values between models were conducted using DeLong’s test to assess statistical significance. The threshold for statistical significance was established at P < 0.05.

In the course of this study, the Pyradiomics library was utilized to derive 1719 radiomics features from CT images. To guarantee the performance and interpretability of the model, Lasso regression was implemented for the selection of these high-dimensional features. The alterations in the model’s performance with varying parameter α iterations were depicted in Figures 3A, B, thereby determining the optimal parameter values and the corresponding number of features. The Lasso regression model was then used to pinpoint features with non-zero coefficients, which were subsequently ranked based on the absolute values of these coefficients, as illustrated in Figure 3C, including texture features (e.g., wavelet HHL_glszm_GrayLevelNonUniformity, original_glcm_ClusterProminence) and shape features (e.g., original_shape_Maximum2DDiameterRow, original_shape_Sphericity).

For clinical features, RF analysis ranked eight predictors by importance (Figure 4): MPD diameter (2.77 mm vs. 4.25 mm in CR-POPF vs non-CR-POPF), lesion location (67.3% vs. 32.7% pancreatic head involvement), bilirubin (129 vs. 78.2 μmol/L), PAR (6.55 vs. 5.53), ALT (95.0 vs. 137 U/L), BMI (mean: 23.2 vs. 22.2 kg/m²), systemic immune-inflammation index (SII) (820 vs. 688), and AST (82.0 vs. 87.0 U/L).Combining radiomics and clinical features, we established a multimodal feature set encompassing 28 variables for subsequent training and validation of ML models.

A total of fifteen CR-POPF prediction models were developed using five distinct ML predictors, incorporating selected radiomics features, selected clinical features, and a combination of both. Among the radiomics-based models, the AdaBoost model demonstrated the highest predictive performance, achieving the highest AUC of 0.87 (Figure 5A), along with the best recall (0.76) and precision (0.71). In contrast, the RF and Extra Trees models exhibited the highest accuracy (0.83); however, the RF model showed lower robustness, with an AUC of 0.74. Regarding clinical models, the Extra Trees model achieved the highest AUC (0.85, Figure 5B) while maintaining a balanced performance in terms of precision (0.75) and recall (0.75).

The CR-POPF prediction model based on the XGBoost predictor with the combination of the selected radiomics and clinical features demonstrated superior performance among all 15 CR-POPF prediction models, achieving an accuracy of 0.85 and an AUC of 0.93 (Figure 5C). Compared to unimodal models, the combined model exhibited statistically significant improvements. Specifically, when compared to the radiomics-only model using the Extra Trees predictor, the AUC difference was 0.17 (p = 0.041), and when compared to the clinical-only model using the Random Forest predictor, the AUC difference was 0.11 (p = 0.041). These results, validated by DeLong’s test (Supplementary Material S1, where red cells denote p < 0.05), highlight the synergistic value of multimodal integration.

Detailed performance metrics are summarized in Table 3, while calibration and decision curve analysis (DCA) curves (Supplementary Material S2) further validated the clinical utility of the combined model across threshold probabilities. These results highlight that the integration of CT radiomics and clinical data significantly enhances preoperative CR-POPF risk stratification.

In the implementation of the XGBoost ensemble model, the SHAP method is utilized to elucidate the final model output through the computation of each variable’s contribution to the prediction. This interpretive strategy yields two categories of explanations: global explanations at the feature level and local explanations at the individual level. Global explanations elucidate the comprehensive behavior of the model and the significance of its features. This is illustrated in the SHAP bar chart and the SHAP summary plots (Figures 6A, C), where the influence of features on the model is assessed via mean SHAP values and presented in a descending sequence, thereby highlighting the top 20 variables that contribute most significantly to the model. The three variables with the highest contribution are wavelet-HHL_glszm_GrayLevelNonUniformity, original_shape_Maximum2DDiameterRow, and lesion location. The SHAP heatmap (Figure 6B) visually represents the direction and magnitude of the effect of each feature across all instances within the model. Additionally, SHAP dependence plots (Figure 7) facilitate comprehension of the manner in which a singular feature influences the output of the XGBoost predictive model. The y-axis denotes the SHAP value of the feature, in contrast to the x-axis, which signifies the value of the feature. The plot provides a visual representation of the fluctuating importance of the feature in relation to its value. A SHAP value exceeding zero corresponds to positive class predictions within the model, signifying an elevated risk of CR-POPF. Local explanations scrutinize the methods by which specific predictions for individual cases are formulated through the amalgamation of personalized input data. Figure 8 delineates instances of four standard positive and negative CR-POPF forecasts. The SHAP Waterfall plot elucidates the contributions of each attribute to the prediction outcome for a singular case. The baseline value symbolizes the model’s fundamental prediction probability, while each feature’s contribution value (also known as the SHAP value) signifies the direction and magnitude of that particular feature’s influence on the prediction. Positive values imply that the feature escalates the likelihood of predicting positive CR-POPF. The final prediction probability, denoted as f(x), is the cumulative sum of the baseline value and all feature contributions.

The integration of radiomics and clinical features through ML offers a transformative approach for preoperative prediction of CR-POPF. In this study, we extracted 1,719 radiomics features from preoperative portal venous phase CT images of 241 PD patients and combined them with clinical variables to develop a multimodal predictive model. The dual application of Lasso regression and RF algorithm for feature selection proved instrumental in balancing dimensionality reduction with biological relevance. Lasso’s regularization properties efficiently distilled 1,719 radiomics features to 20 non-redundant predictors, mitigating overfitting while preserving texture and shape parameters critical for quantifying pancreatic heterogeneity—a strategy validated in pancreatic cancer studies by Kim et al. (2019). Meanwhile, RF’s inherent ability to rank nonlinear interactions among clinical variables identified MPD, lesion location, PAR, and other important features as key contributors, reflecting anatomical risk and systemic inflammation, respectively (Huang et al., 2022; Tan et al., 2022). This hybrid approach harmonizes the strengths of both methods: Lasso’s sparsity induction for dimensionality reduction and RF’s robustness in handling multicollinearity, aligning with methodological frameworks advocating combined techniques for high-dimensional biomedical data (Azodi et al., 2020; Kumarasamy et al., 2021).

The experimental results underscore the differential performance of ML predictors when utilizing single-modal versus multimodal features. Models trained solely on selected radiomics features achieved moderate predictive accuracy (AUC: 0.74–0.87), with texture parameters such as GLSZM and Gray Level Dependence Matrix (GLDM) emerging as pivotal predictors, consistent with studies emphasizing their utility in quantifying tissue heterogeneity and fibrosis—key determinants of pancreatic anastomotic integrity (Lubner et al., 2017; Chitalia and Kontos, 2019; Kim et al., 2019). For instance, Abunahel et al. linked GLSZM features to pancreatic stiffness, a surrogate for soft pancreatic texture widely associated with CR-POPF (Abunahel et al., 2022). Similarly, Capretti et al. reported comparable AUCs (0.75–0.81) using CT texture analysis, underscoring the reproducibility of radiomics in pancreatic risk stratification (Capretti et al., 2022). However, the inherent limitations of unimodal radiomics models—such as their inability to incorporate systemic physiological variables—highlight the necessity of integrating clinical data to enhance generalizability. Our clinical-only model, incorporating variables such as MPD diameter, lesion location, and PAR, achieved an AUC of 0.82–0.85. While this performance aligns with the predictive capacity of established risk scores like the FRS and updated alternative FRS (ua-FRS) (Mungroop et al., 2019; Mungroop et al., 2021), it demonstrates a moderate improvement over their external validation results (AUC: 0.74–0.82), highlighting the potential advantages of integrating modern ML frameworks with preoperative clinical indices. Notably, the significantly smaller MPD diameter in CR-POPF patients (2.77 vs 4.25 mm, P < 0.001) reflects multifactorial pathophysiology involving impaired drainage, reduced fibrosis-mediated anastomotic stability, and elevated duct-to-mucosa tension, synergistically increasing fistula risk (Casciani et al., 2021; Lee et al., 2023). Despite these strengths, clinical models struggle to capture subvisual parenchymal changes, such as microlobular fat infiltration or ductal microcalcifications, which radiomics excels in detecting (Chitalia and Kontos, 2019; Abunahel et al., 2022). This limitation highlights the necessity of integrating multimodal data to address the multifactorial nature of CR-POPF pathogenesis.

The multimodal XGBoost model (AUC: 0.93) outperformed all unimodal approaches, underscoring the synergistic value of combining radiomics and clinical data. This aligns with emerging paradigms in precision oncology, where combined models consistently outperform unimodal approaches by encapsulating both macroscopic pathophysiology and microenvironmental heterogeneity (Shen et al., 2022; Verma et al., 2024). However, ML techniques are frequently characterized as “black boxes,” with limited studies dedicated to elucidating the sources of their predictions. This underscores an additional advantage of our study: following the training and evaluation of the model, we employed SHAP methods to interpret the “black box” nature of the ML model. By presenting the SHAP values, we elucidated the relationship between critical covariates and the estimated risk of CR-POPF: wavelet-HHL_glszm_GrayLevelNonUniformity (reflecting parenchymal disorganization) and MPD diameter jointly drove predictions, mirroring the interplay between ductal anatomy and tissue integrity. Such findings resonate with Lambin et al.’s assertion that radiomics bridges qualitative imaging and quantitative biology, thereby advancing clinical decision-making (Lambin et al., 2017).

Furthermore, case analysis elucidates the contributions of critical features within individual cases and computes the final Shapley values to derive the ultimate prediction probabilities, thereby facilitating personalized predictions. For patients at high risk of CR-POPF, preoperative preventive strategies, including nutritional support, optimization of diabetes and exocrine insufficiency, and respiratory training, may confer substantial benefits (Ausania et al., 2019; Bundred et al., 2020). Additionally, prophylactic medications, such as somatostatin analogs or hydrocortisone, have demonstrated efficacy in reducing complications associated with pancreatic surgery (Allen et al., 2014; Laaninen et al., 2016; Tarvainen et al., 2020). Risk assessment identifies patients best suited for interventions, cutting unnecessary medication costs. Evaluating the risk of CR-POPF also facilitates the management of drainage by enabling the early removal of drains in low-risk patients, consequently diminishing the risks of infection and erosion (Conlon et al., 2001; McMillan et al., 2017). Such personalized preventive measures are essential for mitigating the adverse effects associated with CR-POPF.

Undeniably, our study has several limitations. First, the retrospective study design may introduce selection bias. Second, due to the model being derived from a single center, the sample size is relatively small, and external applicability needs further testing. Third, although our ML model can be used to assess the risk of CR-POPF in precision medicine, too many features limit its clinical application. Future studies should validate this framework prospectively and explore streamlined feature sets to facilitate real-world deployment.