Interpretable Machine Learning to Predict the Malignancy Risk of Follicular Thyroid Neoplasms in Extremely Unbalanced Data: Retrospective Cohort Study and Literature Review

Introduction

Globally, thyroid neoplasms are becoming increasingly prevalent [1]. Among them, follicular thyroid neoplasms (FTNs) represent a major type but have garnered significantly less attention compared to papillary thyroid carcinoma. A key challenge is that over 95% of FTN cases cannot be reliably distinguished as benign (follicular thyroid adenoma [FTA]) or malignant (follicular thyroid carcinoma [FTC]) until diagnostic surgery [2]. This uncertainty often leads to both over- and undertreatment of patients with FTN. On one hand, it is estimated that over 80% of patients who undergo thyroidectomy might ultimately be diagnosed as benign FTN based on postoperative pathology [3]. On the other hand, those with malignant FTN may have already developed distant metastases to the lungs, bones, or other organs by the time they receive surgical treatment.

Several guidelines advocate for enhanced screening, accurate diagnosis, and appropriate treatment for patients with FTN [4,5]. One crucial solution is to develop prediction models to aid clinical decision-making for these patients. To date, machine learning has been proven effective in constructing predictive models for various cancers such as oral, gastrointestinal, and breast cancers [6-8]. Our literature review also indicated that machine learning technology excels at capturing complex, nonlinear relationships and high-dimensional intercorrelations among predictors [9-12].

However, our literature review revealed several limitations among most of the existing studies, mainly including (1) small sample sizes ranging from 18 to 888 participants [13-19], (2) the ratio of FTA to FTC deviating from the real population distributions, (3) reliance on simple linear models unable to capture the complex nonlinearity or interactions underlying predictor-outcome relationships [14,18], (4) using inappropriate metrics to evaluate model performance for the unbalanced data [13-19], (5) lack of assessing the extent to which a predictor influences the model’s prediction (ie, model interpretability) [13,16,19], (6) not evaluating whether the predicted probabilities were consistent with actual outcomes (ie, model calibration) in the development and validation of clinical prediction models [20], and (7) predictors are predominantly confined to sonographic features with limited consideration of other factors such as the presence of Hashimoto thyroiditis (an autoimmune disease that may increase the risk for differentiated thyroid cancer [21]).

To address these limitations, our study has united a multidisciplinary treatment team for thyroid neoplasms and accumulated a cohort of over 1500 patients with FTN over the past decade [22]. This provided us a unique opportunity to develop and validate clinical prediction models to bridge the current research gaps in the field of FTN. Specifically, we aimed to individualize the clinical decision-making for patients with FTN by using interpretable machine learning to not only predict the malignancy risk of FTN but also identify the important predictors that might contribute to the prediction.

Methods

Study Design

This retrospective cohort study followed the suggestions of the TRIPOD (Transparent Reporting of a Multivariable Prediction Model for Individual Prognosis or Diagnosis) statement [23]. Figure 1 shows the framework of our study.

‎

Study Population

Our multidisciplinary research team included experts in the fields of epidemiology, surgery, pathology, ultrasound, and endocrinology. We conducted a retrospective cohort study at Peking University Third Hospital in Beijing, China, from January 2012 to September 2023. Eligible patients were those who underwent surgery and were pathologically diagnosed with FTA or FTC following the procedure. Patients were excluded if they did not undergo ultrasound examinations prior to surgery or if they had nodules classified as follicular tumors of uncertain malignant potential (UMPs). This exclusion was based on two considerations: (1) accurate diagnosis of FTNs required both an experienced pathologist and a complete biopsy sample. The key to distinguishing between benign and malignant FTNs was determining whether the tumor invaded the capsule. Tumors that invaded the capsule or blood vessels were classified as FTC, while those that did not were considered FTA. If the pathologist struggled to assess capsule invasion due to inexperience, or if the sample was inadequately collected during surgery, leading to capsule damage, the tumor might not have been accurately classified as either FTA or FTC. In such cases, it could have been labeled as a UMP. (2) Through a literature review, we found that all previous research had excluded UMPs [13-19]. Therefore, our study also excluded UMPs, enhancing comparability with prior research. We paid close attention to the accuracy of the pathological diagnosis of FTN due to its high professional requirements, which include not only complete sampling but also a thorough examination of all areas of the tumor margin. To ensure this, we invited pathologists with expertise in thyroid tumors to double-check all the postoperational pathological diagnoses in the study population, based on the most recent 2022, 5th edition WHO Classification of Thyroid Neoplasms [24].

It is important to note that our study population reflected the natural distribution of FTNs (ie, the ratio of FTA and FTC), resulting in imbalanced data, with 84.57% (n=1414) of cases was FTA. Specifically, we did not restrict the ratio of FTA to FTC to 1:1 or any other fixed ratio in the main analyses. This approach allowed the results from the developed prediction model to be more readily applicable to external populations with a similar natural distribution.

Data Sources and Processing

The data for this study were sourced from the electronic health records, extracted by professional information management personnel from the hospital’s electronic information system. For critical data sources like thyroid pathology and neck ultrasound reports, a tailored data extraction form was designed using EpiData (EpiData Association), aligning with the study’s research questions. The form was iteratively refined through discussions among researchers, surgeons, ultrasound specialists, and pathologists, followed by trial entries and revisions until finalized. Trained clinical doctors and medical students performed manual data entry, with researchers conducting 2 rounds of random checks to ensure accuracy and consistency. Senior doctors performed a final review to verify data quality. Missing values were imputed using the mean for continuous variables and the mode for categorical variables.

Predictors

We selected the predictors to develop the machine learning–based model based on our systematic review [17,18], domain knowledge [25-27], and data available. The predictors included sonographic features, patients’ age, sex, BMI, whether or not diagnosed as Hashimoto thyroiditis (an autoimmune disease that destroys thyroid cells by cell and antibody-mediated immune processes [28]), and measurements of thyroid hormones. Specifically, the sonographic features included mean diameter, composition (solid, predominantly solid, predominantly cystic, or cystic), echogenicity (hyperechoic, isoechoic, hypoechoic, or anechoic), taller-than-wide (the length in the vertical direction is greater than the width in the horizontal direction: absent or present), margin (circumscribed, ill-defined, irregular, or lobulated), calcifications (microcalcifications, macrocalcifications, peripheral calcifications, punctate echogenic foci of undetermined significance, microcalcifications with comet-tail artifacts, or no echogenic foci), halo (absent halo, even thickness halo, uneven thickness halo, or present halo without evenness of thickness reported), internal blood flow (absent or present), vascularity (mainly central vascularity, mainly peripheral vascularity, mixed vascularity, or avascularity), trabecular formation (typically appears as elongated, band-like, or fibrous echogenicity, arranged in a reticular or cord-like pattern: absent or present), and nodule-in-nodule appearance (a smaller nodule or an area with different echogenic characteristics is present within a larger thyroid nodule: absent or present); the measurements at the latest examination of thyroid hormones included thyroid-stimulating hormone (TSH), free triiodothyronine, and free thyroxine; additionally, TSH-related features derived from all examinations included mean TSH score (interval-adjusted detailed TSH score) [29], time-adjusted root mean square of successive differences of TSH [30], mean TSH (mean value of preoperative TSH), and coefficient of variation of TSH (the ratio of SD of preoperative TSH to mean value of preoperative TSH), and detailed definitions were introduced in previous publication [31]. All selected predictors were carefully checked by both clinicians and researchers to ensure the accuracy and reliability of the study results.

Development and Validation of Machine Learning–Based Models

We established the machine learning–based model as shown in Figure 1. We selected features, trained models, tuned hyperparameters, and validated models, as briefly described below. We used the mlr3 [32] ecosystem in R (version 4.3.3; R Foundation for Statistical Computing), scikit-learn [33], and Shapley Additive Explanations (SHAP) [34] in Python (version 3.11.1; Python Software Foundation) to conduct machine learning.

Feature selection, which aims to reduce the number of features, offers several benefits including minimizing overfitting, enhancing model robustness, and accelerating predictions. Notably, it is particularly advantageous for datasets with a high feature-to-sample ratio, where the number of features exceeds the limited size of data points. To identify a core set of predictors that could effectively predict the outcome without redundancy, we used a novel information-gain approach for feature selection [35]. To finalize the optimal model, we also compared model performance between that with full predictors and that with selected predictors.

We trained 5 classification models including logistic regression, least absolute shrinkage and selection operator (LASSO) regression, random forest, extreme gradient boosting, and support vector machine. We comprehensively considered and weighed (trade-off) the performance, calibration, parsimony, and interpretability of models and selected the most appropriate one as our prediction model.

The random search and cross-validation were combined to select model hyperparameters when training the machine learning model. We performed a random search over more than 45,000 hyperparameter combinations to select the best hyperparameter combination and trained the final classifiers. Additionally, to address the issue of imbalance, both oversampling (increasing the amount of minority class samples with producing new samples or repeating some samples) and undersampling (decreasing the amount of majority class samples) techniques were applied [36].

Evaluation of Model Performance

We evaluated the performance of the developed model in terms of discrimination (the ability of the model to distinguish between those with and without the outcome) and calibration (the consistency or agreement between the observed outcomes and predicted risks from the model). For discrimination, we first showed the confusion matrix including the numbers and percentages of true positive, true negative, false positive, and false negative. We then calculated both the threshold-free and threshold-sensitive metrics. Threshold-free metrics included the area under the receiver operating characteristic curve (AUROC) and the area under the precision-recall curve (AUPRC). While AUROC was a common performance metric for discrimination, AUPRC was considered more useful and informative for handling the unbalanced data in this study [37]. Threshold-sensitive metrics included sensitivity, precision, specificity, and accuracy. For calibration, we first plotted predicted risks (x-axis) against observed outcomes (y-axis) using a smoothed flexible calibration curve based on individual data. We also quantitatively assessed calibration using the calibration slope and calibration-in-the-large.

Interpretation of Model Prediction Results

First, we evaluated the feature importance (ie, the extent of the model depended on the feature) and the feature interaction by using the SHAP summary plot and SHAP interaction value dependence plot (see details in Multimedia Appendix 1) [34]. Second, we figured the partial dependence plots to visualize the direction of predictor-outcome associations, illustrate whether the risk of the outcome increased with a rise or decline in the predictor values, and assess whether this relationship is linear. Third, we plotted individual conditional expectations curves to explore potential modifiers that could influence predictor-outcome associations. Finally, we separated FTC into 2 groups based on whether they were correctly predicted and compared their characteristics. The significance of differences between the groups was tested using the Mann-Whitney U test, as the data did not follow to a normal distribution.

A Systematic Review of the Previous Studies

We conducted a systematic review of previous studies addressing similar topics. According to the PRISMA (Preferred Reporting Items for Systematic Reviews and Meta-Analyses) guidelines [38] (see details in Checklist 1), we searched PubMed, Web of Science, Embase, and IEEE Xplore using the terms “follicular thyroid cancer” and “predict” for papers published up to October 1, 2023 (see details in Multimedia Appendix 2). Eligible studies included those that established prediction models to distinguish FTC from FTA before operation with various preoperative predictors. We included studies with either deep learning models, machine learning models, traditional statistical models, or other relevant methodologies. Studies were excluded if fewer than 50% of patients had FTN or if the papers were not written in English.

We evaluated the Thyroid Imaging Reporting and Data System for Follicular Neoplasm (F-TIRADS) scoring criteria developed by Li et al [18] to differentiate between benign and malignant cases in our dataset. These criteria are based on 6 key features: mean diameter, composition, echogenicity, margin, calcifications, and trabecular formation. Each specific characteristic of these features is assigned a corresponding point value, and the total points across the 6 features indicate the risk level of FTC. For instance, a total score of 12 points or higher suggests an FTC risk exceeding 90% (refer to Figure S1 in Multimedia Appendix 1).

Ethical Considerations

This study was classified as human participant research and was reviewed and approved by the medical research ethics committee of Peking University Third Hospital (IRB00006761-M2023168). As a retrospective analysis, the study was granted a waiver for additional informed consent. During the data extraction process, strict confidentiality measures were implemented to ensure patient privacy and data security. All extracted data were anonymized, with any information that could directly identify patients being removed.

Results

Characteristics of the Study Population

Altogether, we included 1539 patients, 1409 of whom had solitary tumors, and 130 had more than 1 tumor. Thus, a total of 1672 tumors were included and divided into 2 pathological types: FTA (n=1414) and FTC (n=258). The characteristics of the included tumors are listed in Table 1, and the characteristics of the study population are listed in Table 2. The age of the included population was 47.98 (SD 14.15) years (n=1530; missing value=9), the mean BMI was 24.18 (SD 3.66) kg/m² (n=1475; missing value=64), and the female population made up 73.16% (n=1126; male: n=342; missing value=71) of all patients.

Model Performance in Discrimination and Calibration

We compared performance among 5 models (logistic regression, LASSO regression, random forest, extreme gradient boosting, and support vector machine) using the AUROC and AUPRC. As shown in Multimedia Appendices 3 and 4, the random forest model performed better in both AUROC and AUPRC than the other 4 models. With comprehensive consideration of the discrimination, calibration, parsimony, and interpretability of models, we selected the random forest model as the optimal.

We developed a random forest model with a total of 22 features: age, sex, BMI, Hashimoto thyroiditis, thyroid hormones (TSH, free triiodothyronine, and free thyroxine), ultrasonic predictors (mean diameter, composition, echogenicity, taller-than-wide, margin, calcifications, halo, internal blood flow, vascularity, trabecular formation, and nodule-in-nodule appearance), and TSH-related variables (mean TSH score, time-adjusted root mean square of successive differences of TSH, mean TSH, and coefficient of variation of TSH). After 5-fold cross-validation, the AUROC of the prediction model was 0.79 (95% CI 0.77-0.81) and the AUPRC was 0.40 (95% CI 0.37-0.44). When the threshold is gradually lowered from 50%, 40%, 30%, 20%, and finally to 10%, the accuracy, specificity, and precision decreased step by step while the sensitivity increased progressively (Multimedia Appendix 5). The calibration slope and calibration-in-the-large were 1.16 and 0.13, respectively (Multimedia Appendix 6).

In addition, we implemented both oversampling and undersampling techniques to handle the imbalance in our models. However, following oversampling, the AUROC and AUPRC were 0.76 and 0.37, respectively, while after undersampling, the AUROC and AUPRC were 0.77 and 0.39, respectively. Notably, the model performed better before applying these sampling methods, with an AUROC of 0.79 and an AUPRC of 0.40.

Model Performance in Interpretation

The top 5 predictors were the mean TSH score, mean tumor diameter, mean TSH, coefficient of variation of TSH, and TSH level (Figure 2). The 5 top predictors did not show explicit interactions with the sex or other predictors (Multimedia Appendix 7). Due to the strong correlation between TSH level and mean TSH score (correlation coefficient>0.6), we only plotted the partial dependence and individual conditional expectation plots for the top 4 features, excluding the TSH level. The associations between those 4 top continuous features and prediction probability were nonlinear (Multimedia Appendix 8). In general, the risk of malignancy tended to rise as the mean tumor diameter or the coefficient of variation of TSH increased, and the risk of malignancy tended to decrease as the mean TSH score or the mean TSH increased.

‎

Moreover, we compared the characteristics of FTC groups with incorrect and correct predictions. FTC predicted as FTA by the model was classified into the incorrect-predicted group, while FTC predicted as malignant correctly by the model was then classified into the correct group. The mean diameter of the tumor was smaller in the incorrect-predicted group compared to the correct-predicted group (incorrect vs correct mean 2.88, SD cm 1.38 vs mean 3.71, SD 1.36 cm; mean diameter W=1474.5; P=.02).

A Systematic Review of the Previous Studies

After screening citations, we eventually included 7 studies in this systematic review (refer to Figure 3). The characteristics of the included studies are presented in Table 3. The sample sizes of the studies ranged from 18 to 888 patients [13-19]. The ratio of FTA to FTC in previous studies varied from 0.64 to 4.00 [13-19], which was much smaller than the ratio observed in our study (5.50) and in the real population, where the ratio of FTA to FTC can be as high as 5:1 [3]. In total, 3 studies even set the ratio close to 1 to address the imbalance [15,17,19]. As for the model selection, 4 studies developed deep learning models [13,16,17,19], 1 study used a random forest model [15], and the other 2 studies only established linear regression models [14,18], without concerning nonlinear associations or complicated interactions. Previous studies did not use gene mutations and other biomarkers as predictive variables. Except for 1 study from South Korea, which reported an AUROC of just 0.612 [13], the AUROC of the models in the other studies ranged from 0.75 to 0.96. However, none of them used AUPRC as a metric to assess discrimination. As for the interpretation of the models, Lin et al [15] assessed the feature importance, Tang et al [14] drew a nomogram, Li et al [18] developed F-TIRADS scoring criteria, and Yang et al [17] drew a heat map to visualize the importance of pixel regions, but the other 3 studies did not further explore interpretability, including feature importance, or the linear and nonlinear associations and interactions between features and targets.

‎

Additionally, we tested the F-TIRADS scoring criteria developed by Li et al [18] with our dataset. The criteria specify 6 key features for scoring. After filtering our data to include only cases with complete information for these 6 features, we selected 1025 tumors from 993 patients as an external test set. When applying the F-TIRADS scoring criteria, the predictive performance was suboptimal. With a threshold for FTC risk set at >90%, the model achieved an accuracy of 0.82, sensitivity of 0.04, specificity of 0.99, and precision of 0.53. When using a >50% FTC risk threshold, the sensitivity increased to 0.27, while accuracy, specificity, and precision decreased to 0.79, 0.91, and 0.39, respectively. For threshold-independent metrics, the AUROC and AUPRC for the F-TIRADS scoring criteria were 0.47 and 0.59, respectively, in our external test set.

Discussion

Principal Findings

This study systematically established the interpretable machine learning–based model to address the challenge of clinical decision-making for the FTN before surgery. We developed a model using 22 readily available predictors with a preferable AUROC (0.79, 95% CI 0.77‐0.81). Additionally, the model demonstrated excellent interpretability, identifying the mean TSH score, mean tumor diameter, mean TSH, coefficient of variation of TSH, and TSH level as the most important predictors. After comparing groups of incorrect-predicted and correct-predicted FTC, we found that smaller FTCs were more likely to be misclassified as FTA.

Comparison to Prior Work

It is crucial to evaluate the performance of clinical prediction models comprehensively, that is, the models should be well performed in discrimination and calibration. Concerning the discrimination, our developed model was comparable to that of previous studies aimed at predicting the malignancy risk of FTN before surgical treatment. For example, according to Li et al [18], the AUROC reached 0.76 in the LASSO regression model consisting of ultrasound features (the ratio of FTA to FTC: training set: 2.74 and validation set: 3.70); also, in the LASSO regression model, the AUROC reached 0.913 in discriminating FTA from FTC on selected clinical parameters, computed tomography signs, and radiomic features referring to Tang et al [14] (the ratio of FTA to FTC: 4.00). However, neither study used AUPRC as the evaluation metric. We also acknowledge that the model was derived from extremely imbalanced data (the ratio of FTA to FTC: 5.50), and in this context, the AUPRC metric for assessing model performance is more informative and intuitive than the AUROC [37]. For example, the prediction model might perform relatively well when measured by AUROC but may perform unsatisfactorily when measured by AUPRC, in the scenario of imbalanced data. Furthermore, to address the imbalance in our models, we applied oversampling and undersampling techniques, but both failed to improve performance. Our findings were in line with the previous research, which found that oversampling and undersampling generally did not enhance prediction models in large observational health datasets [39]. The possible reasons might be as follows: (1) oversampling and undersampling would modify the outcome proportions in the training data, leading to miscalibration, such as overestimated risks [39]; (2) the synthetic data generated by oversampling may not accurately represent the original distribution of minority class, potentially affecting classification performance [40]; and (3) undersampling reduced the number of majority class samples, limiting the model’s ability to fully use the features of the majority class during training [41].

As one previous systematic review indicated, calibration was commonly overlooked during the development and validation of clinical prediction models [20]. However, calibration metrics are also important to assess the size of the gap between the predicted risk probability and the true risk probability. For instance, grouping can be manipulated to obscure the evaluation of miscalibration in a particular range without a calibration curve and its numerical quantification [42]. In general, our model had relatively good calibration, as the calibration slope was close to 1 and the calibration-in-the-large close to 0.

Based on the results of our systematic review, most of the previous studies show relatively satisfying AUROC, but none of them reported AUPRC. Although Li et al [18] reported a handy score-risking tool for clinicians to assess the malignancy risk of FTN at the diagnosis stage, this tool seemed to not perform ideally in the practice of our data (AUROC 0.47; AUPRC 0.59; sensitivity 0.04 [threshold 90%] and 0.27 [threshold 50%]).

Limitations and Strengths

We should interpret the study findings cautiously. As with other single-center studies, the results from this study were limited in generalizability to patients and clinical settings with distinct characteristics. However, the pathological diagnosis of FTC was highly heterogeneous across different clinical settings due to its challenge in sufficient sampling and accurate diagnosis. Therefore, we advocated for the standardization of FTC diagnosis before the conduction of a multicenter study soon. Additionally, the prediction performance of models, comparable to the previous work with similar predictors, had room to further improve. The clinical utility of the screening stage was also less than ideal. Building on the experiences and lessons learned from this study, we are conducting a prospective cohort study to further optimize the model performance through collecting other costly multidimensional predictors including genomics, ultrasound images, and videos. Besides, our study was retrospective in nature, which may introduce selection bias. Furthermore, we excluded patients with nodules of UMPs, potentially limiting the model’s accuracy in identifying borderline follicular tumors. Moreover, our models did not incorporate other potential predictors, such as genetic markers (eg, BRAF, TRET, and RAS mutations), computed tomography or magnetic resonance imaging characteristics, or family history, due to constraints in data availability.

Our study had several strengths. Our models were advantageous in the large sample size for the present topic, the clinically easy-accessible and clinician-validated predictors, and the comprehensive evaluation with the metrics appropriate for the nature of the data (imbalanced data) [43]. Furthermore, the disease distribution of FTA and FTC in the study population was fully consistent with that of patients with FTN in real-world settings, that is, we did not deliberately over- or undersample patients with any type of disease in the model development, as commonly seen in previous studies [17,18]. As such, findings from our study had theoretically better fidelity and generalizability in real-world settings. In addition, we conducted a systematic review to synthesize findings from previous studies, comprehensively integrate the evidence, and identify research gaps.

Future Directions

Our study paved the way for future research in terms of predictors, models, and targets. Concerning predictors and models, further studies might consider taking advantage of the rapidly developing deep learning models and fully using high-dimension predictors such as ultrasound images and genomics. In terms of targets, it is important to standardize the pathological diagnosis of FTC across multiple centers before conducting a future multicenter study.

Our study is also important for future clinical practice. First, findings from the interpretation of our models indicate that clinicians should comprehensively consider patients’ variables such as thyroid hormones in addition to the ultrasound results. Second, in a natural distribution population with severely unbalanced data (FTA is far more than FTC), preoperative prediction of FTA and FTC by thyroid hormone and ultrasound features alone may face challenges, especially for relatively small-sized FTCs, which are easy to miss detection.

Conclusions

In clinical practice, it remained challenging to sensitively screen, precisely diagnose, and appropriately treat patients with FTN. Interpretation of our developed machine learning–based model suggests that clinicians should also pay attention to patients’ variables such as TSH along with tumor size. However, it may be hard to correctly predict FTNs preoperatively with thyroid hormone and ultrasound features alone, especially for FTCs with small sizes. The findings of our study bridged the gaps of previous work and paved the way for connecting machine learning to interpretation in the field of FTN research. We call for subsequent studies to further examine the generalizability to other contexts.