A total of 300 Korean adults (150 female) from a large online panel using quotas for age, gender, and region of residence participated in this study. We followed the suggestions of Bentler and Chou [21] to determine the appropriate sample size. They suggested a 5:1 ratio of sample size to parameters. This proposal required a minimum of 270 participants. In this study, 300 participants were selected to compensate for missing data. To prevent biased sampling, the numbers of participants were equalized by age, gender, and region. Equal numbers of participants were aged 20–29 years, 30–39 years, 40–49 years, and 50–59 years. In addition, the number of participants by region was proportional to the population, and the gender ratio of the participants was 1:1. The average age of the 300 study participants was 39.28 years (standard deviation (SD) = 10.91), and the age range was 20–59 years. Data were collected using Embrain (https://www.embrain.com/kor/), an online survey company. We collected information about the participants’ gender, age, and living region, but did not collect any identifiable personal information. Participants were provided with information about the survey’s purpose prior to the survey, agreed to confidentiality, and received points that could be used at the survey site in return for their participation.

Data analyses were conducted using SPSS version 27 (IBM, Armonk, NY, USA) and Mplus version 8.8 (Muthén & Muthén, Los Angeles, CA, USA), with a p-value of less than 0.05 deemed statistically significant. Descriptive statistical analysis was performed using SPSS version 27, and factor analysis, reliability assessment, and measurement invariance tests were conducted using Mplus version 8.8. The factor structures of the SIAS-6 and SPS-6 were examined using CFA and ESEM, with ESEM performed using oblique target rotation. Target rotation guides the rotation of factor loadings based on a predefined target matrix, which was specified by the researcher to indicate the expected pattern of loadings. In this matrix, some elements are set to zero or near zero, reflecting the expectation that certain factors will have no or minimal loadings on specific items [19, 20].

Table 1 presents the cut-off criteria for the goodness-of-fit index used in this study [24]. The models were compared according to the fit indices listed in Table 1:

(1) Model 1: A one-factor model in which all 12 items load on the single factor of global social anxiety.

(2) Model 2: A two-factor CFA model in which six items load on the social interaction anxiety factor, while the remaining six items load on the social performance anxiety factor.

(3) Model 3: A bifactor CFA model in which all 12 items load on the general factor (global social anxiety factor) and, at the same time, on one of the specific factors (either the social interaction anxiety factor or the social performance anxiety factor).

(4) Model 4: A two-factor ESEM model in which all 12 items load on the social interaction anxiety factor and the social performance anxiety factor.

(5) Model 5: A bifactor ESEM model in which all 12 items load on the general factor of global social anxiety while simultaneously loading on the specific factors of social interaction anxiety and social performance anxiety.

In addition, the Akaike Information Criterion (AIC) was calculated to compare the fitness of several models and determine the most appropriate model. The AIC addresses the risk of overfitting and underfitting by trading off simplicity and fit of the model. If the difference in the AIC between the two models is six or more, the model with the smaller AIC is considered to have a better fit than the model with the larger AIC [25].

In the comparison of the bifactor models, the average and lowest loading values for the designated factor were calculated along with the fit indices. The average loading values of the designated factor must be greater than 0.50 and the lowest loading values of the designated factor must be greater than 0.30 for a well-defined factor [26]. Coefficient omega hierarchical ($\omega$_H) was calculated as an estimate of the general factor, while coefficient omega hierarchical subscale ($\omega$_H⁢S) was calculated as an estimate of the specific factors to assess the internal consistency coefficients of the models that exhibited the best fit. This is because the omega coefficient is more suitable than the traditional Cronbach’s alpha for measuring the internal consistency of multidimensional constructs [23]. $\omega$_H denotes the variance attributable to the general factor after controlling for all specific factors and $\omega$_H⁢S is the variance attributable to the specific factor after controlling for the general factor. The acceptable cut-off score of $\omega$_H is 0.50 [27]. The substantial cut-off score of $\omega$_H⁢S is above 0.30, the moderate cut-off score is above 0.20 and below 0.29, and the low cut-off score is below 0.19 [28]. Additionally, the explained common variance (ECV) was calculated for the bifactor models. The ECV associated with the general factor represents the share of total common variance explained by that factor, while the ECV linked to the specific factor shows the variance uniquely tied to that specific factor. This distinction allows for a clearer understanding of how much of the variance in the data can be explained by overarching constructs (the general factor) compared with how much is uniquely explained by individual factors (the specific factor). If a general factor explains more than 70% of the common variance, the measure is considered to have a unidimensional structure, rather than multiple distinct factors [29].

Measurement invariance was evaluated using Chen’s criteria [30]. Chen [30] suggests a criterion of –0.01 for changes in comparative fit index (CFI) and a criterion of 0.015 for changes in root mean square error of approximation (RMSEA). The fit indices of each model were compared. The metric invariance model, where loadings were constrained to be equal across genders, was compared with the configural invariance model, where structural parameters were constrained to be equal across genders. A scalar invariance model, where the intercept and measurement residuals were constrained to be invariant across genders, was compared with a configural invariance model.

Two-sided multivariate tests were conducted to assess the fit of skewness and kurtosis in verifying multivariate normality for CFA and ESEM. These tests examine the properties of two distribution curves: multivariate skewness and multivariate kurtosis. Results from these tests that do not support multivariate normality may call into question the validity of subsequent CFA and ESEM analyses.

Drawing on theory and previous studies, we compared five competing models: the one-factor model, two-factor CFA model, bifactor CFA model, two-factor ESEM model, and bifactor ESEM model. Prior to this, we conducted two-sided multivariate tests for skewness and kurtosis. The tests showed that the data deviated from multivariate normality (two-sided multivariate test for skewness fit = 31.97, p 0.001; two-sided multivariate test for kurtosis fit = 244.02, p 0.001). Consequently, a robust maximum likelihood estimator was employed to address non-normality. A robust maximum likelihood estimator can deliver reliable estimates even when the data fail to satisfy the usual assumptions of normality. Table 2 reports the model fit indices of the five different models. The 90% confidence interval (CI) for the RMSEA represents a range of plausible values for the population RMSEA. For the one-factor model and the two-factor CFA model, $\chi$²/df, CFI, TLI, and RMSEA showed a poor model fit and the SRMR displayed an acceptable model fit to the data. For the bifactor CFA model, $\chi$²/df, TLI, and RMSEA exhibited an acceptable fit and the CFI and SRMR indices indicated a good fit. For the two-factor ESEM model, TLI indices showed a poor model fit to the data; $\chi$²/df, CFI, and RMSEA displayed an acceptable model fit; and the SRMR indicated a good model fit. For the bifactor ESEM model, both TLI and RMSEA values suggested an acceptable model fit and the CFI, SRMR, and $\chi$²/df indicated a good fit to the data.

The relative model fit, which refers to determining which of two or more models best represents the data, was assessed using the AIC. Model 2 was a better fit than Model 1 ($\mathrm{\Delta }$AIC = –234.70). Model 4 was a better fit than Model 2 ($\mathrm{\Delta }$AIC = –63.60). Model 3 was a better fit than Model 4 ($\mathrm{\Delta }$AIC = –43.78). Additionally, Model 5 was superior to Model 4 ($\mathrm{\Delta }$AIC = –41.00). Thus, the bifactor CFA and ESEM models were the best alternatives.

The standardized loadings from the bifactor CFA models of the SIAS-6 and SPS-6 are listed in Table 3. In the bifactor CFA model, the average loading values on the general factor were above 0.50 ($\lambda$ = 0.576–0.848, M = 0.662). However, the average loading values were below 0.50 and the lowest loading values were below 0.30 on social interaction anxiety ($\lambda$ = 0.121–0.651, M = 0.426) and social performance anxiety ($\lambda$ = 0.018–0.711, M = 0.399) factors. This suggests that each item had stronger loadings on the general factor compared with their specific factors. For the social interaction anxiety factor, five social interaction anxiety items (items 1, 2, 3, 4, and 5) loaded significantly onto the social interaction anxiety factor ($\lambda$ $\ge$0.30), whereas only one social interaction anxiety item (item 6) failed to load onto the specific factor ($\lambda$ = 0.121). Five social performance anxiety items (items 7, 9, 10, 11, and 12) loaded significantly onto the social performance anxiety factor ($\lambda$ $\ge$0.30), whereas only one social performance anxiety item (item 8) failed to load onto the specific factor ($\lambda$ = –0.018). That is, out of the 12 items, 10 items loaded significantly onto their respective specific factors, suggesting that specific factor cross-loadings remained in the bifactor CFA model even after controlling for the general factor.

Table 4 shows the standardized loadings from the bifactor ESEM models of the SIAS-6 and SPS-6. All items showed significant loadings for the general factors. In the bifactor ESEM, the average loadings on the general factor were above 0.50 ($\lambda$ = 0.536–0.905, M = 0.659). However, the average loadings were below 0.50 and the lowest loadings were below 0.30 for the social interaction anxiety ($\lambda$ = 0.157–0.644, M = 0.444) and social performance anxiety ($\lambda$ = 0.093–0.677, M = 0.404) factors, showing that each item was loaded better on the general factor than on their respective specific factors. For the social interaction anxiety factor, five social interaction anxiety items (items 1, 2, 3, 4, and 5) showed a significant loading on the social interaction anxiety factor ($\lambda$ $\ge$0.30), whereas one item (item 6) failed to load onto the factor ($\lambda$ = 0.157). Five social performance anxiety items (items 7, 9, 10, 11, and 12) had a significant loading on the social performance anxiety factor ($\lambda$ $\ge$0.30), whereas one social performance anxiety item (item 8) did not load onto the target factor ($\lambda$ = –0.093). In other words, 10 of the 12 items exhibited significant factor loadings on their corresponding specific factors, suggesting that specific factor cross-loadings persisted in the bifactor ESEM model, even after accounting for the general factor.

In the bifactor CFA and ESEM models, items related to social interaction anxiety and social performance anxiety loaded significantly onto both the general and specific factors, suggesting that their variances were divided into general and specific components. Consequently, the social interaction and performance anxiety factors were distinct from each other.

For the bifactor models, we computed various indices to assess reliability, including $\omega$_H, $\omega$_H⁢S, and the ECV. As shown in Table 3, $\omega$_H for the general factor (12 items) was 0.795, while $\omega$_H⁢S was 0.282 for the social interaction anxiety factor and 0.224 for the social performance anxiety factor. The general factor’s $\omega$_H was above the acceptable cut-off score of 0.500, while the $\omega$_H⁢S for the specific factors was moderate. The ECV for the general and social interaction anxiety factors were 0.682, 0.159, and 0.158, respectively. Thus, the general factors of the SIAS-3 and SPS-6 explained significantly more of the common variance than the two specific factors of the SIAS-3 and SPS-6 in the bifactor CFA model. However, because the ECV of the general factor did not exceed 0.70, the bifactor CFA model could not be considered unidimensional.

As shown in Table 4, the $\omega$_H for the general factor (12 items) was 0.795, while the $\omega$_H⁢S for the social interaction anxiety factor was 0.206 and for the performance anxiety factor was 0.202. The general factor’s $\omega$_H exceeded the acceptable cut-off score of 0.500 and the $\omega$_H⁢S of the specific factors was moderate. The ECV of the general and social interaction anxiety factors were 0.681, 0.168, and 0.151, respectively. Thus, the general factor of the SIAS-6 and SPS-6 explained approximately two-thirds of the common variance, while one-third of the common variance was accounted for by the two specific factors of the SIAS-6 and SPS-6 in the bifactor ESEM model. The estimated ECV for the general factor was below 0.70, indicating that the bifactor ESEM model could be considered multidimensional rather than unidimensional.

Following the factor analysis results, the measurement invariance analyses focused on the bifactor CFA and bifactor ESEM models, as these provided the best fit for the data. Measurement invariance was examined using multigroup analyses across genders. The findings of the measurement invariance analyses for the bifactor CFA model and for the bifactor ESEM model are shown in Table 5. As seen in the table, for bifactor CFA model, the changes in CFI and RMSEA from configural to metric were –0.004 and 0.001, respectively, which is suitable for invariance, considering Chen’s guidelines [30] for CFI ($\le$–0.010) and RMSEA ($\ge$0.015) changes. The differences in RMSEA and CFI between metric to scholar invariance were $\mathrm{\Delta }$RMSEA = 0.002 and $\mathrm{\Delta }$CFI = –0.004 for bifactor CFA model. These differences were also considered invariant.

Additionally, for the bifactor ESEM model, the changes in CFI and RMSEA from configural to metric were 0.020 and –0.035, respectively, which is suitable for invariance, considering Chen’s guidelines [30] for RMSEA ($\ge$0.015) and CFI changes ($\le$–0.010). The changes in RMSEA and CFI from metric to scholar invariance were $\mathrm{\Delta }$RMSEA = 0.005 and $\mathrm{\Delta }$CFI = –0.007 for the bifactor ESEM model. These changes can also be considered as invariant. Thus, little change in model fit was observed between the different models, indicating overall measurement invariance across genders.