Why Ordinal Variables Can (Almost) Always Be Treated as Continuous Variables: Clarifying Assumptions of Robust Continuous and Ordinal Factor Analysis Estimation Methods

The analysis of factor structures is one of the most critical psychometric applications. Frequently, variables (i.e., items or indicators) resulting from questionnaires using ordinal items with 2–7 categories are used. There are plenty of articles that recommend treating ordinal variables in a factor analysis by default as ordinal and not as continuous imposing a multivariate normal distribution assumption. In this article, we exhibit that the reasoning behind such suggestions is flawed. In our view, findings from simulation studies cannot tell about the right modeling strategy of ordinal variables in factor analysis. Moreover, it is argued that ordinal factor models impose a normality assumption for underlying continuous variables, which might also often be incorrect in empirical applications. However, researchers seldom opt for more flexible modeling strategies that involve correctly specified distributions. Finally, the consequences of modeling choices for validity, reliability, measurement invariance, handling of missing data, and the assessment of global model fit are discussed.


INTRODUCTION
The analysis of factor structures is one of the most critical psychometric applications. Frequently, variables (i.e., items, or indicators) resulting from questionnaires are analyzed. These variables are often assessed with Likert scales that have a finite number of ordered categories. In many applications, the number of categories ranges between 2 and 7. An often posed question by applied researchers is about the most favorable approach for factor analysis in the presence of ordinal variables. First, ordinal variables could be treated as in the case of continuous variables, and the same estimation method would be used. Second, a factor model based on a distributional assumption for ordinal variables could be fitted (i.e., an ordinal factor model). There is a diversity of methodological literature addressing this issue (Dolan, 1994;DiStefano, 2002;Lei, 2009;Flora et al., 2012;Rhemtulla et al., 2012;Sass et al., 2014;Barendse et al., 2015;Asún et al., 2016;Li, 2016;Jia and Wu, 2019). The main message of most of the papers seems to be that ordinal variables should be treated as ordinal (i.e., not being treated as continuous variables) if there are only a few categories or the frequency distributions of variables are skewed (e.g., Rhemtulla et al., 2012). In this article, it is argued that treating ordinal variables as continuous can (almost) always be defended, and the argument for doing so does neither depend on the number of categories nor the marginal distribution of ordinal variables.
The article is structured as follows. In section 2, the two major competing estimation approaches are contrasted, and their assumptions are clarified. In section 3, it is argued that there are many options for modeling ordinal variables, and the choice of the ordinal method that is recommended by default is as arbitrary as treating ordinal variables as continuous variables. Finally, in section 4, the consequences of modeling choices for validity, reliability, measurement invariance, and the assessment of global model fit are discussed.

Alternative Modeling Strategies
Assume that a vector of ordinal variables Y = (Y 1 , . . . , Y I ) is given. For simplicity, we assume that each variable Y i has values 0, 1, . . . , K. Let p i,k = P(Y i ≤ k) denotes the cumulative frequencies of variable i. If the ordinal variables would be treated as continuous, a linear factor model is assumed, where µ are intercepts, is a loading matrix, F is a multidimensional factor variable, and E denotes a vector of residuals. Assume that E(F) = 0 and E(E) = 0. Let us define = Var(F) and = Var(E). Typically, is a diagonal matrix. The covariance matrix = Var(Y) in the factor model given in Equation (1) is modeled as an implied covariance matrix 0 (θ ) Hence, in (2), it is assumed that observed covariances of variables are represented by model parameters (i.e., loadings, factor covariances and variances, and residual variances). For ease of exposition, we interpret as a Pearson product-moment correlation matrix, which is the covariance matrix if ordinal variables Y would have been standardized prior to analysis. Let θ denote the vector of parameters in , , and that are freely estimated.
Two estimation methods for estimating θ can be distinguished (see Jöreskog, 2007, for an overview). First, a multivariate normal distribution can be assumed. Then, the estimated Pearson correlation matrix S for variables Y is a sufficient statistic for θ , and the fitting function (i.e., log-likelihood function in maximum likelihood estimation) is given as The multivariate normal distribution will be misspecified if variables are ordinal. However, ML parameter estimates are consistent and converge to a parameter θ that maximizes the Kullback-Leibler information (White, 1982). The optimal parameter θ is obtained if S is replaced by the true population covariance matrix in (3) (Arminger and Schoenberg, 1989;Olsson et al., 2000;Yuan and Bentler, 2007). Note that the model assumption = 0 (θ ) can be correct when the data is not multivariate normally distributed. The choice of a misspecified ML function F ML must not necessarily result in inconsistent (and, hence, biased) parameter estimates. Moreover, the socalled robust ML estimator (MLR; Savalei, 2014;Yuan and Bentler, 2007) provides valid statistical inference in a misspecified model (Satorra, 1992), and improvements have been proposed (Lai, 2019). Second, weighted least squares estimation based on the estimated correlation matrix S can be employed. Here, we consider diagonally weighted least squares (DWLS), and the fitting function is given as (Jöreskog, 2007) where W is a diagonal weighting matrix. MacCallum et al. (2007) argue that using DWLS instead of ML possesses advantages in the case of misspecified models [i.e., = 0 (θ )]. Alternatively, ordinal variables can be modeled with ordinal marginal distributions Yang-Wallentin et al., 2010). The idea is that there are underlying normally distributed variables The thresholds are given as τ ik = F −1 (p i,k+1 ), where F denotes the standard normal distribution function. In an ordinal treatment of ordinal variables, a linear factor model for the underlying latent normally distributed variables Y * (referred to as latent normality in the sequel) is assumed: Then, the covariance structure of * = Var(Y * ) of latent variables is modeled as * ≈ * 0 (θ * ) = * * * T + * The covariance matrix * can be estimated by employing polychoric correlations (Muthén, 1984). Note that the estimated parameter θ in (2) when treating variables as continuous will typically be different from θ * in (7) when treating variables as ordinal, even if original variables Y would have been standardized prior to analysis (Rhemtulla et al., 2012). Often, DWLS estimation based on the estimated polychoric correlation matrix S * is conducted (Muthén, 1984;Yang-Wallentin et al., 2010): It should be emphasized that the ordinal factor models are also labeled as item response models (IRT). The assumption of latent normality corresponds to the graded response model with a probit link function (Takane and de Leeuw, 1987;Glockner-Rist and Hoijtink, 2003;Kamata and Bauer, 2008). In standard IRT software, ML estimation is typically utilized for estimation by default (e.g., in the R package mirt; Chalmers, 2012) while in standard SEM software, limited information methods, such as DWLS are employed (e.g., in the R package lavaan; Rosseel, 2012). Practitioners often seek advice from methodologists of how to analyze ordinal variables in factor analysis. This amounts to the question of whether Pearson correlations or polychoric correlations should be used in the factor analysis. In the following, we critically discuss the recommendations of some often cited methodological articles.

There Is No "Correct" Modeling Strategy
Many papers recommend that ordinal variables should be modeled as ordinal if there are only a few categories or if the frequency distributions are asymmetrical (Rhemtulla et al., 2012;Li, 2016). These articles often state that biased estimates (i.e., factor loadings) would have been obtained if the (misspecified) model assuming multivariate normality would be applied. Their reasoning is based on simulation studies. We now argue that simulation studies do not help for providing rationale of the appropriate modeling strategy because it would be relatively simple to design simulation studies with ordinal variables that fulfill the linear factor model for Pearson correlations (i.e., Equation 2) instead of the factor model for polychoric correlations (i.e., Equation 7).
Researchers who claim the superiority of ordinal models for ordinal variables presuppose that Equation (7) holds, i.e., the matrix representation * = * * * T + holds for the polychoric correlation matrix (e.g., Li, 2016;Rhemtulla et al., 2012). Hence, they base inference on the latent normal variables Y * . By employing the corresponding estimation function F cat−DWLS , unbiased parameter estimates are obtained. However, if the variables would be analyzed as continuous (i.e., using F ML or F DWLS ), biased parameter estimates would be obtained. This is almost trivial as the true model is ordinal, and the question is whether the incorrect model that treats variables as continuous can also recover true parameters (i.e., whether θ * = θ holds in the population). A typical argument found in several articles is as follows. Under the assumption that continuous variables underlie the observed ordinal variables, the matrix of Pearson correlations underestimates of the correlation matrix among the underlying continuous variables, that is the polychoric correlation matrix (Olsson, 1979;Rhemtulla et al., 2012). Hence, using the Pearson correlation constitutes the incorrect input of the factor analysis, while polychoric correlations would be the correct one, and only the ordinal factor analysis would provide unbiased parameter estimates.
We now sketch the design of a simulation that "shows" that treating ordinal variables as continuous (i.e., using Pearson correlations) result in unbiased estimation while treating them as ordinal in the model assuming latent normality will result in biased estimates. Assume predetermined thresholds τ i,k (k = 0, 1, . . . , K). The Pearson correlation σ ij of variables i and j is a function of thresholds τ i,k and τ i,j , and the polychoric correlation σ * ij . One can always find a polychoric correlation σ * ij such that the covariance of i and j equals σ ij (e.g., by applying some numerical procedure for finding a root). Moreover, assume that the linear factor model = T + holds for observed covariances (i.e., Equation 2 holds true). Therefore, the ordinal factor model (Equation 7) will not be fulfilled, and hence, treating ordinal variables in an ordinal factor model assuming latent normality will result in biased estimates. Of course, a simulation exercise could be additionally carried to demonstrate this reasoning.
One could also argue that true scores (and consequently latent variables) are clearly defined in stochastic measurement theory (Steyer, 1989) as the expected value of an intraindividual distribution of an ordinal variable (see also Holland, 1990). Although observed variables are ordinal, item-specific true scores are bounded but non-integers. A factor model constitutes a model assumption for these true scores.
What can be learned from these observations? At best, findings from the literature that comes with recommendations to practitioners are only useful in identifying data constellations in which the continuous and the ordinal treatment of ordinal variables can result in similar parameter estimates. In our view, simulation studies or empirical data cannot be employed for deciding among the two competitive modeling strategies. Hence, a researcher must decide whether the factor structure should be posed on Pearson correlations or polychoric correlations.
As a cautionary note, one should add that categories of variables are stretched in the ordinal variable model according to their empirical frequencies when representing the latent variables F in the factor model. In contrast, by treating the variables as continuous, no such implicit transformation is carried out. The modeling choice relates to the question about the meaning of distances between categories of an ordinal variable. While it can be almost always be argued that assuming equal distances seem to be implausible in practice, distances that are derived on a sole empirical basis are equally implausible.
As a conclusion, we would not like to argue that dichotomous variables should be treated as continuous (see Maydeu-Olivares, 2005;Tran and Formann, 2010; for differences in parameter estimates for the two modeling strategies). However, we think that for items with 3-6 categories, using the linear factor model by treating variables as continuous is as defensible as treating them as ordinal. However, researchers should be aware that when estimating factor models with a misspecified distribution, the statistical inference should be obtained with the MLR estimator.

THE NORMALITY ASSUMPTION AND THE LATENT NORMALITY ASSUMPTION ARE EQUALLY RESTRICTIVE
In section 2, it was assumed that underlying continuous variables of the ordinal variables are normally distributed (i.e., latent normality holds). Typically, the latent normality assumption has often been taken for granted by applied researchers (Foldnes and Grønneberg, 2020). If there are violations of latent normality, parameter estimates based on the incorrect latent normality assumption can provide substantially biased estimates (Jin and Yang-Wallentin, 2017;Foldnes and Grønneberg, 2020). It has shown that the latent normality assumption can be empirically tested Raykov and Marcoulides, 2015;Foldnes and Grønneberg, 2020). However, it seems that these tests are seldom conducted in practice.
If the ordinal nature of variables would be taken seriously, more sophisticated modeling strategies for factor analysis that try to estimate more flexible distributions are required (Jin and Yang-Wallentin, 2017;Foldnes and Grønneberg, 2019). A particularly attractive distribution class is the factor copula model (Krupskii and Joe, 2013;Nikoloulopoulos and Joe, 2015;Ackerer and Vatter, 2017;Krupskii and Genton, 2018). Copula models decompose a joint distribution for modeling into marginal distributions and modeling the dependency structure. Gaussian copula models pose a multivariate normal distribution for modeling the dependency structure while allowing a semiparametric estimation of the marginal distribution (Hoff, 2007;Gruhl et al., 2013;Murray et al., 2013). As a consequence, underlying latent variables can deviate from the latent normality assumption. More formally, it is assumed that there exists a vector of multivariate normally distributed variables Y * with Var(Y * ) = * . For an ordinal variable Y i , there exists an underlying latent variableỸ i = G −1 i (F(Y * i )), where G i denotes the distribution function of variableỸ i , and F is the normal distribution function. Like in Equation (5), the ordinal variable Y i is obtained by discretizing the underlying continuous variablẽ Y i with respect to thresholds τ i,k . If latent normality is fulfilled, it holds thatỸ i = Y * i . For example, the distributions G i could be, for example, the logistic, skew normal, skew t, or cloglog distribution. These distributions could be fixed or estimated using empirical data (Gruhl et al., 2013).
It should be noted that there is active research in factor analysis for continuous variables with non-normally distributed factors or residuals that can be skewed, bounded, or are mixtures of distributions (Song et al., 2010;Kelava and Brandt, 2014;Zhang et al., 2014;Asparouhov and Muthén, 2016;Lin et al., 2016;Revuelta et al., 2020). Using these more complex distributions would reduce the degree of distributional misspecification in the factor model. Moreover, in a few articles, the estimation of the link functions in item response models is considered (Peress, 2012;Liang and Browne, 2015;Feuerstahler, 2019). Notably, estimating the link function in IRT is equivalent to estimating the marginal distributions of the underlying latent variables Y * i in a Gaussian copula model. To sum up, we believe that also estimating the marginal latent distributions of underlying continuous variables adds a further layer of complexity regarding estimation and interpretation in the analysis. However, if scholars suggest to always model ordinal variables appropriately by a well-fitting model, there is no reason for only fitting the ordinal model based on latent normality. In this regard, multivariate latent normality is a testable assumption as multivariate normality of the ordinal variables, and it can be supposed that both assumptions will typically be violated in practice. As a consequence, both treatments, the continuous (i.e., using Pearson correlations ) and the ordinal (i.e., using polychoric correlations * ), correspond to misspecified models, and it is difficult to speculate about a plausible data-generating model in a concrete application. It can be argued that by choosing a particular fitting function, the parameter of interest is defined, and one could consider the parameter estimate in a sample as an estimator that converges to some optimal parameter in the fitting procedure that would have been obtained in an infinite sample size. As it becomes clear in this case, again, simulation studies cannot decide about choosing an adequate modeling strategy (i.e., choosing a fitting function).

DISCUSSION
In this article, it is argued that the often found recommendations for not treating ordinal variables in factor models as continuous are not justified. The choice for a particular modeling strategy implies that it is assumed whether the linear factor model must hold for Pearson correlations (i.e., the normality assumption) or polychoric correlations (i.e., the latent normality assumption). This choice cannot be derived from an empirical dataset or simulation studies, although some articles argue otherwise.
It should be mentioned that alternative modeling choices are rarely discussed under the perspective of validity (but see Ferrando, 1999, for an exception). Using ordinal factor models implies the use of a non-linear scoring formula for factor scores, which is not the case when variables are treated as continuous. It has to be defended from a validity perspective (i.e., by studying relationships with external criterion variables) why the non-linear scoring rule from an ordinal factor model is preferable to a linear scoring rule from the factor model using Pearson correlations.
We would also like to emphasize that the computation of a reliability measure (e.g., ω) is defensible if ordinal variables would be treated as continuous by analyzing Pearson correlations in factor analysis (see also Chalmers, 2018). As it has been pointed out by Lucke (2005), the reliability of a sum score can be defined on a matrix decomposition into a part for the true score (i.e., T ) and a part for the error (i.e., ). Unbiased estimation of reliability only needs the unbiased estimation of the model parameters that are involved in the computation. It is not required that the factor model is correctly specified with respect to distributional assumptions of variables. In our view, the assessment of reliability based on an ordinal factor model (Green and Yang, 2009) is not necessarily superior to using a factor model that treats ordinal variables as continuous. Moreover, we generally question the strong preference of a model-based assessment of reliability from a factor analysis (e.g., McNeish, 2018;Sijtsma, 2009;Yang and Green, 2011) compared to designbased reliability measures (i.e., internal consistency, see Meyer, 2010, for assumptions of Cronbach's alpha).
In factor analysis, the investigation of measurement invariance (Millsap, 2011) is heavily discussed. It has been recommended that invariance analysis for ordinal variables should also treat variables as ordinal (Chen et al., 2020;Svetina et al., 2020). By continuing our arguments, we think that the assessment of invariance can be equally defended by treating ordinal variables as continuous. Effect sizes of non-invariance could be defined in the metric of raw scores (i.e., based on means and covariances; see Gunn et al., 2020).
The assessment of the global model fit in ordinal factor models utilizing fit statistics is more challenging than in models that presuppose multivariate normality. Recently, it has been proposed to define fit statistics for ordinal factor models that have the same rationale as in the normality case by replacing Pearson correlations (i.e., matrix ) by polychoric correlations (i.e., matrix * ) in the definition of population effect sizes of global model fit (Savalei, 2020). Pursuing such a strategy implies that the latent normality assumption is taken for granted, and potential misspecification of latent normality is not quantified in the assessment of model fit. Hence, we think that the assessment of model fit should rely on bivariate distributions of observed variables (Maydeu-Olivares, 2013). In this respect, the assessment of model fit is independent of using a particular fitting function, a property which, however, might see some scholars as a disadvantage.
We argued that a researcher might be interested in interpreting model parameters of a factor model with a misspecified distribution. However, in the presence of missing data, fitting a misspecified factor model with ML or weighted least squares approaches only result in consistent parameter estimates if data are missing completely at random (MCAR; see Yuan, 2009). In the case of missing at random (MAR) data, in general, inconsistent and biased parameter estimates will be obtained except for linearly related variables (Yuan, 2009;Yuan and Bentler, 2010;Yuan et al., 2012). For ordinal variables, linear relations among variables are implausible. However, it could be argued that there are typically only mild deviations of the MCAR assumption in empirical datasets (Newman, 2014). If substantial deviations of MCAR are suspected, and a misspecified factor model should be estimated, imputation based approaches might be preferable (Gottschall et al., 2012;Jia and Wu, 2019). The use of sufficiently complex imputation models, such as the Gaussian copula model (Hollenbach et al., 2018), mixture models (Murray and Reiter, 2016), or latent class models (Vermunt et al., 2008;Si and Reiter, 2013) are advantageous to minimize possible distributional misspecifications for MAR data. Appropriate imputation models can also treat specific deviations from MAR (missing not at random; MNAR; Harel and Schafer, 2009;Jung et al., 2011;Kano and Takai, 2011;Zhang and Reiser, 2015;Bartolucci et al., 2018;Kuha et al., 2018;Pohl and Becker, 2020).
As pointed out by a reviewer, the application of factor models based on the normal distribution requires fewer methodological skills than an ordinal factor model. Forcing the application of the more sophisticated ordinal factor model might have a biasing effect against research produced in developing countries because, in those countries, the use of factor models under the normal distribution still prevails. The reviewer remarked that our recommendation for also considering the more straightforward approach might enhance inclusion and promote diversity in academia. This aspect was indeed not the primary motivation for writing this article, but it could be a pleasant side effect.
To conclude, we tried to argue that simple suggestions derived from simulation studies cannot tell what the appropriate modeling strategy for handling ordinal variables in the factor analysis is. It is maybe improbable to reach a consensus about this issue. We tend to prefer the modeling approach that possesses higher validity. As a comprehensive assessment of validity cannot be reduced to a single quantitative measure, the choice of an appropriate factor modeling approach is not a purely statistical issue. However, it must be defended by a researcher with a concrete question at hand. Finally, we would like to emphasize that modeling ordinal variables with Pearson correlations as well as polychoric correlations can provide complementary information. Moreover, the results of alternative model specifications are worth to be reported (see also Steegen et al., 2016;Hoffmann et al., 2020).

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.