The linear factor model (FA) for a vector of variables y* of J elements is

where ξ is a vector of latent factor, e is a vector of random measurement errors¹, and Λ is the matrix of factor loadings. Under the additional assumption of normality, the factors, ξ, follow a standard normal distribution, and the errors, e, are distributed as normal (0, φ), where φ is the standard deviation. In consequence, the distribution of y* is multivariate normal with variance-covariance matrix:

where Ψ is the diagonal variance-covariance matrix of e.

The categorical factor analysis model (FAC) assumes that y* is an unobservable variable and that the manifest responses are obtained from y*through a discretization process (Muthén, 1984; Ferrando and Lorenzo-Seva, 2013). For example, the binary factor model assumes that the response to item j is 1 if yj* is sufficiently high and 0 otherwise. More specifically:

where τ_j is a threshold parameter. In general, a model for K ordered response categories is based on a discretization process for y* based on K-1 threshold parameters.

Mathematically, the skewness of the distribution of the manifest variable, y_j, depends on the placement of the thresholds. In the case of dichotomous data, the distribution will be skewed when the threshold parameters have an extreme value. On the contrary, the distribution of responses will be roughly symmetric when the threshold (the population mean of ξ) is close to zero. If the threshold is far away from zero, one of the categories will be more probable than the other. In the general case of an item with K response categories, skewness is associated to an uneven location of thresholds around the factor mean.

In general, the variance-covariance matrix of y, say Σ, will not conform to the structure implied by the linear model of Equation (2). However, Muthén and Kaplan (1985) identified some circumstances where the covariance structure of (2) holds for both y* and y. In particular, when a single factor model holds for y* and all the items have the same slope, error variance and thresholds, all correlations between y* are equal. In this case, all the pairs of y_j have the same correlation, ρy=a⁢ρy*, where a is the attenuation factor.

Figure 1 provides details about the effect of the skewness on the Pearson correlation (Jorgensen and Johnson, 2022), showing the correlation for two variables that measure a common factor. Factor loadings are set to 0.71 and the correlation between y₁ and y₂ is computed as a function of the threshold.

The most evident result in Figure 1 is that the correlation depends on the threshold and is higher for threshold values close to zero. This is a noteworthy difference between FAC and FA, because in the later model the correlation is mathematically independent of the mean of the variables. A more subtle result shown by Figure 1 is that when the variables are skewed in opposite directions, the attenuation of the correlation is more pronounced and the overfactoring will be potentially more severe. This can be appreciated by comparing the correlations in the corners of Figure 1. For example, the correlations in the corner (−2, −2) are about 0.20 whereas, in the corner (−2, 2), they are about 0.11.

The maximum-likelihood estimation of the FA model consists of finding the values of the parameters that minimize the discrepancy function (Browne and Cudeck, 1992).

where Σ is a variance-covariance matrix and p is the number of manifest variables. The population value of the discrepancy is computed by setting Σ to the true variance-covariance matrix (for example, those reproduced by a FAC model) and minimizing F with respect to Σ*. In estimation, Σ is set to the sample variance-covariance matrix (Σ = S). The likelihood-ratio chi-square goodness of fit statistic is X2=NF(S,Σ)*, and E(X²) = df (Cudeck and Browne, 1985).

Fitting a FA model to a population that conforms to a FAC model has two unwilling consequences. First, in the special circumstances where the FA fits the true variance-covariance matrix [that is, F(Σ,Σ)*=0 in the population], the FA contains the correct number of latent factors but the slopes will be attenuated. Second, in most cases the FA do not fit the population model [the population value is F(Σ,Σ)*>0]; in that circumstances, the FA cannot reproduce the true variance-covariance matrix, the chi-square has an expectation of E(X²) = δ + df (where δ is a positive non-centrality parameter), and the probability of rejecting the FA model with the correct number of factors will be inflated. The psychometric literature has demonstrated that this bias may lead to distortions such as overfactoring and the emergence of difficulty factors (McDonald and Ahlawat, 1974).

This article compares several widely used approaches to estimate the number of latent factors from ordinal data. These approaches consist of comparing the chi-square associated with a particular model (FA or FAC) and selecting the most parsimonious model that obtains a non-significant statistic. The alternatives are the likelihood-ratio chi-square for the FA model, the modified chi-square statistics that compensate for the lack of normality (Satorra and Bentler, 2010; Bryant and Satorra, 2012), the chi-square associated with FAC under two estimation methods [maximum likelihood (ML) and WLSMV], and the (GRM; Samejima, 1969) from IRT. The purpose is to evaluate the precision of the different methods and the potential bias associated with fitting a theoretically incorrect model (the FA). The comparison is conducted via simulation studies.

The simulated conditions include the sample size (N), the number of observed variables (p), the number of thresholds (τ), the magnitude of the factor loadings (λ), and the skewness as measured by the γ coefficient (Joanes and Gill, 1998). The levels of these independent variables are summarized in Table 1. In summary, the number of conditions examined was 120 = 2 (sample size levels) × 2 (number of observed variables levels) × 3 (number of thresholds levels) × 2 (loading levels in the factors) × 5 (skewness levels). The number of simulated samples for each condition is 1000.

We have manipulated N because of its direct relation to chi-square and the potential effect on the results. The number of factors may have an impact on the results of FA given that a model for continuous variables might better approximate the data when the number of response categories increases. The number of variables is manipulated because the size of the fitted model is another condition that affects the population value of the discrepancy function and might influence the results of the chi-square.

The simulating model is a common-factor model (one latent factor) for categorical variables that realistically represents data from rating scales. We have simulated only the one-factor model to keep the number of conditions at a manageable level. Because the true model contains a single factor, we conclude that overfactoring occurs whenever the number of estimated factors is more than one.

The following models were fitted to each sample:

The R code for fitting these models appears in Appendix E of the Supplementary materials.

The procedure for simulating categorical data with a prescribed value of skewness is due to Olsson (1979). We are working in the specific case of ordinal data with binomial marginals. The binominal marginal has one free parameter, π, which can be chosen to get the desired skewness. Let γ be the coefficient of skewness and r be the number of thresholds. Using Equation (7) in Olsson (1979), the parameter π of the binomial distribution with given values of r and γ is:

Once that π has been computed from Equation (5), the second step of the procedure consists of computing the cumulative probabilities of the binomial (r, π) distribution, say pr = [pr₁, pr₂, …, pr_(r+1)]. Finally, the thresholds τ₁, …, τ_k are the standard normal deviates associated to the probabilities given by pr; that is, the thresholds are τ_t = F⁻¹(pr_t) for t = 1, …, r.

Table 2 shows the threshold values associated with the different values of γ and r manipulated in the simulation study. The same thresholds were used for all the variables that have the same value of γ. Notice that it is not guaranteed that the set of thresholds associated with a given γ is unique. Other algorithms may remove the assumption of binomial marginals and render different thresholds, which could have an impact on the results of the simulation. The R code for implementing this procedure is given in Appendix E.

In all simulated conditions, the population model is FAC. However, there are conditions where the FA is capable of reproducing the true population variance-covariance matrix and there are others where it is not. As explained in Section “Simulation study,” the FA can reproduce the variance-covariance matrix when all the items have the same parameters (slopes, error variance, and thresholds; that is, in the conditions with no mixed skewness because mixed skewness means that the items have different thresholds). In the conditions with no mixed skewness, the FA and the FAC models are equivalent in the population, the two models have the same number of factors, and attenuation coefficients for the Pearson correlations (Muthén and Kaplan, 1985) take the values shown in Supplementary Table 1. In these cases, any discrepancy between the number of factors estimated by FA and FAC is not attributable to the population model and the simulation provides information about the relative performance of the estimation methods.

In the conditions with mixed skewness (mild and mixed skewness, strong and mixed, none + strong, and mixed), the maximum-likelihood discrepancy function for the FA model attains a nonzero value in the population. The population value of the discrepancy and the RMSEA (Browne and Cudeck, 1992) appear in Supplementary Table 2. The consequence is an inflated chi-square statistic and an increased probability of rejecting a FA model that contains the correct number of factors. In these cases, the simulations inform about the magnitude of the overfactoring effect due to fitting an incorrect model (the FA). Notice that the conditions of mixed skewness are more realistic and representative of practical applications, where the variables rarely are perfectly parallel measures.

The data were analyzed using the empirical proportion of samples (EPS) in which the one-factor model is retained for each condition. The logit of the p-value was analyzed as a function of the independent variables manipulated in the study using an ANOVA model to compute the effect size (partial-η²) associated with each condition. Notice that the EPS is the proportion of correct decisions, and 1—EPS is the proportion of samples in which the latent dimensionality is overestimated. Thus, the EPS can be interpreted as the statistical specificity of the testing procedure and 1—EPS as the Type-I statistical error.

All the analyses were conducted using the R programming language (R Core Team, 2021) and the aforementioned packages. Some other R packages were used for the analysis of the results: the skewness and kurtosis were estimated using e1071 (Meyer et al., 2021), and the ANOVA table of the results of the simulation was computed using sjstats (Lüdecke, 2021). All the methods used in this article are included in the R language and the indicated packages, so they are readily available and free of charge for any practitioner.

The success of a method cannot be based solely on the estimation of the number of factors. A method that correctly estimates the latent dimensionality whilst producing wildly incorrect parameter estimates is of limited usefulness. Parameter estimates under the different models were analyzed to assess their relative merit. Recovery was evaluated by the root mean squared error between the true and estimated parameters (RMSE).

The estimated parameters for FA have two specific sources of error that are not present in the other methods: Attenuation and structural error. As explained previously, attenuation occurs when a FA model is able to reproduce the true variance-covariance parameters (which conform to a FAC model) but the FA parameters have smaller values than the true FAC model. Attenuation occurs in conditions with no skewness and strong positive skewness, and its magnitude depends on the coefficients shown in Supplementary Table 1. In the conditions where the FA cannot reproduce the population variance-covariance matrix of the FAC model (those that involve mixed skewness), the FA estimates have the additional bias of structural error. These are population biases that will not vanish no matter what the sampling size is. The simulations are informative about the magnitude of the three sources of error combined (sampling error, attenuation, and structural error), and the increase of error that is specific to FA in comparison to the other methods.

The mirt function estimates the GRM using the IRT parameterization, whereas the true parameter values for the simulation are in the FA parameterization. Previously to computing the RMSE, the IRT parameters were converted to the FA parameterization using the equations (A15) and (A16) of Paek et al. (2018).

The recovery of lambda in relation to the manipulated conditions appears in Supplementary Table 8 and Supplementary Figures 4–6. WLSMV and FAC performed similarly and obtained the higher RMSE in the conditions of a large number of thresholds, small loadings, and strong positive skewness. The results indicate that GRM outperforms these methods and has a smaller RMSE that remains stable across conditions. The FA is affected by attenuation and inconsistency biases and is generally poorer than the other methods. The results for ψ follow a similar pattern and are summarized in Supplementary Table 9.

Recovery of τ is summarized in Supplementary Table 10 and Supplementary Figures 7–9. The results are similar to those found for the recovery of λ. There are small differences between WLSMV and FAC, and the two methods render an increased RMSE in the conditions of a large number of thresholds and strong positive skewness. The GRM obtains a smaller RMSE than the other methods and is stable across conditions.

The results of our simulation have shown that, in the present context, the chi-square statistic for FA leads to overfactoring. However, there are remarkable differences between the models and the study conditions. The magnitude of the factor loadings and the skewness of the factor indicators have the largest effect size in relation to the logit of the p-value, but this effect varies between the models. The empirical proportion of rejection for the one-factor FA model is generally high and sharply increases in the conditions where this model is wrong in the population. The PA shows better performance than the chi-square, but the overfactoring is still high in the presence of skewness. Notice that the effect of skewness might depend on the choice of thresholds values. We have used the algorithm by Olsson (1979) to generate the threshold values associated with the gamma coefficient. However, other algorithms may render different thresholds associated with the same gamma and might have a different effect on the performance of the chi-square. One topic for further investigation would be the development of other algorithms for generating thresholds and the study of the differential impact of the choice of thresholds while keeping gamma fixed to a constant value. The pattern of the factor loadings also has a noticeable effect, and smaller rejection rates are found for the tau-equivalent than for the congeneric model.

Interestingly, dimensionality recovery for the FAC model is tightly related to the estimation method, as recovery is poor for ML and reliable for WLSMV and dT_M.

The effect sizes of all independent variables are negligible for the GRM, which indicates that the model is robust to the manipulated conditions. The Type-I error is about 30%, regardless of skewness, which is still far from the nominal levels, but it improves the results of the other models. The number of factor indicators also has an effect, and more rejections of a correct model are found with a larger number of factor indicators. This is an expected result because bigger samples are necessary for the likelihood-ration statistic to approximate its chi-square reference distribution when there are a large number of factor indicators. On the counter side, the small effect size associated with N means that the improvement in EPS decreases as N increases. From a theoretical point of view, the GRM is not exactly true in this simulation because data were generated assuming a normal distribution for the error (normal link) whereas the implementation of the GRM in the mirt package assumes a logit link. The logit link was developed as a convenient approximation to the normal link and the difference between them is small. However, this difference does not vanish as sample size increases, which converts the logit link to an inconsistent procedure. One problem for future research is to investigate if the results of the GRM improve when estimated with a normal link.

Finally, the WLSMV and the Satorra-Bentler dT_M are the best alternatives to estimate dimensionality, with rejection rates close to the nominal 5% level that remains relatively constant across all the simulated conditions. The recovery of parameters for GRM was better than for WLSML, and less affected by those conditions in which the RMSE increases for WLSMV. Thus, a good combination of methods could be to estimate the dimensionality using a Satorra-Bentler fit statistic, as in WLSMV, and estimate parameters using a marginal ML algorithm as in the GRM once the dimensionality of the model has been fixed.

A real data sample was analyzed to illustrate the different approaches in an applied context. The data consists of responses to the scale of attitudes toward censorship, originally published by Rosander and Thurstone (1931), and reprinted by Shaw and Wright (1967). The scale is composed of 20 Likert-type items with six response alternatives. The responses are labeled: Strongly Disagree, Disagree, Slightly Disagree, Slightly Agree, Agree, and Strongly Agree. The sample size is 223 observations. The items and the data are publicly available in the Georgia Tech Psychometric Research and Development Lab² and were retrieved on 03/17/2022. We have applied the same methods used in the simulation study (FA, FAC, WLSMV, dT_M, and GRM) and some new methods for completeness. The new methods are:

The R codes for all the analyses appear in Appendix G.

In the simulation study, we have assumed that the population model is an ordinal factor model and we have manipulated the parameters of the model to investigate the performance of several estimation methods, but the correctness of the model has not been put into question. In a real data analysis, the validity of an ordinal factor model cannot be taken for granted as it might not be correct in the population, which introduces another source of bias in the analysis. Some of the above analyses (normality of latent responses, normality of latent factors, and the MNCM) evaluate part of the assumptions of the ordinal factor model to determine if they may be maintained for the present sample.

Table 9 contains the descriptive statistics for the items, including the response frequency of the categories, the mean, the standard deviation, and the kurtosis. The variables have a pattern of mild and mixed skewness. Most of them have a mean value close to the midpoint of the scale, and a skewness value close to zero. However, there are also several variables negatively skewed because of the mean being close to the higher point of the response scale, and one of the variables has a low mean and positive skew. The multivariate skewness and kurtosis coefficient have a p-value smaller than 0.01 and the hypothesis of normality is rejected.

The hypothesis of normality of the latent responses is discarded (p-value < 0.01). This result threats the validity of the methods based on polychoric correlations since these correlations are highly sensitive to the validity of the normality assumption for the latent responses (Foldnes and Grønneberg, 2020, 2022). However, as we do not have information about what the distribution of the latent responses might be, we are not in a position of estimating polychoric correlations under a different distribution.

Regarding the analysis of normality for the latent factor, the HQ for the normal, Davidian-2, Davidian-3, and empirical histogram models are 13782.5, 13789.3, 13788.2, and 14175.1, respectively. Similar to AIC and BIC, the HQ supports the model that minimizes it. Since the normal distribution minimizes the HQ, the Gaussian distribution can be retained.

The goodness-of-fit statistics for the FA models with one to four factors appear in Table 10, including the chi-square against a saturated model and the modified chi-squares: T_M and T_MV. The chi-square estimates three factors whereas the robust statistics estimates two factors. Since the data show non-negligible skewness and kurtosis, these results point in the direction of relying in robust statistics and retaining two factors. Parallel analysis gives a similar result and the number of estimated factors is two for PA and three for CPA.

Table 11 contains the goodness-of-fit statistics for FAC with the WLSMV estimator, the chi-square statistics corrected for non-normality (the T_M, the T_MV), and the Satorra-Bentler correction of the chi-square difference (dT_M). The statistics differ in their conclusions. The chi-square and the chi-square difference give support to a three-factor model. On the other hand, the modified chi-square statistics (T_M, T_MV, and dT_M) suggest models with higher dimensionality.

The FAC model under ML results in five or six estimated factors, depending on the statistic used to test model fit (chi-square of absolute model testing or chi-square between pairs of nested models). The results appear in Table 12. However, the FAC method based on ML apparently is questionable in light of the simulation results.

The MNCM was estimated to evaluate the assumption that responses follow an ordinal scale. The difference between an ordinal and a nominal model is that there is only one slope for each item in the former model and one slope for each response category in the later model (except for one of the categories of the item, which has a slope of zero for identifiability). The results appear in Table 13 and indicate that the MNCM does not significantly fit better than the GRM. Thus the ordinality assumption involved in the GRM can be maintained in light of this analysis.

This empirical example shows mixed results. The most parsimonious results are those of FA and GRM, whereas the modified chi-squares lead to models of higher dimensionality. The large number of factors resulting from polychoric-based methods (FAC in its different variants) may be due to an unprecise estimation of polychoric correlations because of the lack of normality of the latent responses. In the simulation, we considered only the case in which latent normality holds. This assumption is not tenable in the empirical application, and the psychometric literature shows that this might be a source of imprecision in the estimates of polychoric correlations (Grønneberg and Foldnes, 2022). All in all, we would recommend the solution with a small number of factors for parsimony and the questionability of polychoric-based methods in this sample.