In this study, 409 pelvic patients (220 cervical cancer, 73 rectal cancer, 91 uterine cancer, and 25 others) who received IMRT or VMAT at the Radiotherapy Centre of Hunan Cancer Hospital from November 2020 to December 2022 were retrospectively collected. The QA plan was calculated with the Pinnacle³ treatment planning system (Version 9.2, Philips) for 196 cases of IMRT. The dose grid was 3 mm, and Delta4 device (ScandiDos, Sweden) was used to perform dose verification on the Varian linear accelerator. The QA plan was calculated with Eclipse (Version 13.6, Varian) and Monaco (Version 5.11.03, Elekta) treatment planning systems for 213 cases of VMAT. The dose grid was 3 mm, and the ArcCHECK device (Sun Nuclear, United States) was used to perform dose verification on the Varian and Elekta linear accelerators. The linear accelerators and the measuring devices were regularly calibrated during the measurement period to ensure that the equipment is in a good performance state.

As recommended by the American Association of Physicists in Medicine (AAPM) TG 218 report [19], the mean value of GPR was 96.2% ± 3.2% (the range was 78.8%–100%) based on the criterion of absolute dose, 3%/2 mm, global normalization, and a 10% dose threshold. For the classification model, the setting of the tolerance value of the treatment plan “pass” and “fail” affects the performance of the model [20]. In order to build a classification prediction model with better performance for the data of this institution, 99.5% of the average measured GPR, i.e. 95.7%, was used as the threshold of the GPR classification. When the GPR was greater than this threshold, the measured GPR was expressed as “pass” and recorded as “1”, otherwise, “fail” and recorded as “0”.

Radiomics features refer to the semi-quantitative and/or quantitative features extracted from radiography (medical images), which, combined with artificial intelligence technology, play important roles in radiotherapy [21]. The 10% isodose line of the maximum dose was included as the area where the radiomics features were extracted in this study. Batch extraction of the features was performed using the radiomics library in Python 3.7. The image types included the original image (Original), the wavelet transform image (Wavelet), and the Gaussian filter image (LoG). There were seven different types of features: shape features 2D/3D, first-order features, gray-level co-occurrence matrix (GLCM), gray-level size zone matrix (GLSZM), gray-level run length matrix (GLRLM), neighboring gray tone difference matrix (NGTDM), and gray-level dependence matrix (GLDM); altogether, 1,130 features (Original:107, Wavelet:744, and LoG:279) were extracted in total.

The whole dataset is randomly divided, with 90% of the data (368 cases) being used as the training set and 10% being used for the test set. Due to the imbalance of the data, the stratified sampling method was used, making the proportion of all kinds of data in the training set and the test sets consistent with the original data. Feature selection technology is a key link in building a machine learning prediction model based on radiomics, which can avoid high-dimensional data disaster problems, reduce the training time, increase the interpretability of the model, and enhance the prediction performance of the model [22]. In addition to interpreting the output of machine learning models, the SHAP value can also be used as a feature selection method for processing high-dimensional data [23]. In this work, the SHAP value, combined with the XGBoost algorithm, was used in feature selection.

The training set was input into the XGBoost model, and the SHAP value of each feature in the sample was subsequently calculated to measure feature importance, which was obtained by averaging the contributions in all possible per-mutations of the special collection [24]. The SHAP value of feature i was defined, as shown in Eq. 1. φ i = ∑ S ⊆ N \ i S ! N − S − 1 ! N ! ν S ∪ i − ν S , (1) where N denotes the feature sets of the original data and S represents any feature subset in N. S ⊆ N \ i represents a subset of all elements in the sequence before feature i, ν S represents the output of a machine learning model for a feature subset S, and ν S ∪ i − ν S denotes the cumulative contribution of feature i. The sequence number of features started from 0 by default, and the top 45 features (See Supplementary Data Sheet S1) were ultimately selected as the optimal feature subset to be input into the four classification prediction models.

Normalization was performed on the training set, and this transformation was subsequently applied to the testing set to prevent information leakage from the testing data. Four tree-based machine learning algorithms, namely, RF, AdaBoost, XGBoost, and LightGBM, were selected to fit the training data. Grid searching [25] and five-fold cross-validation were used to obtain the model with the highest performance parameters applied to the test data.

RF is a special bagging method, which finds the optimal solution among randomly selected features of a decision tree to split each node and integrates these predictions of the decision tree to avoid overfitting of the model [26]. AdaBoost is a practical boosting algorithm, which creates a highly accurate classifier by adjusting a relatively weak and inaccurate combination of weights for the same training set [27]. XGBoost is an improved algorithm based on the gradient boosted decision tree (GBDT), where the entire dataset is used to generate each decision tree and the residuals between the predicted and true results of the previous decision tree model are taken into account in the generation of the latter decision tree [28]. LightGBM mainly proposed the gradient-based one-side sampling algorithm, mutually exclusive feature binding algorithm, parallel features, and data to solve the multi-feature problem of large data encountered in practical applications [29].

The performance of the binary classification model was evaluated using precision, sensitivity, specificity, F1 score, and the area under the curve (AUC). The curve is the receiver operating characteristic curve (ROC). Precision indicates the ratio of correctly predicted positive instances to the total number of instances predicted as positive (Eq. 2). Sensitivity represents the ratio of correctly predicted positive instances to the actual number of positive instances (Eq. 3). Specificity represents the ratio of correctly predicted negative instances to the actual number of negative instances (Eq. 4). The F1 score is a measure of a model’s accuracy that takes into account both precision and recall (Eq. 5). The ROC is a curve with the false positive rate at different thresholds as the horizontal coordinate and the true rate as the vertical coordinate, and the AUC value represents the area of the region below the ROC curve. TP and FP represent the number of positive and negative samples which are predicted as positive. TN and FN represent the number of positive and negative samples which are predicted as negative. The modeling and analysis procedures were performed using Python 3.7. p r e c i s i o n = T P T P + F P , (2) s e n s i t i v i t y = T P T P + F N , (3) s p e c i f i c i t y = T N F P + T N , (4) F 1 − s c o r e = 2 * p r e c i s i o n * s e n s i t i v i t y p r e c i s i o n + s e n s i t i v i t y . (5)

Figure 1 shows the confusion matrix of the four classification prediction models. The classification performance of each model on the testing set can be calculated based on the confusion matrix, as shown in Table 1. The results show that the RF model achieves a sensitivity of 0.96, and the specificity of the XGBoost and LightGBM models was 0.62. The precision and F1 score of the XGBoost model were 0.84 and 0.88, respectively. Figure 2 shows the ROC curves of the four classification prediction models, where the AUC values of RF, AdaBoost, XGBoost, and LightGBM were 0.81, 0.77, 0.85, and 0.83, respectively.

Figure 3 shows the SHAP summary plots for the four different models on the test set. The importance of the input features was ranked by SHAP values, where the most important features in the RF, AdaBoost, XGBoost, and LightGBM models were features 41, 19, 3, and 2, respectively. The higher the ranking of features, the greater the influence on the model output, and the overall influence of each feature on the model output can be observed. Different colors represent feature values (high values are in red, and low values are in blue), and wide areas indicate large sample clusters. As shown in Figure 3C, most of the blue points of feature 3 were distributed in regions with positive SHAP values. The lower value of feature 3 will have a positive driving effect on the model output and improve the probability of plan passing. Most of the red points in feature 1 were distributed in regions with positive SHAP values, indicating that the higher value of feature 1 increases the probability of plan passing. Table 2 shows the names of the top 10 significant features of the XGBoost model.

The SHAP force plot for two samples under the XGBoost model is shown in Figure 4. The length of the arrow in the figure indicates the magnitude of the feature’s impact, where red represents a positive effect and blue represents a negative effect on the final output of the sample. For example, Figure 4A shows the force plot of sample 0 in the test set under the XGBoost model, demonstrating the influence of each input feature on this sample’s predicted output. The model’s base_value was 1.074, and the predicted output_value f(x) for this sample was 0.915. The difference between the total length of the red arrows and that of the blue arrows equals the distance between the base_value and the output_value. Figure 4B shows the SHAP force plot for sample 1 that has the same base_value, and its predicted output_value f(x) was −0.853.