Descriptive statistics

Table 1 presents the descriptive analysis of the 103 cadavers in this study, comprising 62 males (mean age 47.3 ± 16.8 years) and 41 females (mean age 40.1 ± 18.8 years). The overall age range was 18–86 years, with an average age of 44.5 ± 17.9 years and a median age of 41 years. In terms of BMI, 47.6% were of normal weight (18.5–24.9 kg/m2, consisting of 26 males and 23 females), 17.5% were underweight (< 18.5 kg/m2, consisting of 8 males and 10 females), and 34.9% were overweight (> 25 kg/m2, consisting of 28 males and 8 females) (Supplementary File 2).

Table 1 Descriptive statisticsInter-observer analysis results

TEM analysis was performed to assess the inter-observer reliability across 36 dependent and 21 independent variables. Approximately 90% of the variables exhibited TEM values below or around 2 to 2.5 mm, indicating high measurement consistency [26,27,28]. Only two variables, FH-NP (4.45) and 31-Occlusal Plane (3.09), prominently exceeded the threshold (Fig. 3).

Fig. 3

TEM results of the measured variables

Correlation analysis between independent and dependent variables

Figure 4 presents the results of the Pearson correlation analysis between the independent and dependent variables. The analysis revealed that 11 independent variables were positively correlated with FSTT across nearly all landmarks: age, BMI, FMIA, IMPA, SNA, SNB, S-N-Gn, 21-NS, 31-GoGn, and 21-NP. Conversely, 10 variables were negatively correlated with the FSTT: sex, CI, SC, FMA, SN-GoGn, FH-NP, S-Gn-FH, 21–31, 31-occlusal plane, and AB-occlusal plane.

Fig. 4

Visualization of correlation analysis between the independent variables and FSTT using a heat map. Bold numbers indicate that the correlation is significant (p < 0.05)

BMI exhibited the strongest correlation with FSTT compared to the other independent variables. Among the 36 landmarks evaluated in this study, 33 showed a strong relationship with BMI. The next variable that correlated closely with FSTT across nearly all landmarks was sex, with males displaying a higher mean FSTT than females, especially at midline landmarks, the upper third of the face, and the subnasal area.

For the orthodontic profiles, SNB, SN-GoGn, and S-N-Gn emerged as the three variables most strongly correlated with FSTT. All three variables showed significant correlations primarily in the upper third of the face, midface, and cheek regions.

Correlation analysis among independent variables

Figure 5 shows the correlations among the independent variables. The orthodontic profile variables, particularly SC, Tweed, and Northwestern, exhibited strong intercorrelations, either positive or negative correlations. The highest positive correlation was observed between SNB and S-N-Gn (r = 0.94), while the highest negative correlation occurred between SC and NAP (r = -0.94). Additionally, FMIA and 21-NS had the highest number of significant correlations (p < 0.05) with the other variables, with a total of 14 significant intercorrelations.

Fig. 5

Visualization of correlation analysis among the independent variables using a heat map. Bold numbers indicate that the correlation is significant (p < 0.05)

Multicollinearity analysis and principal component analysis

Considering the high correlations among independent variables, we conducted VIF analyses to identify the effect of multicollinearity within each variable. The results showed high VIF values for the orthodontic profile variables, particularly for SC, Tweed, and Northwestern (Supplementary File 1).

To address high multicollinearity, we performed PCA on those variables with high VIF values. Through 10 iterations of dimension reduction, eliminating variables with a measure of sampling adequacy (MSA) around 0.5 in each anti-image matrix, we obtained two PCs (PC1 and PC2) with a Kaiser-Meyer-Olkin (KMO) value of 0.79 and a total eigenvalue/variance of 76.6%. The KMO value indicates strong intervariable correlations, supporting the reliability of PCA and confirming that multicollinearity is no longer a significant issue. Meanwhile, the eigenvalue was ideal, while further variable reduction was avoided for two reasons: no MSA values in the anti-image table were below 0.5, and further reduction decreased the total variance to below 70%, compromising data representativeness [29, 31, 32].

These two PCs included the variables SC, FMIA, IMPA, 21–31, 31-GoGn, 31-Occlusal Plane, and 21-NP (Supplementary File 1). The first PC was influenced by Tweed and Northwestern variables, while the second was driven by SC and Northwestern variables. PC1 was predominantly influenced by broader variations, such as FMIA, IMPA, 21–31, 31-GoGn, and the 31-occlusal plane, while PC2 captured narrower variables driven by SC, FMIA, and 21-NP. This suggests that PC1 reflects patterns specific to the Tweed and Northwestern profiles, whereas PC2 captures patterns specific to the SC and Northwestern profiles.

Multiple linear regression analysis for PCA model

PC1 and PC2 were used as new independent variables in the multiple linear regression analysis. Therefore, the predictive formula we developed consisted of baseline variables (age, sex, and BMI), CI, and the two PCs. We also conducted RMSE, MAE, and R2 metric analyses to evaluate the performance of the model (Table 2). Lower values of RMSE and MAE indicate better model performance, whereas higher R2 values reflect a greater proportion of variance explained by the model. From the metric analysis results, the Rhi landmark showed the lowest RMSE and MAE values, at 0.83 and 0.59, respectively. Meanwhile, the highest R² value was achieved at the Mm-L landmark, with a value of 0.61.

Table 2 Regression analysis of PCA models (age, sex, BMI, CI, PC1, and PC2)

The PC values (PC1 and PC2) were calculated by summing the product of the loadings for each independent variable with the actual value of that variable for a given subject. This calculation yielded the weighted contribution of each independent variable to the respective PC. Therefore, the formula for calculating PCs is expressed as follows:

PC1 =– (SC × 0.04) + (FMIA × 0.72)– (IMPA × 0.84) + (21–31 × 0.80)– (31-GoGn × 0.88) + (31-Occl. plane × 0.86)– (21-NP × 0.38) (1).

PC2 = (SC x 0.89)– (FMIA × 0.55) + (IMPA × 0.14)– (21–31 × 0.20) + (31-GoGn × 0.15)– (31-Occl. plane × 0.25) + (21-NP × 0.77) (2).

After obtaining the PC1 and PC2 values for each landmark, these values were incorporated into the regression formula, which is structured as follows:

FSTT = B0 + (age × B1) + (sex × B2) + (BMI × B3) + (CI × B4) + (PC1 × B5) + (PC2 × B6) (3).

Multiple linear regression analysis for baseline model

We also developed the baseline model using three variables (age, sex, and BMI) for each landmark. These models also used RMSE, MAE, and R2 as evaluation metrics (Supplementary File 1). The results were similar to those of the PCA models, with the lowest RMSE and MAE values observed at the Rhi landmark (0.86 and 0.64, respectively), whereas the highest R2 value was recorded at the Mm-L landmark (0.61). BMI was correlated with all landmarks, sex with 10 landmarks, and age with 7 landmarks. The regression equation structure was similar to that of Eq. (3) but without CI, PC1, and PC2.

Comparison of FSTT prediction accuracy: regression models vs. holdout samples BMI-based mean

A comparative analysis was conducted between the FSTT predictions derived from the primary dataset’s PCA-based regression model and those estimated using mean BMI values from the holdout dataset (Supplementary File 2). The RMSE and MAE values consistently showed a modest improvement in the PCA-based regression model predictions across the majority of landmarks. The most notable improvement in RMSE and MAE was observed at the Point-A landmark, showing a reduction of 2.55 mm and 1.56 mm, respectively.

However, middle third landmark (Apc-L/R), some lower third landmarks (Id, Point-B, Pog, and Gn), and four cheek landmarks (Spm2-L, Spm2-R, Sm2-L, and Sm2-R), demonstrated higher values of either RMSE, MAE, or both standard errors in the PCA-based regression model, indicating reduced predictive accuracy for these specific areas using the regression approach.

However, overall, these results suggest that PCA-based regression models are slightly more accurate in FSTT estimates compared to predictions based solely on BMI-based mean values.

Comparison of FSTT prediction accuracy: between the two linear regression models

The comparison of RMSE, MAE, and R² showed marginal differences between the PCA-based and baseline regression models (Supplementary File 2). In the RMSE analysis, the PCA model had slightly smaller values across almost all landmarks except for Ft-L, Zm-R, and Mm-L, which had equal values, and Pog, Fe-R, Zy-L, and Zy-R, which had larger values. For MAE, the PCA model also generally showed smaller values, except for Ft-L (equal values) and Pog, Fe-R, Zm-R, Zy-L, Kdl-R, Mm-L, Go-L, and Smp2-L (larger values). In the R² analysis, the PCA model showed larger values for half of the landmarks. Meanwhile, others, such as Fe-L, Spo-R, Mm-L, and Go-L, had equal values, and Sg, G, Pog, Fe-R, Ft-L, Ft-R, So-L, So-R, Zm-R, Zy-L, Zy-R, Kdl-R, Mm-R, and Go-R had smaller values.

After adding the orthodontic profile variables to the regression formula, landmarks in the cheek area (Spm2 and Sm2 on both sides) showed the most significant improvements. In the RMSE analysis, Spm2-L, Spm2-R, Sm2-L, and Sm2-R exhibited reductions of 1.00, 0.22, 0.24, and 0.14, respectively. In the MAE analysis, the values that decreased for these landmarks were 0.72, 0.21, 0.29, and 0.06, respectively. For R², the observed increases were 0.21, 0.06, 0.07, and 0.03, respectively. We also observed a large metric difference in the upper lip region, especially at the Point-A and Pr landmarks. For Point-A, the RMSE change was 0.17, MAE was 0.14, and R2 was 0.12. For Pr, the metric changes were 0.19, 0.17, and 0.11 for RMSE, MAE, and R2, respectively. These values exceeded the average changes for each metric for all landmarks, where the change in mean RMSE was 0.09, MAE was 0.07, and R² was 0.02.