We retrospectively analyzed data from 455 patients who underwent lumbar discectomy at our institution between January 2003 and September 2023. Of 455 patients, we selected 195 patients who underwent follow-up outpatient visits 6 months or later to determine if there was a recurrence. In total, 126 patients were finally included in the analyses after 69 patients were excluded due to the following reasons : 1) 41 patients had no preoperative magnetic resonance imaging (MRI) within 1 month including those who had primary surgery at local hospitals [11], 2) 12 patients herniated within 6 months after lumbar discectomy [5], which suggests that the pain persisted because the surgery was less effective rather than a recurrence, and 3) 16 patients whose preoperative MRI had artifacts that obstructed the measurement of muscle mass. Patients were included in the RLDH group if they had both clinical symptoms and disc herniation confirmed at the surgical level on follow-up MRI after 6 months or later [5]. Thus, 38 patients were included in the RLDH group, and the remaining 88 were included in the non-RLDH group (Fig. 1).

We selected 15 variables that can be easily obtained from the majority of patients prior to lumbar discectomy from electronic medical record data, including demographic and clinical variables, and medical history such as age, sex, height, weight, BMI, level of LDH, smoking, alcohol, DM, hypertension, cerebrovascular accident, and modified frailty scale. We additionally included three LDH operation-related variables such as total operation time, length of stay (LOS), and blood loss.

To investigate the relationship between paraspinal muscles and RLDH, we calculated the volume of muscles using preoperative MRI data. We transformed T2-weighted axial images into the Digital Imaging and Communications in Medicine (DICOM) format. By using the ANOTHEN (NOH THINKING Co., Ltd., Seongnam, Korea) program, we preprocessed individual MRI data by windowing and labeling the location area, and then volume measurements were made. The axial CSA was divided into seven types : PM, MM, QL, ES, vertebra, fat, and spinal cord. In this study, we considered volume measurements for vertebra, fat, and spinal cord in addition to paraspinal muscles.

To further evaluate the generalizability of the final model, we conducted external validation using an independent dataset obtained from another institution. This dataset was not involved in model training or internal cross-validation (CV). For external validation, we used data from 62 patients who underwent lumbar discectomy at a different tertiary care hospital between 2015 and 2022. Inclusion and exclusion criteria were consistent with those of the internal cohort. Of these, 18 patients developed recurrent LDH and 44 did not.

Descriptive statistics for baseline characteristics of demographic variables were calculated by considering their data types. We conducted Mann-Whitney U-test for numerical variables and implemented the chi-squared test and Fisher’s exact test for categorical variables. A significance level of p<0.05 was set as the threshold for statistical significance. In the case of MRI measurements, we conducted pre-test odds ratio analysis to determine the statistical association of paraspinal muscle volumes with LDH. The volume of paraspinal muscles was converted from a continuous variable to a categorical variable based on a cutoff value to proceed with multivariable logistic regression (LR). The optimal cutoff was calculated based on receiver operating characteristic (ROC) curve (Supplementary Table 1). For LR analysis, muscle volumes were dichotomized into categorical variables using optimal cut-off values derived from ROC analysis. In contrast, for ML model training, all muscle volume variables were entered as continuous values to preserve the full granularity of the data.

We considered feature selection for dimension reduction to boost the accuracy of the RLDH prediction model. For feature selection, we used the combination of two methods, Boruta selection [8] and forward selection methods to effectively select significant features. The Boruta selection algorithm is an ensemble-based feature selection algorithm that identifies all relevant variables in a dataset by iteratively comparing the importance of real features to that of randomly shuffled copies. To reduce bias from feature importance calculation, we incorporated the Boruta selection and forward selection methods using a nested 5-fold CV scheme. The initial features as accepted by the Boruta algorithm were selected, and then subsequently, additional features were incorporated into the forward selection process, with each feature’s influence on model development being recorded. The Boruta and forward selection methods identified the following features as significant for ML model training: weight, QL volume, PM volume, vertebral body volume, fat volume, LOS, and operation time. These features were consistently retained across CV folds and were used as input variables in the final ML models (Supplementary Fig. 1).

Using selected features from training sets, five supervised ML algorithms were used for RLDH risk prediction : LR, extreme gradient boosting (XGBoost), light gradient boosting machine (LightGBM), categorical boosting (CatBoost), and multilayer perceptron (MLP). The task of RLDH risk prediction is a binary classification task, so we used the area under ROC curve (AUROC) as the primary model evaluation metric for this study. Grid search was used to optimize the hyperparameters of each model until an optimal model was identified to maximize the area under the curve (AUC) score. The average AUC of the five test results was reported to indicate the final performance of the model. The overall performance for the AUC score of each test data was interpreted based on a 95% confidence interval (CI). Additional evaluation metrics such as accuracy, F1-score, sensitivity, and specificity were employed to assess the robust performance of each model. All statistical analyses were performed with R version 4.3.0 (R Foundation for Statistical Computing, Vienna, Austria), and ML analyses were performed with Python version 3.11.8 (Python Software Foundation, Wilmington, DE, USA).

A total of 126 patients underwent lumbar discectomy and met the inclusion criteria of our study (64 females [50.8%] and 62 males [49.2%]). The demographic data are presented in Table 1. The majority of patients who underwent successful lumbar discectomy did not undergo outpatient follow-up, which led to difficulties in obtaining a sufficient sample size. The patients’ average age was 47.9±15 years. Between the RLDH and non-RLDH groups, only LOS was significantly different.

The results of the LR analyses are shown in Table 2. In the univariate LR, the volumes of QL and fat were significant. Only QL was significant in multivariate LR with an odds ratio of 7.894, which implies that patients with QL greater than 55 mL were 7.894 times more likely to have RLDH than patients with QL less than 55 mL.

The comparative performance of five ML models for RLDH risk prediction is summarized in Table 3. XGBoost demonstrated the best overall performance, with the second highest AUC (0.811; 95% CI, 0.791-0.832), accuracy (0.794; 95% CI, 0.782-0.807), balanced sensitivity (0.658; 95% CI, 0.630-0.6987), and specificity (0.894; 95% CI, 0.883-0.904). CatBoost also performed well, with the highest AUC of 0.813, but its accuracy and F1-score were lower than those of XGBoost. Light-GBM had the lowest AUC, accuracy, and F1-score compared to the other boosting-based models, indicating poorer performance overall. Additionally, the deep learning model MLP and the linear model LR, both showed the poorest overall performance. These results suggest that XGBoost provides the most reliable and balanced predictive performance for RLDH risk among the models tested, particularly in contexts where both sensitivity and specificity are crucial. Fig. 2 illustrates the ROC plots prepared using ML-based predictive models.

For RLDH risk prediction, the XGBoost model was selected as the final predictive model due to its superior performance across multiple evaluation metrics. To interpret the predictions extracted by XGBoost and to deeply understand the contribution of individual features to the risk of RLDH, we employed Shapley Additive exPlanations (SHAP) values [12]. SHAP values offer a unified measure of feature importance by quantifying the contribution of each feature to the model’s output; it allows for a transparent post-hoc analysis, providing insights into how each feature influences the prediction of RLDH risk. Fig. 3 compares the feature importance rankings derived from XGBoost’s internal methods with those obtained using SHAP values. Despite both methods showing significant differences in terms of importance, the most common influential variables of both methods were LOS and L5/S1. In terms of SHAP values (Fig. 3B), the Weight variable was the most important risk factor for RLDH prediction, and this variable had the greatest impact on the outcome of XGBoost (Fig. 3A), followed by vertebra and fat, which were also significant for the predictions. Additionally, clinical and operational factors such as LOS, specific spinal levels, and operation time contributed meaningfully to the model’s decisions. Overall, the results emphasize that predictions of XGBoost are driven by a combination of anatomical features, clinical history, and procedural variables, underscoring the complexity and multifactorial nature of the predictive process (Supplementary Fig. 2).

A sample-wise force plot was then prepared to visualize and assess the positive-negative correlation between the features and each prediction value (Fig. 4). In terms of sample-wise analysis, Fig. 4A shows the model’s contribution of each variable toward predicting a value of 1. In this case, variables such as Weight and QL significantly increase the predicted value, thereby playing a crucial role in steering the model toward a prediction of 1. As shown in Fig. 4B, the contributions were analyzed when the model predicted a value of 0. Here, Weight and LOS substantially reduced the predicted value, serving as key determinants in guiding the model toward a prediction of 0. Overall, the Weight variable considerably influenced the model’s predictions, demonstrating its importance across different prediction outcomes.

Table 4 presents the external validation results performed on an independent dataset obtained from another institution, which was not used for model training. The final model selected was used for RLDH prediction. The performance indicators show the generalizability of the model, maintaining consistent accuracy, precision, recall, and AUROC compared to the internal validation results.

First, the sample size of the study was small, which led to insufficient statistical power and caused data imbalance and poor generalization performance in our ML models. This is a problem caused by the loss of follow-up of lumbar discectomy patients. Out of 455 patients, 260 were lost to follow-up. This is likely because the majority of patients with successful lumbar discectomy and resolved symptoms do not return to the hospital. Therefore, if we increase the sample size by combining data from other institutions, we can expect better performance of the ML models. Second, the retrospective design of the study necessitates the addition and updating of data. Moreover, we were unable to analyze variables for which we had no data. Analysis of the values of more risk factors and their use in training the ML models would have improved the model performance.