This single-center retrospective study consecutively included patients who underwent US-FNA of TNs between March 2018 and March 2023 (n = 4,210). To evaluate temporal generalizability, the dataset was divided chronologically into two cohorts: March 2018 to February 2022 (training dataset) and March 2022 to March 2023 (validation dataset). Exclusions criteria included: absence of ultrasound images, pathology-confirmed non-thyroid lesions, operator experience <3 years, multiple punctures (only the last result retained), and missing biochemical data. After exclusions, the final study population consisted of 1,573 patients in the training dataset and 672 patients in the validation dataset, with an approximate ratio of 7:3 between the two cohorts. The overall study design and patient selection flow are illustrated in Figure 1.

Clinical information and biochemical results for all patients were obtained from an electronic medical data management system. The following clinical features were recorded: patient’s age and sex. Biochemical results included free triiodothyronine (FT3), free thyroxine (FT4), thyroid-stimulating hormone (TSH), antithyroid peroxidase autoantibody (A-TPO), thyroglobulin antibody (A-TG), thyroglobulin (TG), and thyrotropin receptor antibody (TRAb). All biochemical tests were conducted within 1 month of the FNA.

All cytopathologic examinations were performed by two pathologists with >8 years of thyroid cytopathology experience and subsequently reviewed by a senior pathologist with >15 years of experience. Findings were classified according to the 2023 revision of BSRTC into four groups: Group 1 (BSRTC II), Group 2 (BSRTC I/III), Group 3 (BSRTC IV), and Group 4 (BSRTC V/VI).

Ultrasound data were retrieved from the institutional imaging system. Two US radiologists (>7 years of thyroid imaging experience) independently assessed thyroid echotexture, nodule position, capsule distance, size, volume, composition, echogenicity, echotexture, margin, shape, orientation, calcifications, posterior features, halo and Adler’s semiquantitative grading for nodule blood flow (Grades 0–3). Discrepancies were resolved by consensus with a senior radiologist (>20 years of experience).

SPSS statistical software (version 20.0; IBM Corporation, Armonk, NY, USA) was used for the statistical analysis. Baseline characteristics between the training and validation datasets were compared using the Mann–Whitney U test for continuous variables and the chi-square or Fisher’s exact test for categorical variables. Univariate analyses were further performed within the training dataset to identify factors associated with pathological classification, applying the Kruskal–Wallis test for continuous variables and the chi-square or Fisher’s exact test for categorical variables across the four groups. A p-value of <0.05 was considered statistically significant.

The prediction model (B-Model) was developed using multivariable logistic regression in SPSS based on training dataset, and it adopted a three-stage architecture as illustrated in Figure 1: (1) distinguished benign from non-benign nodules (Group 1 vs. non-Group 1) by Equation P₁; (2) differentiated malignant from non-malignant nodules (Group 4 vs. non-Group 4) by Equation P₂; (3) separated follicular neoplasm from indeterminate/atypia nodules (Group 3 vs. Group 2) by Equation P₃. Each equation had two versions: one that included biochemical indicators as independent variables P(w), and another that did not include biochemical indicators as independent variables P(w/o). For other special circumstances, a supplementary version was designed P(c). Multivariable logistic regression analyses with backward stepwise selection were applied to identify independent variables x_1-i. Based on clinical significance or published reports, we graded each risk factor, selected an appropriate grade as the baseline risk reference value, and recorded the score as 0 (1, 6, 13). β_0-i is the regression coefficient of each independent variable. Using these parameters, we calculated P as the dependent variable corresponding to each risk factor classification using the following formula, where exp denotes the natural exponential function:

The dependent variable P in the equation above uses 0.5 as a threshold value. Similar cascaded/sequential logistic regression approaches have been applied in recent medical prediction studies to improve classification performance and manage class imbalance (15–17).

The data in the validation dataset were used to select the equations and validate the performance of the prediction models. By substituting the data into previously established equations and considering the actual pathologic results as the gold standard, the sensitivity (SEN), specificity (SPE), accuracy (ACC), positive predictive rate (PPV), negative predictive rate (NPV) and area under the receiver operating characteristic curve (AUC-ROC) of each equation were evaluated. Finally, the rate of reduction in rFNAs after the B-Model implementation was calculated using the following equation:

(FN: True Group 2/3 cases incorrectly classified as Group 1/4 by B-Model).

In the training dataset, the final cohort included 1,573 patients [median age: 48 years (IQR: 38–57)] of the initial 2,818 patients, after the exclusion of 1,245 patients. In the validation dataset, the final cohort included 672 patients [median age: 50 years (IQR: 40–58)] of the initial 1,392 patients, after excluding 720 patients. The patient characteristics, US features, and biochemical results are shown in Table 1. Overall, no significant statistical differences were observed between two cohorts for most baseline characteristics except three laboratory indicators (FT4, A-TG, and A-TPO; p = 0.047, <0.001, and 0.002, respectively). These differences likely reflect case-mix shifts from time-based cohort division and variability in laboratory assays.

In the training dataset, univariate analysis identified significant differences (p < 0.05) in 2 patient characteristics, 15 US features, and 4 biochemical markers across the groups (Table 2). Specifically, thyroid echogenicity and A-TG levels were significantly different between Groups 1 and 3 (p = 0.047 and p = 0.046, respectively) whereas FT4 levels were significantly different between Groups 2 and 4 (p = 0.032). All significant variables were included as independent covariates in the subsequent multivariate analysis.

There versions of Equation P₁ were derived: P₁(w/o) (χ² = 457.323, p < 0.001), P₁(w) (χ² = 300.627, p < 0.001), and P₁(c) (χ² = 300.627, p < 0.001). P₁(c) was generated by cross-validation to address the absence of biochemical indicators in P₁(w). Two versions of Equation P₂ were developed: P₂ (w/o) (χ² = 324.479, p < 0.001) and P₂ (w) (χ² = 198.300, p < 0.001). Two versions of Equation P₃ were established: P₃ (w/o) (χ² = 148.499, p < 0.001) and P₃ (w) (χ² = 98.663, p < 0.001).

The validation results showed that among the three Equation P₁ variants, P₁(c) demonstrated the highest SEN (88.3%), SPE (68.0%), ACC (83.1%), PPV (89.2%), and NPV (66.1%), while maintaining comparable ROC-AUC (0.830 vs. 0.842/0.842 in P₁(w/o)/P₁(w), all p < 0.001). The reduced variable count (from 10 to 6) enhanced clinical utility. In the final selected Equation P₁, significant predictors included markedly hypoechoic feature (OR: 10.286, 95% CI: 6.118–17.296), hypoechoic feature (OR: 4.703, 95% CI: 3.190–6.932), irregular/extra-thyroidal extension (OR: 1.705, 95% CI: 1.180–2.463), enhanced posterior features (OR: 1.853, 95% CI: 1.265–2.715), and shadowing (OR: 2.809, 95% CI: 1.220–5.031), whereas lobulated shape showed nonsignificant association (OR: 1.122, 95% CI: 0.636–1.980). Isoechoic/hyperechoic pattern, oval-to-round shape, and absent posterior features were identified as independent protective factors for benign nodules.

Among the 498 non-Group 1 cases predicted by Equation P₁, Equation P₂ (w) demonstrated higher SEN (80.6% vs. 74.4%) and NPV (74.0% vs. 73.3%) compared to P₂ (w/o), with comparable ROC-AUC (0.735 vs. 0.759, both p = 0.000). Thus, P₂ (w) was selected to reduce missed diagnoses of malignancy. Key risk factors in Equation P₂ included isthmus location (OR: 4.000, 95% CI: 1.475–10.843), size > 5 mm (highest risk at 5–10 mm; OR: 3.058, 95% CI: 1.671–5.596), markedly hypoechoic/hypoechoic features (OR: 20.203, 95% CI: 5.203–81.179), taller-than-wide shape (OR: 5.165, 95% CI: 2.889–9.235), microcalcifications/complex calcifications (OR: 1.199, 95% CI: 0.626–2.296), and elevated TRAb (OR: 1.628, 95% CI: 1.119–2.368). These were independent predictors of malignant nodules.

Among the 181 cases predicted as neither Group 1 nor 4 by Equation P₁ and P₂, Equation P₃(w) showed higher SPE (96.0% vs. 95.4%) than P₃(w/o) with similar SEN (both 37.5%), ACC (93.4% vs. 92.7%), PPV (both 2.9%), NPV (70.0% vs. 72.7%), and ROC-AUC (0.814 vs. 0.837, both p = 0.000). The predictive performance of Equations P₁, P₂ and P₃ in the validation dataset are presented in Table 3 and Figure 2.

For the validation dataset, the number of cases correctly predicted by the B-Model were 115, 91, 3, and 195 in Groups 1, 2, 3, and 4, respectively. The prediction results of B-Model in the validation dataset are presented in Table 4. True Group 2/3 cases were 255, and true Group 2/3 cases incorrectly classified as Group 1/4 by B-Model was 153. The rFNA reduction rate was 40%.