This retrospective observational study was approved by the 960th Hospital of the PLA Joint Logistics Support Force Research Ethics Committee (No. 2024-171). Research participants were performed in accordance with the Declaration of Helsinki, and each patient provided written informed consent. All eligible patients were informed about the use of their data for study and had the option to decline to participate (22).

This study included patients with thyroid nodules diagnosed as Bethesda III nodules based on initial FNA between January 2020 and June 2024. Final diagnoses were confirmed through repeat biopsies, 12-gene testing, or surgery (23). Inclusion criteria included: (1) thyroid nodules classified as Bethesda III on FNA cytology; (2) preoperative ultrasound, molecular testing, and thyroid function tests; (3) patients underwent thyroidectomy; (4) at least one repeat FNA on the same Bethesda III nodule within one year; (5) 12-gene panel testing. Exclusion criteria included: (1) other thyroid cancer types, such as follicular carcinoma or medullary thyroid carcinoma; (2) patients lacking laboratory, imaging, or pathological data; (3) uncertain repeat FNA results; (4) any mutations detected in the 12-gene test; (5) patients without at least six months of follow-up after repeat FNA or 12-gene testing. Patients were randomly assigned to a training set and a validation set in a 7:3 ratio. The training set was used to develop models with various machine learning algorithms, and the validation set was used to evaluate model performance.

Data for the variables assessed in this study were collected from patients’ hospitalization electronic medical records (EMRs), including basic patient information, thyroid imaging reporting and data system, cytological assessments, molecular testing following US-FNA, laboratory indicators (within one month before surgery) and postoperative pathological data. Basic patient information included age and sex. Ultrasound features of thyroid nodules were analyzed using TI-RADS terminology, including size, composition, echogenicity, margin, shape, echogenic foci, halo, orientation, location, color Doppler flow imaging (CDFI) pattern, echotexture, posterior features, cervical lymph nodes and solitary nodule (24). The largest transverse, anteroposterior, and vertical diameters of all nodules were recorded, with the largest of these three measurements used to assess nodule size. Additional features of the nodule were recorded according to C-TIRADS. C-TIRADS classified the nodules by assigning points for composition, echogenicity, shape, margin, and echogenic foci to determine the TIRADS level. Cervical lymph nodes were considered abnormal if they exhibited any of the following features: (a) loss of central hilar echo, (b) cystic change, (c) calcification, (d) cortical hyperechogenicity, (e) increased and irregular vascularity, or (f) cervical lymphadenopathy. Two radiologists reviewed and documented all sonographic characteristics of the thyroid nodules. Disagreements were resolved through discussion or consultation with a third radiologist. In cytological assessments and molecular testing, cytologists examined the cells and cellular structures according to the second edition of the Bethesda System for Reporting Thyroid Cytopathology (3). Bethesda III cases were reviewed and categorized into the following subtypes: (1) AUS-nuclear atypia: Includes focal nuclear atypia or mild but extensive nuclear atypia; (2) AUS-other: Includes architectural atypia (often sparsely cellular samples predominantly composed of microfollicles) and Hurthle cell atypia (oncocytic features). Two molecular tests were available at our institution: BRAF and a twelve-gene molecular panel. Patients underwent twelve-gene panel testing, which included BRAF, RAS (NRAS, HRAS, KRAS), PIK3CA, AKT1, RET, CCDC6-RET, NCOA4-RET, TP53, ETV6-NTRK3, TPM3-NTRK1, TERT, and PAX8-PPARG promoters. Laboratory indicators obtained included thyroid stimulating hormone (TSH), free thyroxine (FT4), free triiodothyronine (FT3), thyroglobulin antibody (TgAb), thyroid peroxidase antibody (TPOAb), and thyroglobulin (Tg) levels. All predictive factors were derived from objective data and image archiving in the EMRs.

The definitive diagnosis of malignant tumors was confirmed through histopathological analysis of surgically excised tissue or repeated fine-needle aspiration biopsies with malignant results (Bethesda categories V or VI). In contrast, benign nodules were diagnosed based on one of the following criteria: (1) histopathological confirmation from surgically excised specimens, (2) absence of mutations in a 12-gene panel, or (3) at least one benign FNA result (Bethesda II) from repeated biopsies, with no subsequent malignant findings. All patients who did not undergo surgical intervention were followed up with ultrasonography for at least six months to assess any suspicious changes in the nodules. Histopathological results were independently reviewed by two experienced pathologists, and any discrepancies were resolved through consultation with a third pathologist.

We calculated the required sample size using the events per variable (EPV) metric, a widely recognized method in statistical analysis (25–27). We also followed the 4-step procedure proposed by Riley et al. (28) to calculate the required sample size (29). According to guidelines, the malignancy rate of Bethesda III nodules typically ranges from 13% to 30% based on follow-up of surgically resected nodules (3). Since we intended to include patients who had not undergone surgery, the expected proportion of the endpoint event was estimated at 0.15. The calculation process, formulas, and results are shown in Supplementary Table S2 .

Statistical analyses were conducted using SPSS 25.0 software and R (version 4.4.2; R Foundation, Vienna, Austria), with a significance level set at p < 0.05. The “mice” package in R was used to assess the missing data mechanism, and the VIM package was used for data visualization. Multiple imputation was used to handle missing data, with 10 imputations performed based on the established MAR mechanism, as recommended by standard guidelines. The Multiple Imputation by Chained Equations (MICE) method was used for imputation, implemented through the “mice” package in R. The imputation model with the lowest Bayesian Information Criterion (BIC) was chosen for data optimization. Continuous variables were imputed using predictive mean matching, categorical variables using logistic regression, and multinomial variables using multinomial logistic regression. All relevant covariates, including predictors, outcome variables, and other variables not included in the predictive model, were incorporated into the imputation model to capture the relationships among the variables.

Supplementary Figure S1 presents the complete study flowchart. Before constructing the predictive model, six feature selection methods were applied to mitigate the high correlation between predictor variables and improve both the predictive accuracy and interpretability of the model. These methods were chosen based on their proven effectiveness in handling high-dimensional data and identifying the most relevant predictor variables in predictive modeling. In clinical prediction models, the combined use of multiple variable selection methods can compensate for the limitations of individual approaches, thereby enhancing the stability of variable selection and the predictive performance of the model (30, 31). The feature selection methods used include:

We summarized 29 candidate predictors selected through various filtering methods and took their intersection. Finally, we consulted with clinical experts and combined them with clinical reality to determine the final predictors for constructing Bethesda III nodules malignancy prediction model.

Six machine learning models were used to predict the malignancy of Bethesda III nodules: LR (glm stats package), DT (rpart package), RF (randomForest package), XGB (xgboost package), SVM (e1071 package), and KNN (kknn package). The selected clinical features were fed into the six models using various machine learning algorithms. To evaluate the models’ performance on unseen data, five-fold cross-validation was applied to obtain the parameters. Performance metrics from the confusion matrix were used to assess model efficacy in both the training and validation cohorts, including the Receiver Operating Characteristic curve, sensitivity (SEN), specificity (SPE), positive predictive value (PPV), negative predictive value (NPV), accuracy (ACC), F1 score, and Brier score. A calibration curve was used to compare predicted probabilities with actual outcomes. Decision Curve Analysis was used to evaluate the net benefit of the models at various thresholds. Additionally, the DeLong test was used to determine whether significant differences existed in the AUC values of the models. Based on the evaluation of these metrics in both the training and testing sets, the optimal model was selected.

LR was selected as the optimal model based on its performance. The model’s performance was demonstrated using 500 bootstrap samples, including the ROC curve, calibration curve, and DCA results. To assess model interpretability and evaluation, we calculated the AUC to evaluate the model’s discriminative ability and examined its calibration using the Hosmer-Lemeshow (HL) test. DCA was used to evaluate the net benefit at different thresholds, providing a comprehensive assessment of the model’s effectiveness in real-world medical decision-making scenarios. Bootstrapping, a statistical method, was applied to estimate model accuracy by repeatedly sampling with replacement from the original dataset. The model was trained on these new datasets and evaluated using out-of-bag data. This process was repeated several times to obtain a distribution of performance metrics, yielding a robust estimate of model reliability and variance. The nomogram score is the sum of the scores assigned to each risk factor, where higher scores indicate a greater risk of Bethesda III nodules malignancy. This graphical tool simplifies the estimation of individual risk or probability based on the variables predicted by the model.

All analyses were conducted in R. For continuous variables, the Shapiro-Wilk test was used to assess their normality. This test was selected because it generally performs better than other tests, such as the Kolmogorov-Smirnov test, particularly with small sample sizes. Variables following a normal distribution were described using the mean ± standard deviation, while non-normally distributed variables were described using the median and interquartile range, providing more robust measures of central tendency and variability in the presence of outliers. Frequency and percentage were reported for categorical variables. To assess collinearity between variables, the variance inflation factor (VIF) was calculated. A VIF value below 5 and tolerance greater than 0.1 typically indicate no significant collinearity between variables. This metric helps ensure that our regression models are not unduly affected by multicollinearity, which can distort the estimated relationships between predictors and outcomes. Fisher’s exact test was used for categorical variables with low expected frequencies to ensure accurate significance testing. This test provides a precise method for determining the likelihood of observing a given set of frequencies in categorical variables, particularly useful with small sample sizes. Several R packages were used for specific analyses: comparegroups for baseline description, glmnet for LASSO regression, forestmodel for forest plots, pROC, ggROC, and fbroc for discriminative analysis, PRROC for PR curves, rms for calibration using val.prob and calibrate functions, ResourceSelection for the Hosmer-Lemeshow test, dcurves for DCA, and rms for nomogram construction.

Between January 2020 and June 2024, 7,371 patients underwent ultrasound-guided fine needle aspiration, while 2,761 patients underwent thyroidectomy. Among the thyroid nodules assessed by US-FNA, 10.2% (754/7,398) were diagnosed with Bethesda III nodules. Of these Bethesda III nodules, 31.6% (238/754) underwent surgical resection, 7.3% (55/754) underwent repeat biopsies, and 6.9% (52/754) were subjected to 12-gene panel testing. Notably, 36.4% (20/55) of the Bethesda III nodules that underwent repeat biopsies had indeterminate results. A total of 69 patients were excluded due to incomplete records (16), other pathologically confirmed thyroid cancer types (2), uncertain repeat biopsy outcomes (20), mutations found by 12-gene panel testing (1), or insufficient follow-up for at least six months after repeat biopsy or 12-gene panel testing (30). Consequently, the current study included 271 patients: 213 females and 58 males, with a total of 276 Bethesda III nodules. Approximately 65.2% (180/276) of the nodules were malignant, while 34.8% (96/276) were benign. BRAF mutations were detected in 41% (113/276) of the cases. In our dataset, only 36 cases (13.04%) were classified as AUS and 25 (9.06%) as FLUS. Neither AUS nor FLUS showed statistical significance in the univariate analysis. The median age of patients was 51 years (IQR 41.00–57.25), and the median nodule size was 5.2 mm (IQR 4.00–7.70). In the current dataset, five variables were missing, all related to thyroid function tests, resulting in 232 missing values, or approximately 2.8% of the total data points. These missing values were distributed among multiple variables, as detailed in Supplementary Table S1 and Supplementary Figure S2 . The cohort was randomly divided into a training set (n = 193) and a validation set (n = 83). As shown in Table 1 , no significant differences were observed in the baseline characteristics between the two groups (P > 0.05).

The predictor variables selected by the six methods are presented in Table 2 . Detailed parameters for all methods (SR-FS, SR-BS, SR-BE, LASSO, Boruta and RF-RFE) can be found in Supplementary Table S4 and Supplementary Figures S3 - S5 . A total of 10, 6, 6, 5, 8, and 10 predictors were identified using the six selection methods. As shown in Supplementary Table S4 , variables such as TSH, FT4, and Echotexture did not show a significant association with malignancy in the initial univariate screening and were therefore not carried forward into the final model building stage. Figure 1 illustrates the intersection of predictors selected by the six methods. Predictors with more than six intersections were selected as final predictors, including BRAF, composition, shape, orientation, and TI-RADS. Ultimately, five predictors were included in the development of the model.

In our comparative analysis, we evaluated five machine learning models against the LR model. Although complex models such as RF achieved better training set performance, the LR model exhibited the most robust and clinically meaningful results in the validation set, as evidenced by Table 3 and Figure 2 . The six models were developed using six distinct machine learning algorithms. The estimated odds ratios for the logistic regression model are presented in Supplementary Table S5 and visualized in a forest plot ( Supplementary Figure S6 ). We further illustrate visualizations for the other models, including the relative importance of potential features and heatmaps of confusion matrices for the decision tree, random forest, extreme gradient boosting, support vector machine, and k-nearest neighbors models derived from the training cohort. However, due to the nature of the KNN model, the ranking of feature importance is not applicable in this case.

The RF model demonstrated the highest area under the ROC curve (AUC, 0.923; 95% CI: 0.887–0.959) in the training set ( Figure 2A ), whereas the LR model achieved the highest AUC (0.823; 95% CI: 0.732–0.915) in the validation set ( Figure 2B ). The AUC values for the remaining models in the validation set were as follows: (1) DT: 0.758, 95% CI: 0.657–0.859; (2) RF: 0.792, 95% CI: 0.694–0.890; (3) XGB: 0.817, 95% CI: 0.722–0.912; (4) SVM: 0.822, 95% CI: 0.730–0.913; (5) KNN: 0.771, 95% CI: 0.668–0.873 ( Figure 2B ). Based on the DeLong test ( Supplementary Tables S7 , S8 ), no significant difference in AUC was observed in the validation cohort. This may be due to the limited sample size or the similar performance of the models.

The DT model showed the best consistency between observed and predicted results in both the training and validation sets ( Figure 2C ). In contrast, the LR, SVM, KNN, and XGB models exhibited similarly good consistency between observed and predicted results in the validation set ( Figure 2D ). The consistency between observed and predicted results in the RF model was less stable in the validation set ( Figure 2D ).

In the training set, all models demonstrated similar DCA results ( Figure 2E ), whereas the LR, SVM, and XGB models achieved the best DCA outcomes in the validation set ( Figure 2F ). In the validation set, using alternative models resulted in a greater net benefit compared to no treatment or full treatment strategies when the threshold probability was <80%.

The RF model exhibited the highest precision and specificity in the training set, whereas the LR model achieved the highest recall and F1 score, as well as the second lowest Brier score. In the validation set, the LR model achieved the highest recall and F1 score, whereas the SVM model demonstrated the highest precision. Detailed information can be found in Table 3 . Model performance in the validation set served as the primary criterion for identifying the optimal model. Although the LR model exhibited moderate specificity (0.55), its superior AUC (0.823), recall (sensitivity, 0.85), and F1 score (0.81) justified its selection as the optimal model. For malignancy risk assessment, minimizing false negatives through high recall remains clinically critical.

In summary, while the RF model may exhibit potential overfitting, the LR model demonstrates not only strong interpretability but also achieves the highest AUC in the validation set. Taking into account model performance, complexity, generalization capability, and practicality, the LR model was ultimately chosen. The forest plot illustrates variables with P-values greater than 0.05 (composition), and due to its clinical significance in daily medical practice, this variable was retained in the logistic regression model for interpretation of its effect size.

The nomogram was constructed by incorporating five variables—BRAF, composition, shape, orientation, and TI-RADS—into the predictive model ( Figure 3 ). Following 500 bootstrap iterations, the LR model demonstrated an AUC of 0.871 (95% CI: 0.817–0.926) in the training set and 0.846 (95% CI: 0.762–0.930) in the validation set ( Figures 4A, B ). After 500 bootstrap iterations for calibration, the calibration curve closely aligned with the ideal diagonal. The average absolute errors (MAEs) for the two datasets were 0.013 and 0.031, respectively, indicating that the predicted probabilities closely aligned with the actual outcomes across different samples ( Figures 4C, D ). The Hosmer-Lemeshow test confirmed good consistency between the predicted and observed results (P = 0.414). DCA showed that applying the LR nomogram provided greater net benefit than no treatment or full treatment strategies when the threshold probability ranged from 16% to 96% in the training set and from 15% to 95% in the validation set ( Figures 4E, F ). Using this model, clinicians can more accurately assess the risk of malignancy in Bethesda III nodules, providing optimized management and treatment options. The integration of this predictive model significantly enhances the precision of patient management. This highlights the importance of leveraging precise, data-driven decisions in clinical practice. We further assessed the independence of the variables in the model by calculating the VIF and found that all VIF values were well below the common threshold of 5, indicating low multicollinearity and confirming the model’s stability and reliability. The detailed VIF values are provided in Supplementary Table S6 . Our analysis compared the ROC curves of the nomogram with those of the five individual predictors for malignancy risk assessment. The nomogram exhibited superior diagnostic performance, as illustrated in Supplementary Figure S14 .

Compared to the complex logistic regression formula, the nomogram is simpler, more intuitive, and clinically practical. To use the nomogram, a line is drawn from the value of each feature to the “points” axis, obtaining the corresponding score. After summing these points, the total score is located on the “total points” axis. Finally, a line is drawn downward from the total score to the “probability of malignancy” axis to determine the corresponding risk. For example, when a patient is BRAF-positive with an ultrasound indicating a solid nodule, a regular shape, a vertical-to-horizontal ratio <1, and a TI-RADS score of 4a, their “BRAF” score is 91, “composition” is 32, “shape” is 0, “orientation” is 0, and “TI-RADS” is 24. The total score is 147, corresponding to a malignancy probability of 0.72 (72%). Therefore, the surgeon may consider this patient at high risk for Bethesda III nodules malignancy and recommend diagnostic lobectomy ( Figure 3A ).

To support clinical implementation of this predictive model, we developed a web-based calculator that enables physicians to input patient-specific clinical data and immediately obtain malignancy risk estimates for Bethesda III nodules, thereby supporting more informed clinical decisions. The tool streamlines risk assessment while facilitating personalized patient management through individualized risk stratification. The calculator is publicly accessible online at https://13583155338-l.shinyapps.io/PredictionofAUS/ ( Figure 3B ).