Enhancing the Diagnosis of Discogenic Low Back Pain: A SHAP-Explained Random Forest Model Combined with Lumbar Spine MRI T2 Radiomics

Summary

This study enhances discogenic low back pain (DLBP) diagnosis using lumbar spine MRI T2 data, radiomics, and a Random Forest model with SHAP analysis, achieving high accuracy and interpretability for improved clinical decision-making.

Abstract

Low back pain (LBP) is a leading cause of disability and reduced quality of life globally, with discogenic low back pain (DLBP) accounting for 39% of cases. Accurate diagnosis of LBP etiology is challenging due to the lack of reliable methods. This study aims to improve DLBP diagnostic efficiency using lumbar spine MRI T2 data combined with radiomics and machine learning. This retrospective study analyzed MRI data from 81 DLBP patients and 162 healthy controls. Radiomics features, clinical data, and high-intensity zone (HIZ) imaging features were extracted. The data were divided into four groups (d0, d1, d2, D), and 20 predictive models were built using Random Forest (RF), Decision Tree (TREE), Support Vector Machine (SVM), K-Nearest Neighbors (KNN), and Logistic Regression (LOG). Model performance was evaluated using Receiver Operating Characteristic (ROC) area under the curve (AUC), precision recall (PR) AUC, accuracy, sensitivity, specificity, positive predictive value (PPV), negative predictive value (NPV), and F1 score. SHapley Additive exPlanations (SHAP) were applied to interpret the most significant features. The Random Forest in group D showed the best performance, with ROC AUCs of 0.9861 (train) and 0.9580 (test), PR AUCs of 0.9813 and 0.9179, and F1 scores of 0.9254 and 0.8148, respectively. SHAP analysis identified first-order kurtosis as the top feature contributing to DLBP diagnosis. The Random Forest model with SHAP analysis significantly improved DLBP diagnosis, offering high performance and interpretability to enhance clinical decision-making.

Introduction

With the aggravation of population aging and the widespread occurrence of improper sitting postures in daily work and life, the number of patients with low back pain (LBP) is increasing year by year¹. The average prevalence of LBP in the general adult population is approximately 12%, with higher prevalence observed in individuals aged 40 or older and in women, with a lifetime prevalence of about 40%². Among all LBP cases, discogenic low back pain (DLBP) accounts for approximately 39%, with its etiology being complex and diverse, often associated with disc degeneration, particularly annular tears that reach the outer annulus of the disc^3,4,5. Increased disc pressure stimulates nerve endings sensitized by inflammatory mediators, leading to discogenic pain^5,6.

Currently, discography is considered the gold standard for DLBP diagnosis, confirming etiology by provoking discogenic pain⁷. However, its validity and safety are questioned due to potential false positives and the risk of accelerating disc degeneration, limiting its clinical application^5,8. The high-intensity zone (HIZ) refers to high-intensity signals in the posterior annulus on T2-weighted MRI images, pathologically characterized by annular tears or ruptures^5,8,9. Although the relationship between HIZ and DLBP has been widely studied, HIZ has low specificity and is also common in asymptomatic individuals, often leading to misdiagnosis^10,11,12,13. In T2-weighted images of lumbar degenerative disc disease (LDDD) and chronic LBP patients, nearly all degenerated discs show reduced signal intensity¹⁴. According to Pfirrmann et al.'s classification, disc degeneration is graded from I (normal disc) to V (severe degeneration)¹⁴. T2-weighted imaging is widely used in multimodal studies for the diagnosis and evaluation of LDDD and LBP, providing critical insights into disc degeneration and its clinical relevance^15,16. However, low signal intensity does not accurately reflect morphological changes in the disc and has little correlation with the degree of pain caused by DLBP^17,18,19.

Currently, clinical decision-making for DLBP primarily relies on the subjective assessment of clinicians, who focus on factors such as HIZ, disc height, black discs, or signal changes in the disc^20,21. Additional considerations include patient age, gender, height, weight, and body mass index (BMI). The duration and intensity of LBP are also within the scope of consideration. Such a diagnosis method, dependent on subjective assessment, may lead to erroneous treatment plans. Therefore, developing an effective diagnostic tool to improve the diagnostic efficiency of DLBP is urgently needed⁴.

MRI, as a non-invasive imaging technique, plays an important role in the diagnosis of LBP. In particular, axial T2-weighted imaging (T2WI), with its high tissue contrast, provides more accurate physiological information and pathological status of the disc compared to traditional computed tomography (CT)^22,23. Radiomics feature analysis, by converting image data from MRI, CT, etc., into quantitative features for advanced mathematical analysis, provides a powerful tool for the evaluation and interpretation of medical images²⁴. This method can reduce biases from subjective analysis, thereby enhancing the accuracy and consistency of diagnostic results²⁵.

In recent years, artificial intelligence and machine learning techniques in medical imaging analysis have brought new hope for disease diagnosis²⁶. By enhancing diagnostic accuracy, automating clinical workflows, and enabling personalized treatment strategies, AI technologies are profoundly transforming the landscape of intelligent healthcare²⁷. In DLBP diagnosis, the integration of radiomics with machine learning algorithms holds great potential to improve diagnostic accuracy and efficiency²⁸. However, the "black-box" nature of complex machine learning models limits their clinical applicability, making explainable AI (XAI) a critical component of clinical decision support systems²⁹. SHapley Additive exPlanations (SHAP), a game theory-based interpretive method, quantifies each feature's contribution and impact on model diagnostic predictions, making complex machine learning models more transparent and interpretable³⁰. It has become one of the most widely used XAI techniques in clinical decision support systems, effectively enhancing clinicians' trust and acceptance of AI systems²⁹.

This study developed a Random Forest algorithm incorporating radiomics and SHAP interpretability analysis, significantly improving the diagnostic efficiency of DLBP compared to traditional subjective assessments. Furthermore, the SHAP interpretation method enhances the transparency of the model's diagnostic predictions, boosting clinicians' confidence in the diagnosis and improving the accuracy of clinical decision-making. This study posits that this model will reduce unnecessary suffering caused by misdiagnosis and delayed treatment, promote the advancement of personalized treatment, and ultimately improve the overall treatment outcomes for patients.

Protocol

Ethical considerations and study population
This retrospective study was approved by the institutional ethics committee. Informed consent was waived as all protected health information was anonymized. The study population consisted of patients retrieved from the imaging database of the First People's Hospital of Nantong, who underwent lumbar spine MRI due to low back pain between January 2022 and December 2023. Clinical features of the patients were also collected (Table 1).

Inclusion and exclusion criteria
The current diagnostic standard for DLBP follows the 1995 International Association for the Study of Pain discography method, which provokes pain through increased pressure but is not widely accepted due to its high invasiveness⁸. Therefore, this study adopted a restrictive intraoperative disc pain provocation test as the diagnostic method for DLBP, consistent with discography principles while minimizing further disc damage and complications³¹. The procedure involved patients in a prone position under spinal anesthesia, with the target lumbar disc located using C-arm X-ray guidance. An 18 G (<22 G) needle was inserted via a posterolateral approach into the central nucleus pulposus, avoiding nerve roots and the dural sac. After confirming needle placement, saline or non-ionic contrast (e.g., iohexol) was injected at a rate not exceeding 0.5 mL/min, with pressure below 50 psi and a total injection volume of less than 3 mL, to simulate increased intradiscal pressure and provoke familiar low back pain symptoms, with a visual analog scale (VAS) score ≥ 7.

DLBP group inclusion criteria: The patients included here underwent MRI examination, have recurrent low back pain for over 3 months, with failed conservative treatment, with or without lower limb numbness or radiating pain, and are positive for the intraoperative disc pain provocation test.

Non-DLBP group inclusion criteria: Patients included here underwent MRI; have no history of low back pain within 3 months or are healthy individuals undergoing physical examination; have no MRI abnormalities; and have a standardized assessment with Oswestry Disability Index (ODI) <10 and VAS ≤2.

Exclusion Criteria: Exclude patients who have other causes of low back pain, such as significant disc herniation compressing nerves, fractures, spinal infections, spondylolisthesis, tumors, osteoporosis, or metabolic bone diseases; history of surgery prior to examination; unclear or poor-quality images; inability to identify the region of interest (ROI).

Ultimately, 243 patients were included, comprising 81 DLBP patients and 162 controls.

MRI parameters
All patients included in this study underwent 3.0 T MRI scans, using sequences that included sagittal and axial T1-weighted imaging (T1WI) and T2-weighted imaging (T2WI). Three different MRI machines were used in the study: Siemens Verio, Siemens Prisma, and Philips Ingenia CX. The scan parameters were set as follows: for sagittal T2WI, TR ranged from 2000 to 4597 ms, TE from 90 to 120 ms, slice thickness from 4.0 to 4.8 mm, with 15 slices included, bandwidth from 250 Hz to 340 Hz, matrix size of 384 × 384 or 512 × 512, phase field of view percentage of 100%, and a readout field of view of 300 mm.

HIZ imaging features measurement and statistics
This study used 3D Slicer software (version 5.6.1, https://download.slicer.org/?version=5.6.1) to manually analyze lumbar MRI T2-weighted images from 243 patients and controls. The data processing workflow was as follows: Images were imported using the Add DICOM Data function in 3D Slicer. Measurements were performed using the Markups function by two spine research graduate students under the guidance of radiologists and spine surgeons with over 10 years of clinical experience. Discrepancies in measurements were resolved through consultation with the two doctors. All measurement values were averaged by two statisticians to ensure accuracy. On sagittal images, the two mutually perpendicular maximum diameters of the high intensity zone (HIZ) were measured, defined as Hizh (near-vertical direction) and Hizw (near-horizontal direction). The Segmentation function was used to delineate the HIZ region, and the area of the most prominent HIZ on the sagittal plane was calculated, denoted as Hizarea. On axial images, the maximum length of the HIZ was measured, named Hizl. Additionally, three binary variables were defined: Hiz (presence of HIZ), Other (presence of multisegmental HIZ), and Position (whether HIZ crosses the posterior midline) (Table 2).

Radiomics feature extraction and standardization
To integrate the specific operation of Resampling Scalar Volume into the context of processing lumbar MRI T2-weighted images from 243 patients and controls using 3D Slicer, the following detailed steps were followed: The process began with importing lumbar MRI T2-weighted images into 3D Slicer. To ensure consistency and reduce heterogeneity bias, all images were resampled to a voxel size of 0.6 × 0.6 × 0.6 mm using the Resample Scalar Volume module. The specific steps were: Navigate to the Modules section and select Resample Scalar Volume. Under Parameter Set, ensure Resample Scalar Volume is selected. In the Resampling Parameters section, set Spacing to 0.6, 0.6, 0.6 to define the target voxel dimensions. Select an appropriate Interpolation method from options such as linear, nearest neighbor, bspline, hamming, cosine, welch, lanczos, or blackman, with linear as the default. For Output Volume, select or create a new volume to store the resampled data. After verifying settings, click Apply to execute the resampling process.

After resampling, two spine research graduate students, guided by radiologists and spine surgeons with over 10 years of experience, performed semi-automatic delineation of the region of interest (ROI). Discrepancies were resolved through consultation. The specific steps were: Select the layer (resampled new volume), create a new segmentation layer, and use the Draw function for layer-by-layer delineation. Use the Fill between slices function for interlayer filling. Smooth the formed ROI using the Median smoothing method with a kernel size of 3.0 mm, 5 x 5 pixels. Using the PyRadiomics library, 107 radiomic features were extracted from the resampled images, including shape features, first-order statistical features, gray-level co-occurrence matrix (GLCM), gray-level run-length matrix (GLRLM), gray-level size zone matrix (GLSZM), neighboring gray-tone difference matrix (NGTDM), and gray-level dependence matrix (GLDM) features (Supplementary Table 1). Finally, the extracted feature data were standardized using the Z-score method to transform their natural range into a standardized range.

Grouping based on features
The groups were assigned as follows:
d0 (Base group): Included clinical features, n(d0) = 5.
d1 (Base finetune group): Included clinical features and HIZ imaging features, n(d1) = 12.
d2 (Model group): Included clinical features and radiomics features, n(d2) = 112.
D (Model finetune group): Included clinical features, HIZ imaging features, and radiomics features, n(D) = 119.

Data reading and preprocessing
Data for the respective groups were read using the read_excel function in R software (version 4.3.1, https://www.r-project.org/, platform: x86_64-w64-mingw32/x64 [64-bit]) with UTF-8 encoding (system default), which supports most written languages globally. The select function was used to separate the target variable from the feature variables.

Feature selection
To ensure reproducibility of data splitting, the four datasets (D, d2, d1, d0) followed the same processing workflow, with a fixed random seed of 80. A label column was added to the data frame, converted to a factor type, and random labels (1 and 0) were generated using a binomial distribution to split the data into training and test sets in an 8:2 ratio (80% probability for the training set). Due to the limited number of features, d0 and d1 groups did not require feature selection, while d2 and D groups underwent feature selection:

First, the Mann-Whitney U test was applied to the training set to select features with p-values <0.05. Then, Lasso regression with 10-fold cross-validation and L1 regularization was performed to determine the optimal penalty coefficient and select features from a single iteration. A random seed from 1 to 100 was set, and the above steps (data splitting, Mann-Whitney U test, and Lasso regression) were repeated in a loop. Features appearing more than 50 times in 100 iterations were selected for subsequent modeling to ensure reliability in model construction and diagnostic prediction. Modeling features for each group are listed in Table 3.

Grid search and model tuning
To optimize the models, grid search was applied to systematically explore and evaluate different hyperparameter combinations. Through iterative tuning and performance assessment, this study identified the set of hyperparameters that provided the best ROC AUC values for the train and test sets, thereby optimizing the model's performance.

Model development and assessment
All model development and data analysis were conducted in R (version 4.3.1). This study applied various machine learning algorithms, including Random Forest (RF), Support Vector Machine (SVM), Decision Tree (TREE), K-Nearest Neighbors (KNN), and Logistic Regression (LOG), to develop 20 models across the four groups (d0, d1, d2, and D) for predicting DLBP diagnosis. Model performance was evaluated using the following metrics: ROC AUC, PR AUC, accuracy, sensitivity, specificity, positive predictive value (PPV), negative predictive value (NPV), and F1 score.

The ROC AUC and PR AUC values were directly generated from the ROC and PR curves, while the remaining metrics (accuracy, sensitivity, specificity, PPV, NPV, and F1 score) were calculated and statistically analyzed using corresponding functions in R. Among the 20 developed models, this study selected 8 representative models for assessment based on the above performance metrics.

The formulas for each metric are as follows:

Q

TP: True Positive; TN: True Negative; FP: False Positive; FN: False Negative.

SHAP interpretability analysis
The Random Forest model in the D group demonstrated the best performance in model assessment. To further interpret the model's predictions, this study rebuilt the model in Python (version 3.7.9) and performed SHAP interpretability analysis. To maintain consistency, the train hyperparameters of the Random Forest model in Python were identical to those used in R, and the data processing workflow was also consistent. Specifically, we used the tree model interpreter to analyze the contribution of each feature to the model's diagnostic predictions in the train and test sets. For better visualization, we generated SHAP value distribution plots, feature importance plots, and SHAP force plots for individual predictions.

Statistical analysis
All data analyses were conducted using R (version 4.3.1, https://www.r-project.org/, platform: x86_64-w64-mingw32/x64 (64-bit)) and Python (version 3.7.9, https://www.python.org/downloads/release/python-379/). Continuous variables were described as mean ± standard deviation, while categorical variables were expressed as frequency and percentage. The comparison between DLBP and Non-DLBP was analyzed using the U-test, with a P-value of less than 0.05 considered statistically significant.

Results

This study was based on MRI T2WI data from 81 surgically confirmed DLBP patients and 162 controls, totaling 243 patients. Radiomic features were extracted and combined with clinical and HIZ imaging features, with data split into training and test sets in an 8:2 ratio to build predictive DLBP diagnostic models. Specifically, seven HIZ imaging features were obtained after measuring the target disc using 3D Slicer and were included in the d1 and D group models along with clinical features. The 107 radiomic features extracte...

Discussion

The D group model, integrating clinical features, HIZ imaging features, and radiomic features, stabilized at 15 key features through feature selection and multiple iterations, demonstrating superior average performance compared to other groups, with higher stability and predictive ability. The random forest model outperformed other models in all groups (except d0), particularly in the D group, with ROC AUC values close to 1 and minimal differences between training and test sets, establishing it as the best model. The d0 ...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
3D Slicer 5.6.1	3D Slicer	https://download.slicer.org/?version=5.6.1
Data Storage	Zenodo	https://doi.org/10.5281/zenodo.17365220.
GeForce RTX 3060 Laptop GPU	NVIDIA	N/A
Ingenia CX 3.0T MRI machine	Philips	N/A
Prisma 3.0T MRI machine	Siemens	N/A
Python 3.7.9	Python	https://www.python.org/downloads/release/python-379/
R 4.3.1	R	https://www.r-project.org/
Verio 3.0T MRI machine	Siemens	N/A