Prognostic factors identification using machine learning in temporomandibular disorder treatment responders

2.1 Study design and population

This retrospective cohort study was conducted through a review of electronic medical records retrieved from the Hospital Information System for patients diagnosed with TMD and treated at the Orofacial Pain and Dental Sleep Medicine clinic, Faculty of Dentistry, Khon Kaen University, between 01 January 2018 and 30 June 2023. The inclusion criteria were: (1) patients diagnosed with TMD based on the Diagnostic Criteria for Temporomandibular Disorders (DC/TMD) [13] by qualified dental professionals, using the International Classification of Diseases, Tenth Revision (ICD-10) [14]; (2) patients who initiated treatment during the defined study period; and (3) patients with available treatment records for a minimum of one year of follow-up, with final documentation completed no later than 30 June 2024. The exclusion criteria were: (1) patients presenting with concurrent neuropathic or idiopathic orofacial pain conditions such as trigeminal neuralgia or burning mouth syndrome; (2) individuals with TMD secondary to movement disorders, muscular hypertrophy, fractures, or soft tissue lesions; and (3) cases with incomplete demographic or clinical data or patients who were lost to follow-up or could not establish further contact.

2.2 Data collection

A total of 1237 patient records were initially retrieved based on the following ICD-10 diagnostic codes: M26.62 (TMJ pain), K07.60 (TMD), K07.61 (TMJ disc-condyle complex disorders), K07.62 (TMJ subluxation), M19.9 (degenerative joint disease), M79.0 (myofascial pain), R25.2 (muscle cramp and spasm), and M27.8 (coronoid hyperplasia). After applying the predefined inclusion and exclusion criteria, a final cohort comprising 1050 eligible patients was established. Clinical information was systematically extracted and manually text-mined from the electronic medical records corresponding to the baseline visit only. All data were subsequently de-identified using a standardized case record form created in Microsoft Excel 365 (Version 2406, Microsoft Corp., Redmond, WA, USA). Initially, 26 demographic and clinical variables were collected. Categorical variables were processed using one-hot encoding, which expanded the number of variables to 36. Variables demonstrating high intercorrelation (correlation coefficient >0.7) were removed to reduce multicollinearity and improve the robustness of subsequent predictive analyses. All demographic, clinical, diagnostic, and outcome variables that were collected and analyzed are summarized in Table 1 (Ref. [4, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]).

2.3 Data preprocessing

Prior to analysis, all collected data were subjected to exploratory data analysis (EDA) to evaluate variable distributions and detect patterns within the dataset. Bivariate analyses were conducted to identify correlations among variables, highlight key features, and assess potential regression relationships. Descriptive statistics were calculated to determine the extent to which variables conformed to a normal distribution. Outliers and missing values were identified during this process and were subsequently excluded from further analysis. Data attributes were transformed into appropriate analytical formats. Specifically, numeric features were standardized using a standard scaler (min–max scaling), while categorical and mixed-type variables were processed using one-hot encoding. For inferential statistical analysis, the Chi-square test was employed to compare categorical variables, and the Mann-Whitney U test was used for continuous variables. A two-sided p-value < 0.05 was considered to indicate statistical significance.

2.4 Machine learning algorithms, model development, and evaluation

A final dataset comprising 1050 patients was included in the analysis. Prior to model development, variables exhibiting high intercorrelation (Pearson’s correlation coefficient >0.7) were removed to minimize multicollinearity and enhance predictive accuracy.

A nested cross-validation (Nested CV) framework was implemented to ensure robust model evaluation and generalizability. This framework involved 10 independent random trials. The outer loop employed a 4-fold K-Fold cross-validation strategy, partitioning the dataset into 75% training and 25% testing subsets. Within each outer fold, an inner cross-validation loop was used to perform hyperparameter optimization via GridSearchCV (scikit-learn library [26]), with the objective of maximizing the receiver operating characteristic-area under the curve (ROC-AUC) on the training set.

Data imbalance was addressed using built-in balancing mechanisms specific to each algorithm, as well as through hyperparameter tuning strategies that incorporated the synthetic minority oversampling technique (SMOTE) and class-weighting adjustments. Seven machine learning algorithms were developed and evaluated: Decision Tree, Random Forest, Logistic Regression, K-Nearest Neighbors (KNN), Support Vector Classifier (SVC), Gradient Boosting, and Extreme Gradient Boosting (XGBoost). Model performance was assessed using multiple metrics, including accuracy, precision, recall, F1 score (micro and macro), and ROC-AUC.

For comparative model analysis, receiver operating characteristic (ROC) curves were generated for each algorithm. Mean true positive rates (TPRs) and false positive rates (FPRs) were calculated and averaged across all trials. The resulting aggregated mean ROC curves and their standard deviations were plotted to allow visual comparison of each model’s predictive performance. Furthermore, normalized confusion matrices were computed and visualized to evaluate the classification accuracy in greater detail.

To assess the relative importance of predictive features, two interpretability methods were employed: information gain and SHAP plots. For information gain, feature importance scores were computed by averaging the mean results obtained from 10 repeated trials. SHAP analysis was performed on the optimal model using the optimal hyperparameter configurations identified during model tuning, thereby offering interpretable insights into the contribution of each feature to the model’s predictions. Three types of SHAP plots were used to present the results: (1) bar plots showing the mean absolute SHAP value for each feature across all predictions, indicating the average magnitude of impact and enabling feature ranking; (2) SHAP summary plots (beeswarm plots), which depicted both the magnitude and direction of each feature’s effect; and (3) heatmaps that provided a visual overview of SHAP values for individual predictions across multiple features [12, 27].

2.5 Software and tools

Machine learning algorithms were implemented using the Python programming language (version 3.10.12, Python Software Foundation, Wilmington, DE, USA) within the Google Collaboratory notebook environment (Version 4, California, USA) [26]. All electronic devices utilized during the research process, including computers and storage drives, were secured with password protection. Data encryption protocols were employed to ensure the protection of sensitive information in accordance with standard data security practices, thereby maintaining participant confidentiality throughout the study.

3.1 Patient demographics and clinical characteristics

From a retrospective review of 1237 patients newly diagnosed with TMD, 65 patients were excluded due to incomplete data, 49 patients were excluded based on the presence of other comorbid pain conditions, and 73 patients classified as improperly discharged were also excluded. This resulted in a final dataset of 1050 patients eligible for EDA, model development, and evaluation, as illustrated in Fig. 1.

Fig. 1.Flowchart of all TMD patients for model analysis. EDA: Exploratory Data Analysis; TMD: Temporomandibular Disorders.

The mean age of the included patients was 33.02 ± 16.67 years. The cohort exhibited a female predominance, with 779 females (74.19%) and 271 males (25.81%). Additional clinical and demographic characteristics, stratified by discharge status, are summarized in Table 2. Quantitative variables are presented as mean ± standard deviation (SD), while categorical variables are reported as frequencies and percentages. Given the non-normal distribution of the data, the Mann-Whitney U test was applied for comparisons involving continuous variables, and the Chi-square test was used for categorical variables, as shown in Table 2.

3.2 Correlation and multicollinearity analysis of variables

EDA revealed correlations among both continuous and binary variables in the dataset. As illustrated in Fig. 2A, two prominent correlations were observed among the continuous variables: anxiety and stress levels (correlation coefficient = 0.78), and anxiety and depression levels (correlation coefficient = 0.68). Among the binary variables, shown in Fig. 2B, the strongest association was identified between back pain and neck pain (correlation coefficient = 0.51). Due to the high correlation between stress and anxiety levels (correlation coefficient >0.7), the stress variable was excluded from further analysis to minimize redundancy. This resulted in a refined dataset consisting of 36 variables for model development. Despite these observed correlations, multicollinearity was not deemed problematic. Variance inflation factor (VIF) analysis indicated that all variables remained within acceptable limits, with the highest VIF value being 3.59 for the neck pain variable. Since VIF values exceeding 10 are typically indicative of problematic multicollinearity, the results confirm that multicollinearity was within acceptable thresholds in the dataset.

Fig. 2.Correlation matrix of variables. The gradient color scale represents the strength of the correlation coefficients for (A) continuous variables and (B) binary variables. Stronger correlations are indicated by darker shades. BMI: Body Mass Index; NRS: Numeric Rating Scale.

3.3 Evaluation of model performance

The performance metrics of all machine learning algorithms are summarized in Table 3. Model evaluation primarily relied on ROC-AUC, due to its robustness in assessing classifier performance in the presence of class imbalance. Among the evaluated models, the Random Forest algorithm demonstrated the highest overall predictive performance based on ROC-AUC. It achieved an accuracy of 0.7901 ± 0.0383. The precision for predicting undischarged and discharged patients was 0.3512 ± 0.1239 and 0.8611 ± 0.0193, respectively. Recall values were 0.2459 ± 0.0891 for the undischarged group and 0.8943 ± 0.0549 for the discharged group. The model yielded an F1 score (micro) of 0.7901 ± 0.0383 and an F1 score (macro) of 0.5709 ± 0.0254. Class-specific F1 scores were 0.2654 ± 0.0584 for the undischarged group and 0.8764 ± 0.0263 for the discharged group. The corresponding ROC-AUC was 0.6442 ± 0.0455. In contrast, the KNN model exhibited the lowest ROC-AUC despite a relatively higher accuracy. It achieved an accuracy of 0.8336 ± 0.0170, with precision values of 0.4321 ± 0.2042 for the undischarged group and 0.8456 ± 0.0192 for the discharged group. However, recall for the undischarged group was notably low at 0.0655 ± 0.0338, compared to 0.9810 ± 0.0123 for the discharged group. The F1 score (micro) was 0.8336 ± 0.0170, and the F1 score (macro) was substantially lower at 0.5092 ± 0.0268. Class-specific F1 scores were 0.1103 ± 0.0519 for the undischarged group and 0.9081 ± 0.0102 for the discharged group. The model’s ROC-AUC was 0.5717 ± 0.0490.

The normalized confusion matrices presented in Fig. 3 illustrate the average classification performance of each machine learning model across 10 repeated experiments, highlighting their ability to differentiate between class 0 (undischarged group) and class 1 (discharged group). The Gradient Boosting and KNN classifiers showed excellent performance in identifying discharged patients, achieving true positive rates of 99% and 98%, respectively. However, both models demonstrated poor sensitivity in detecting the undischarged group, with correct classification rates of only 7%. The Random Forest model also performed well for the discharged group, correctly identifying 89% of those cases, while achieving a higher recall for the undischarged group at 25%. Notably, the XGBoost model demonstrated more balanced classification, with 73% of discharged patients and 45% of undischarged patients correctly identified.

Fig. 3.Normalized confusion matrices of all machine learning models for classifying discharged (class 1) and undischarged (class 0) TMD patients. The values represent the average of 10 repeated experiments for each model. Each matrix displays the proportion of correctly and incorrectly classified instances within each class, providing a visual summary of classification performance. KNN: K-Nearest Neighbors; SVC: Support Vector Classifier; XGBoost: Extreme Gradient Boosting.

In contrast to models with high-class imbalance, SVC exhibited more balanced classification performance across both outcome groups. Logistic Regression correctly identified 54% of patients in the undischarged group and 65% in the discharged group. Similarly, SVC achieved classification rates of 52% and 67% for the undischarged and discharged groups, respectively. The Decision Tree model demonstrated the weakest performance overall, correctly classifying only 50% of undischarged and 62% of discharged patients.

Fig. 4 presents the average ROC-AUC values for all machine learning models. Based on results from the nested CV framework, Logistic Regression yielded a mean ROC-AUC of 0.64 ± 0.04, while Random Forest achieved a comparable mean ROC-AUC of 0.64 ± 0.05, albeit with a slightly higher standard deviation. Among the models evaluated, these two demonstrated the most consistent and reliable predictive performance across the outer validation fold.

Fig. 4.Comparison of average receiver operating characteristic (ROC) curves for all machine learning models, derived from nested cross-validation (Nested CV). Each curve represents the mean performance across 10 repeated experiments. The area under the curve (AUC) and corresponding standard deviation are indicated in the legend. The diagonal line denotes the performance of a random classifier (AUC = 0.5). KNN: K-Nearest Neighbors; SVC: Support Vector Classifier; XGBoost: Extreme Gradient Boosting.

3.4 Feature important analysis

3.4.2 SHAP value of feature importance

The SHAP-based feature importance results for the Random Forest model are shown in Fig. 5B. In this bar plot, features with higher absolute SHAP values, reflecting greater average impact on model predictions, are ranked at the top [27]. The most influential feature was onset duration, which had the highest mean SHAP value of 0.05. This was followed by pain severity (mean SHAP = 0.03), anxiety level (mean SHAP = 0.03), presence of neck pain (mean SHAP = 0.02), and BMI (mean SHAP = 0.02). Positive SHAP values for these variables indicated that they contributed to classifying a patient into the discharged group. Notably, while these top-ranked features showed substantial individual importance, the aggregate contribution of the remaining lower-ranked features collectively accounted for the highest cumulative impact, with a combined mean SHAP value of 0.06. Furthermore, there was complete overlap between the top nine features identified by SHAP analysis and those ranked highly in the information gain method, which included anxiety level, onset duration, BMI, pain severity, age, depression level, presence of psychiatric problems, sleep problems, and neck pain.

Fig. 5.Visualization of feature importance in the Random Forest model. (A) Bar chart displaying the ranking of predictive features based on information gain scores. (B) Bar chart showing the mean absolute SHAP values for each predictor, indicating their relative contribution to classification into the discharged group. BMI: Body Mass Index; NRS: Numeric Rating Scale; SHAP: Shapley Additive exPlanations.

SHAP beeswarm plot shown in Fig. 6A illustrates the distribution and individual contributions of features to the predictive outcomes generated by the Random Forest model. Each dot represents a Shapley value for a single feature in one patient case. Feature names are presented on the y-axis, while the x-axis denotes the magnitude and direction of each Shapley value. The color gradient of each point reflects the actual feature value, ranging from low (blue) to high (pink), and the vertical dispersion highlights the variability of feature effects across the dataset. Features are ranked in descending order of overall importance [27]. Notably, shorter onset duration, lower pain severity, and reduced anxiety levels were predominantly associated with positive SHAP values, indicating a higher likelihood of discharge. In addition, the absence of neck pain and higher BMI, ranked fourth and fifth, respectively, also contributed positively to discharge predictions. The SHAP summary heatmap in Fig. 6B provides further insight into the influence of key predictors across individual patients. In this visualization, each row corresponds to a feature, while each column represents a single patient instance. The color intensity reflects the magnitude and direction of the SHAP value, with pink indicating a positive contribution to the prediction (i.e., discharged group) and blue indicating a negative contribution (i.e., undischarged group). This heatmap highlights the consistent importance of onset duration, pain severity, and anxiety level, which exhibit strong contributions across multiple cases [27]. In contrast, features such as neck pain and BMI showed consistent but more moderate influence on model predictions.

Fig. 6.SHAP value distribution and feature contributions in the Random Forest model for predicting discharge outcomes. (A) Beeswarm plot showing the direction and magnitude of each feature’s impact on the model output for the discharged group. Each dot represents a single patient; the horizontal axis indicates the SHAP value, while the color denotes the actual feature value (pink for high and blue for low). Features are listed in order of their overall importance. (B) SHAP summary heatmap displaying feature contributions to model predictions across individual patient cases. Each row corresponds to a feature, and each column to a patient. Pink indicates a positive SHAP value (favoring discharge prediction), while blue indicates a negative SHAP value (favoring undischarged prediction). BMI: Body Mass Index; NRS: Numeric Rating Scale; SHAP: Shapley Additive exPlanations.

4.1 Prognostic factors influencing TMD patient readiness for discharge

To our knowledge, this is the first study to investigate prognostic factors influencing treatment response based on readiness for discharge in patients with TMD managed at a tertiary orofacial pain clinic, utilizing advanced machine learning methodologies. By using SHAP value-based interpretation within the machine learning framework, we identified key predictors that contribute to a patient’s likelihood of achieving discharge readiness [28]. Previous studies have highlighted the role of psychological factors [9], pain intensity [29], and biopsychosocial factors [5] as significant predictors of treatment response in TMD patients. In our analysis, onset duration emerged as the most influential predictor across multiple machine learning models in determining patient readiness for discharge. Pain severity was identified as the second most important variable, likely reflecting its direct impact on functional status and overall recovery trajectory. Anxiety level, as an indicator of psychological well-being, also played a key role in capturing both the emotional and physical dimensions of TMD. Additionally, the presence of neck pain, which is frequently associated with orofacial and masticatory dysfunction, emphasized the interconnected nature of symptom presentation in TMD patients [23]. BMI also emerged as a relevant predictor, suggesting a potential association between lifestyle or systemic factors and the chronicity of TMD symptoms. In addition, the cumulative influence of numerous additional features, quantified through SHAP values, further reinforces the multifactorial and heterogeneous nature of TMD, which collectively shapes patients’ readiness for discharge.

4.2 Toward a holistic model of TMD management

In this tertiary care setting, discharge decisions were jointly determined by the attending clinician and the patient, and documented in the medical records. Criteria influencing discharge often included tolerable pain levels, meaningful improvement in jaw function, or the patient’s ability to adapt to residual symptoms such as TMJ clicking. This study aimed to capture the patient’s subjective response to treatment, while simultaneously evaluating objective clinical parameters as prognostic indicators. Although objective metrics have traditionally dominated assessments of treatment success or failure, emerging evidence suggests that subjective satisfaction may persist despite minimal change in quantifiable outcomes [38, 39].

The findings of this study underscore the multidimensional nature of TMD, wherein multiple factors, rather than any single variable, contribute to a patient’s readiness for discharge. From a neurobiological perspective, pain is shaped by both physiological and psychological processes, reinforcing the need for a comprehensive, biopsychosocial approach to clinical management [1, 23]. Historically, chronic pain treatment models have prioritized the resolution of primary symptoms, often at the expense of adequately addressing psychosocial dimensions during routine consultations [31]. However, there is an increasing number of studies that support the use of integrative care models that concurrently address emotional, social, and physical health domains within a single clinical encounter [3, 6]. A stronger doctor-patient relationship, coupled with a focus on functional well-being and psychological coping capacity rather than symptom severity alone, may foster more effective and collaborative chronic pain management. This patient-centered approach promotes self-awareness, resilience, and active engagement in self-management practices [30]. Adopting such a holistic framework in TMD care may thus contribute to more sustainable outcomes and improved quality of life.

Our findings indicated that 16.10% of patients experienced persistent TMD, preventing them from being discharged, suggesting that most (i.e., nearly 84%) were able to manage or recover from their symptoms, aligning with the fluctuating nature of TMD, in which pain-free intervals often exceed periods of active symptoms [2]. As a result, any subsequent reduction in symptoms may partly reflect natural fluctuations rather than treatment effects, a phenomenon known as regression to the mean. This interpretation is further supported by epidemiological studies, which report that 15–20% of TMD cases exhibit severe symptoms as measured by the Graded Chronic Pain Scale, while only approximately 10% of adolescents experience worsening pain over time [3]. These findings suggest that TMD symptoms often resolve when patients are able to clearly identify symptom triggers and engage in effective self-care strategies [3, 9]. A key element in this process is self-confidence, as evidenced by studies demonstrating that clear explanations of TMD pathophysiology, combined with the identification and reduction of aggravating factors, can alleviate symptom burden [6]. Therefore, by providing patients with knowledge and reinforcing self-belief, clinicians may facilitate more effective self-management and increase the likelihood of long-term recovery.

4.3 Enhancing model performance

Some models in this study, including Random Forest, KNN, and Gradient Boosting, demonstrated high accuracy but only moderate ROC-AUC values. This discrepancy is commonly observed in the context of imbalanced datasets, where prediction bias favors the majority class. A major challenge encountered in the machine learning pipeline was overfitting, which reduced the generalizability and predictive accuracy of certain models. Specifically, Gradient Boosting and KNN showed a pronounced tendency to favor predictions toward the discharged group, likely due to class imbalance or overfitting to the dominant class. While these models performed well in identifying discharged patients, they frequently misclassified undischarged cases—a pattern similarly noted with the Random Forest model.

In contrast, the XGBoost classifier yielded a more balanced performance between groups, although its overall predictive capacity remained suboptimal. Logistic Regression and SVC provided the most balanced classification results, demonstrating less bias and potentially greater suitability in scenarios where equal attention is required for both outcome groups. Meanwhile, the Decision Tree model exhibited limited discriminative power, possibly due to overfitting or inadequate model complexity. To improve model performance, several strategies were employed. Hyperparameter tuning was applied to optimize model-specific settings; for example, varying the depth of tree-based models allowed identification of the appropriate level of complexity to avoid overfitting [26]. Another essential technique was feature selection, which was used to eliminate irrelevant or redundant variables, thereby reducing noise and enhancing overall accuracy [9]. Finally, managing data imbalance was essential to ensuring reliable predictions. Techniques such as oversampling, undersampling, or the use of algorithms designed to address imbalance, such as SMOTE or class-weight adjustments, were effective in reducing the risks of biased predictions. In this study, class-weight adjustment was implemented to account for data imbalance, contributing to more reliable and equitable model performance.

4.4 Feature importance illustrations using information gain or SHAP

Information gain analysis ranks features based on their capacity to reduce uncertainty in predicting the target variable. A higher information gain indicates a greater contribution of that feature to the prediction, and the metric is calculated using entropy, which measures the level of uncertainty within the dataset. However, a known limitation of this approach is its tendency to favor features with many unique values, which may appear predictive even if their actual contribution is minimal. In addition, in datasets containing noise, entropy estimates may become unstable, potentially leading to suboptimal decision splits [43].

To overcome these limitations and obtain a more nuanced understanding of feature contributions, SHAP analysis was also applied. SHAP is grounded in the cooperative game theory and provides contrastive explanations by comparing individual predictions to the average model output, distributing attribution fairly among input features [44]. One of SHAP’s strengths lies in the interpretability of its values, which are expressed in the same units as the model output. SHAP facilitates both global interpretation, through feature importance rankings, dependence plots, interaction effects, clustering, and summary plots, and local interpretability at the level of individual predictions. Despite these advantages, SHAP also has limitations. Calculating exact Shapley values across many instances can be computationally demanding, and global SHAP analyses require aggregation over a large number of samples to ensure stability. Moreover, in the presence of multicollinearity, SHAP may assign disproportionately high importance to statistically improbable combinations of features, potentially complicating interpretation [27].

4.5 Strengths and limitations

A major strength of this study was the setting in a tertiary care center, which enabled access to a substantial number of patients with chronic TMD. This environment facilitated specialist-led clinical assessments, thereby ensuring a higher level of diagnostic accuracy compared to studies relying solely on self-reported symptoms. The application of machine learning further enhanced the study’s analytical depth, allowing for the identification of complex relationships between variables that may extend beyond conventional expert assumptions while reducing subjective bias. Additionally, machine learning techniques are often superior to traditional statistical methods in capturing non-linear associations, thus providing a more nuanced understanding of the dataset.

However, this study also has several limitations. First, the findings may not be generalizable to broader TMD populations, particularly those of non-Asian ethnicities. Furthermore, patients in tertiary care settings often present with more complex conditions and higher expectations for treatment, which may not reflect the clinical profiles typically seen in primary care. The inherent class imbalance in TMD outcomes, where most patients are ultimately discharged, may also bias the predictive models toward the majority class, thereby resulting in moderate overall model performance. Nonetheless, the predictors identified in this study may still serve as a useful screening tool for general practitioners to identify potentially complicated cases. These findings are consistent with those of the prospective OPPERA study, which similarly investigated prognostic indicators in persistent TMD patients [7]. The retrospective design posed additional challenges, including susceptibility to recall bias, variability in documentation quality, and missing data, all of which may compromise data integrity. Moreover, because formal discharge criteria were not uniformly established, discharge status in this study was determined based on documentation in the medical records, incorporating both professional judgment and patient preference. Future prospective or multi-center studies utilizing standardized data collection protocols would help mitigate these limitations and improve the external validity of TMD research.

4.6 Future recommendations

Future studies should aim to improve both the generalizability and predictive performance of the models by expanding the dataset to include patients from primary care settings, thereby capturing a more heterogeneous patient population and broader treatment contexts. Additionally, increasing the representation of non-discharged patients is essential to address class imbalance and enhance model discrimination. Clinical variables with stronger associations to discharge readiness, such as structure-specific pain (e.g., myalgia or arthralgia) and the 20-item checklist of general symptoms, should be incorporated, as supported by findings from the prospective OPPERA study [7]. Model performance may also benefit from a two-stage approach that integrates expert clinical judgment with the best-performing, well-balanced algorithm, such as the Logistic Regression model.