Factor Analysis of Ordinal Items: Old Questions, Modern Solutions?

1. Introduction

In 1904, Charles Spearman introduced a groundbreaking statistical methodology aimed at uncovering latent variables, those not directly measurable but estimable through observed correlations. Spearman’s approach, inspired by the pioneering work of Galton and Pearson, sought to operationalize the existence of two intelligence factors inferred from children’s school performance. He posited a general intelligence factor (G factor) and a specific intelligence factor (S factor), suggesting that positive correlations across seemingly unrelated subjects indicated these underlying factors (Spearman, 1904a).

Although Spearman’s theory of intelligence factors faced challenges from researchers like Thomson (1916) and Thurstone ^(1938),, his method, Exploratory Factor Analysis, became a cornerstone in psychometry. This approach facilitated the estimation of intangible attributes, contributing significantly to the development of psychological assessment.

Around the same period, Rensis Likert (1932) introduced his influential scale to measure opinions and attitudes. This scale, composed of 7-point response items coded from 1 to 7, quickly gained popularity in psychology and remains widely used today for assessing latent psychometric variables.

Consequently, factor analysis is frequently employed to uncover the latent constructs assessed by Likert and similar rating scales. However, the validity of results obtained through factor analysis of ordinal items, which comprise a significant portion of psychometric scales used in daily psychological assessments, remains controversial. Ongoing debates, such as those discussed by Liddell and Kruschke (2018), underscore the need to consider both traditional and modern solutions when addressing the challenges associated with factor analysis of ordinal items.

1.1. Factor Analysis

Factor analysis aims to uncover latent variables through observed correlations. For instance, in the Beck Depression Inventory (BDI), responses to mood-related items reflect latent depression, with correlations among these items indicating its presence. Spearman’s factor analytic model, explains observed variables’ behavior through latent factors. Thurstone (1938) and Cattell (1943) expanded this model, introducing exploratory factor analysis (EFA) with a reduced set of common factors and specific factors for each variable. Here a set of observerved variables (X₁, X₂,… X_p) can be explained by a reduced m common factors (m<< p) and p specific factors as¹:

$$\begin{array}{l}{x}_{1j}={\mu }_1+{\lambda }_{11}{f}_{1j}+{\lambda }_{12}{f}_{2j}+\mathrm{...}+{\lambda }_{1k}{f}_{kj}+{\epsilon }_{1j}\\ x_{2j}={\mu }_2+{\lambda }_{21}{f}_{1j}+{\lambda }_{22}{f}_{2j}+\mathrm{...}+{\lambda }_{2k}{f}_{kj}+{\epsilon }_{2j}\\ ⋮\\ x_{pj}={\mu }_p+{\lambda }_{p1}{f}_{1j}+{\lambda }_{p2}{f}_{2j}+\mathrm{...}+{\lambda }_{pk}{f}_{kj}+{\epsilon }_{pj}\end{array}$$

(1)

where m_i represents the expected value of variable x_i, l_ik represents the loading of factor k (k ∈{1, …, m}) on variable i, f_kj represents the common factor k for individual j (j ∈{1,…, n}), and e_ij represents the residual or specific factor associated with observed variable i (i ∈{1, …, p}) on individual j . Common factors are assumed to be independent and specific factors are assumed independent and normally distributed.

The factor analysis can be done on the covariance matrix of the observed variables estimated as:

$$S=\left[\begin{array}{lll}{s}_{11}\hfill & \cdots \hfill & s_{1p}\hfill \\ \hspace{0.17em}⋮\hfill & \hspace{0.17em}⋮\hfill & \hspace{0.17em}⋮\hfill \\ s_{p1}\hfill & \cdots \hfill & s_{pp}\hfill \end{array}\right]$$

(2)

whose generic element is Galton’s covariance:

$$s_{ij}={\displaystyle \frac{{\displaystyle \sum _{l=1}^n(X_{il}-{\overline{X}}_i)(X_{jl}-{\overline{X}}_j)}}{n-1}}$$

(3)

Since variables in the social sciences and humanities often have different scales and ranges of measurement, it is common to standardize them, centering the mean to 0 and reducing the standard deviation to 1 as $z=(x-\overline{x})/s$. The variance-covariance matrix is then reduced to the matrix of the product-moment correlations, also known as Pearson’s correlations:

$$R=\left[\begin{array}{lll}1\hfill & \cdots \hfill & r_{1p}\hfill \\ \hspace{0.17em}⋮\hfill & \hspace{0.17em}⋮\hfill & \hspace{0.17em}⋮\hfill \\ r_{p1}\hfill & \cdots \hfill & 1\hfill \end{array}\right]$$

(4)

whose generic element is Pearson’s well-known correlation coefficient:

$$r_{ij}={\displaystyle \frac{{\displaystyle \sum _{l=1}^n(X_{il}-{\overline{X}}_i)(X_{jl}-{\overline{X}}_j)}}{\sqrt{{{\displaystyle \sum _{l=1}^n(X_{il}-{\overline{X}}_i)}}^2}\sqrt{{{\displaystyle \sum _{l=1}^n(X_{jl}-{\overline{X}}_j)}}^2}}}$$

(5)

Using standardized variables, model (1.1) can now be written in matrix form, as:

$$z=\Lambda f+\mathsf{\epsilon }$$

(6)

where L is a matrix of factor loadings, f is a vector of common factors, and e is a vector of specific factors. From model 1.7 the correlation matrix can then be written as:

$$R=\Lambda \Lambda ’+\psi$$

(7)

where y is the diagonal covariance matrix of the specific factors. Common factor loadings can be estimated from eq. 8. However, this equation is nonidentified², so Spearman, in his “method of correlations” followed Pearson’s Principal Components method for factorizing the correlation matrix. The variables’ communality (see eq. 2) are therefore estimated as the eigenvalues of the principal components retained from the correlation matrix. More recent methods of extraction of the common factors include the principal axis factoring method and the maximum likelihood method (de Winter and Dodou, 2011). In short, the principal axis factorization method considers that the initial communalities are the coefficient of determination of each variable regressed on all others (and not 1, as in the principal components method) proceeding with an analysis of principal components in a optimization algorithm that stops when the extracted communalities reach their maxima. The maximum likelihood method, as proposed by Lawley and Maxwell (1962), estimates the factor loadings that maximize the likelihood of sampling the matrix of observed correlations in the population. The solution is obtained by minimizing the maximum likelihood (ML) function³

$$F_{ML}=tr\left[{(\mathrm{\Lambda }\mathrm{\Lambda }\text{'}+{\psi }^2)}^{-1}R\right]-\log \left|{(\mathrm{\Lambda }\mathrm{\Lambda }\text{'}+{\psi }^2)}^{-1}R\right|-p$$

(8)

for  and y by a two-step algorithm that estimates the eigenvalues (common factor variance) and specific factor variances of the correlation matrix R.

Since the late ‘60s, the maximum likelihood method has been widely used in confirmatory factor analysis (CFA) (Jöreskog, 2007). In CFA, a predefined factor structure is validated, imposing theoretical constraints on the model, in contrast to the “free” estimation in EFA. The key difference between the two models lies in the specific loading of common factors on items in CFA, while EFA allows all common factors to potentially load on all items as illustrated in Figure 1.

This section provides a condensed overview of factor model parameter estimation methods, as detailed in Jöreskog (2007) and Marôco ^(2021b). Iterative techniques minimize a discrepancy function f (Equation 10) using sample variance-covariance matrix (S) and expected covariance matrix ($\Sigma (\widehat{\theta })$):

$$f=F(S,\mathsf{\Sigma }(\widehat{\theta }))$$

(9)

For maximum likelihood estimation (ML), under the normality assumption for the manifest variables, the discrepancy function to minimize is

$$f_{ML}=\log |\mathsf{\Sigma }(\widehat{\theta })|+\hspace{0.17em}\mathrm{tr}(S\mathsf{\Sigma }{(\widehat{\theta })}^{-1})-\log \left|S\right|-\hspace{0.17em}p$$

(10)

Non-normal data analysis recommends the Asymptotic Distribution Free (ADF) (Arbuckle, 2006) or Weighted Least Squares (WLS) (Li, 2014) estimators:

$$f_{WLS}=(S-\Sigma (\widehat{\theta }))\text{'}{W}^{-1}(S-\Sigma (\widehat{\theta }))$$

(11)

The generic element of W, calculated from the sample covariance elements and fourth-order moments of the observed x variables, is (Finney & DiStefano, 2013):

$$\begin{array}{l}{w}_{ij,kl}=s_{ijkl}-s_{ij}{s}_{kl}\\ s_{ijkl}={\displaystyle \frac{1}{n}}{\displaystyle \sum _{a=1}^n(x_{ai}-{\overline{x}}_i)(x_{aj}-{\overline{x}}_j)(x_{ak}-{\overline{x}}_k)(x_{al}-{\overline{x}}_l)}\end{array}$$

The WLS estimator, considering sample covariances and fourth-order moments, may face convergence issues with small samples (<200), addressed by Diagonally Weighted Least Squares (DWLS) or Weighted Least Squares Mean and Variance (WLSMV), both special cases of WLS. However, ignoring item covariances reduces statistical efficiency, impacting model fit. Ordinal factor analysis utilizes the polychoric correlation matrix, with WLS estimator incorporating sample polychoric correlations and thresholds to estimate the asymptotic covariance matrix (W) (Finney & DiStefano, 2013; Muthén, 1984):

$$f_{WLS}=(R-\hat{\mathsf{\rho }})\text{'}{W}^{-1}(R-\hat{\mathsf{\rho }})$$

(12)

where R is the matrix of the polychoric correlations and thresholds estimated from the sample; $\mathsf{\rho } ^$ is the model implied correlation and thresholds.

Stanley Stevens and the Scales of Measurement

During the 1930s-40s, Stanley Stevens introduced a taxonomy for “scales of measurement,” influenced by the work of a committee from the British Association for the Advancement of Science (N. Campbell, p. 340 of the Final Report, cit in Stevens, 1946). Stevens proposed four scales: Nominal, Ordinal, Interval, and Ratio, which unified qualitative and quantitative measurements. Nominal scales evaluate attributes without magnitude, permitting only frequency and classification. Ordinal scales sort categories without quantifiable order, allowing classification and sorting operations. Interval scales record attributes with fixed intervals but lack an absolute zero, allowing distance comparison but not ratios. Ratio scales assess attributes with interpretable and comparable ratios. Stevens also outlined permissible central tendency and dispersion measures for each scale: nominal scales report modes; ordinal scales report modes, medians, and percentiles but not means or standard deviations; interval or ratio scales allow means, standard deviations, and correlation coefficients. Ratio-based statistics, like the coefficient of variation, are only admissible for ratio scales.

1.1.1. Is it Permissible to do Factor Analysis with Ordinal Items?

Stevens’s taxonomy profoundly influenced the social sciences, particularly psychology, guiding statistical analysis for qualitative and quantitative variables. Many textbooks and statistical software rely on his taxonomy for recommendations on appropriate statistical analysis methods. However, criticisms of Stevens’s taxonomy have emerged (Lord, 1953; Velleman & Wilkinson, 1993). Controversy remains regarding the treatment of variables measured on an ordinal scale (Feuerstahler, 2023). Stevens recognizes that most scales used effectively by psychologists are ordinal scales and thus “In the strictest propriety the ordinary statistics involving means and standard deviations ought not to be used with these scales” (Stevens, 1946, p. 679). However, on the same page, Stevens recognizes that “On the other hand, for this ‘illegal’ statisticizing there can be invoked a kind of pragmatic sanction: In numerous instances, it leads to fruitful results.” Despite pragmatic justification for their use, some methodologists advocate against using such statistics for ordinal variables (Jamieson, 2004; Liddell & Kruschke, 2018). According to Stevens, only the mode, median, and percentiles are permissible for ordinal variables. Prohibiting statistics like mean and correlation coefficients raises concerns for psychometry, where Likert items and factor analyses are prevalent. If statistics not permissible for ordinal scales are banned, could it compromise the validity of “correlational psychology” (Spearman, 1904b, p. 205) and countless psychometric scales derived from ordinal measures through factor analysis on the Pearson correlation matrix? The validity of decades of psychometric research may be at risk?

2. Other Methods to Analyze Ordinal Items

Confirmatory factor analysis (CFA) emerged in the latter half of the 20th century, driven by the need to empirically validate theoretical models across diverse populations. Jöreskorg and Sörbom developed LisRel, the first commercial software for structural equation modeling, laying the groundwork for CFA. Babakus, Ferguson, & Jöreskog (1987) conducted a seminal simulation study on CFA of ordinal items, finding that factor loadings estimated with polychoric correlations surpassed those obtained with Pearson’s, Spearman’s, and Kendall’s correlations. This study confirmed previous observations that Pearson correlations tend to underestimate factor loadings, particularly with ordinal items having fewer than 5 categories and considerable asymmetry. Independent simulation studies have corroborated these findings, noting that as the number of categories increases and distribution bias decreases, Pearson correlation estimates approach those of quantitative variables (Bollen, 1989; Rhemtulla et al., 2012).

2.1. The Spearman Correlation Coefficient

Spearman developed his correlation coefficient as an alternative to Pearson’s, specifically for variables that couldn’t be quantitatively measured. Spearman’s rho, based on rank differences, is an adaptation of Pearson’s coefficient applied to the ranks (r_x1i, r_x2i ) rather than the original observations (x_1i, x_2i) (Spearman, 1904a):

$$R_S={\displaystyle \frac{{\displaystyle \sum _{i=1}^n(r_{x_{1i}}-{\overline{r}}_{x1})(r_{x_{2i}}-{\overline{r}}_{x_2})}}{\sqrt{{{\displaystyle \sum _{i=1}^n\left(r_{x_{1i}}-{\overline{r}}_{x_1}\right)}}^2{\displaystyle \sum _{i=1}^n{\left(r_{x_{2i}}-{\overline{r}}_{x_2}\right)}^2}}}}$$

(13)

2.2. The Polychoric Correlation Coefficient

The polychoric correlation coefficient, introduced by Karl Pearson in 1900, evaluates correlations between non-quantitatively measurable characteristics. This coefficient gained attention in factor analysis after simulation studies demonstrated its superiority in estimating associations between ordinal items (Babakus et al., 1987; Finney & DiStefano, 2013; Finney et al., 2016a; Holgado-Tello et al., 2009; Li, 2016; Muthén, 1984). It assumes that qualitative characteristics on ordinal or nominal scales result from discretized continuous latent variables with bivariate normal distributions. Thus, the polychoric correlation coefficient estimates the linear association between two latent variables operationalized by ordinal variables (e.g., Likert-type items). For instance, an item with 5 categories, such as an agreement scale, is conceptualized as dividing a latent variable into categories defined by thresholds on the latent z:

$$X=\{\begin{array}{l}1-\mathrm{Completely} \mathrm{disagree}\hspace{0.17em}\hspace{0.17em}if \xi \le {\xi }_1\hfill \\ 2-\mathrm{Disagree}\hspace{0.17em}\hspace{0.17em}if {\xi }_1<\xi \le {\xi }_2\hfill \\ 3-\mathrm{Neither} \mathrm{disagree} \mathrm{or} \mathrm{agree} \hspace{0.17em}if {\xi }_2<\xi \le {\xi }_3\hfill \\ 4-\mathrm{Agree} \hspace{0.17em}if {\xi }_3<\xi \le {\xi }_4\hfill \\ 5-\mathrm{Completely} \mathrm{agree}\hspace{0.17em}\hspace{0.17em}if \xi >{\tau }_4\hfill \end{array}$$

(14)

The observed variable’s five response categories represent an approximation of the five intervals at which the latent continuous variable, denoted as ξ, has been divided using four cutoff points or thresholds: ξ₁, ξ₂, ξ₃, and ξ₄. Each of these k-1 thresholds can be determined based on the relative frequency of the k categories, assuming normality:

$${\xi }_k={\Phi }^{-1}\left({\displaystyle \sum _{i=1}^k\frac{n_k}{N}}\right)$$

(15)

where F^-1 is the inverse of the normal distribution, n_k is the absolute frequency of the item’s k category, and N is the total of the item’s responses. By convention ξ₀=-∞ and ξ₅ =+∞. The polychoric correlation coefficient is generally estimated by maximum likelihood (see e.g. Drasgow, 1986). In a two-step algorithm, the procedure begins by estimating the tresholds for x₁ and for x₂. If the category i of variable x₁ is observed when x_1i-1 ≤x₁< x_1i and if category j of variable x₂ is observed when x_2j-1 ≤ x₂ < x_2j then the joint probability (P_ij ) of observing the category x_1i and x_2j is:

$$P_{ij}={\displaystyle \underset{{\xi }_{1i-1}}{\overset{{\xi }_{1i}}{\int }}{\displaystyle \underset{{\xi }_{2j-1}}{\overset{\xi 2j}{\int }}\varphi ({\xi }_1,{\xi }_2;\rho )\partial {\xi }_2\partial {\xi }_1}}$$

(16)

where

$$\varphi ({\xi }_1,{\xi }_2;\rho )={\displaystyle \frac{1}{2\pi \sqrt{(1-{\rho }^2)}}}\times e^{\left({\displaystyle \frac{-1}{2(1-{\rho }^2)}}({\xi }_1^2-2\rho {\xi }_1{\xi }_2+{\xi }_2^2\right)}$$

(17)

is the bivariate normal density function of x₁ and x₂ with m =0 and s = 1 (without loss of generality, since the original latent variables are origin- and scale-free), and Pearson correlation r. If n_ij is the number of observations x_1i of X₁ and x_2j of X₂, the sample likelihood is:

$$L=k{\displaystyle \prod _{i=1}^r{\displaystyle \prod _{j=1}^s{P}_{ij}^{n_{ij}}}}$$

(18)

where k is a constant and r and s are the number of categories of X₁ and X_2, respectively. In the second step, the maximum likelihood estimate of the polychoric correlation coefficient (r) is obtained by deriving the Ln(L) in respect to r:

$$\begin{array}{l}{\displaystyle \frac{\partial Ln(L)}{\partial \rho }}={\displaystyle \sum _{i=1}^r{\displaystyle \sum _{j=1}^s\frac{n_{ij}}{P_{ij}}}}\left[\varphi ({\xi }_{1i},{\xi }_{2j};\rho )-\varphi ({\xi }_{1i-1},{\xi }_{2j};\rho )\right]\hspace{0.17em}-\\ -\varphi ({\xi }_{1i},{\xi }_{2j-1};\rho )+\varphi ({\xi }_{1i-1},{\xi }_{2j};\rho )\end{array}$$

(19)

Equating the partial derivatives of (19), in respect to all model parameters, to 0 and solving it in respect to r gives the polychoric correlation estimate (Drasgow, 1986). The asymptotic covariance matrix of the estimated polychoric correlations can then be estimated using computational methods derived by Jöreskog (1994) after Olsson (1979). The maximum likelihood solution is, generally, obtained iteratively in several commercial software (e.g., Mplus, LisRel, EQS), free software ( R psych and polycor packages) but not in IBM SPSS Statistics (up to version 29).

Simulation studies suggest that the polychoric correlation coefficient outperforms Pearson, Spearman, and Kendall correlations in reproducing the correlational structure of latent variables when operationalized as ordinal observed variables However, its usefulness as an alternative to the Pearson correlation coefficient is debated. For small samples (n < 300 - 500) with many variables, polychoric correlations may not yield better results than Pearson correlations (Rigdon & Ferguson, 2006). Additionally, assumptions of bivariate normality and large sample sizes (> 2000 observations) are needed for accurate estimation, posing challenges in practical structural equation modeling (Arbuckle, 2006; Yung & Bentler, 1994). Moreover, the polychoric correlation lacks statistical robustness when the underlying distribution deviates strongly from bivariate normality, leading to potential inaccuracies in ordinal factor analysis (Ekström, 2011; Robitzsch, 2020). It may also introduce bias when ordinal factor models impose a normality assumption for underlying continuous variables, which might also often be true in empirical applications (Grønneberg & Foldnes, 2024; Robitzsch, 2020, 2022).

This paper compares Pearson’s, Spearman’s, and Polychoric correlation coefficients using real and simulated data. It demonstrates calculating these coefficients using statistical software and applying the matrices in exploratory and confirmatory factor analysis (EFA and CFA). The findings are discussed regarding their implications for practical factor analysis applications.

3. Results

3.1. Case Study 1. Exploratory Factor Analysis of Ordinal Items with 3 Points

To demonstrate factor analysis with ordinal items, data from a study on voters’ perceptions (Marôco, 2021a) using a 3-point ordinal scale is considered. Items include “People like me have no voice” (WithoutVoice), “Politics and government are too complex” (Complex), “Politicians are not interested in people like me” (Desinter), “Politicians once elected quickly lose contact with their voters” (Contact), and “Political parties are only interested in people’s votes” (Parties). Using IBM SPSS Statistics, classical EFA can be conducted on the covariance matrix or Pearson correlation matrix. Syntax for factor extraction using principal components with Varimax rotation on the correlation matrix is provided as supplementary material. Furthermore, EFA can be performed on the Spearman correlation matrix in SPSS Statistics. Syntax involves converting variables to ranks before conducting factor analysis. Lastly, EFA on the polychoric correlations matrix is discussed. Since SPSS Statistics does not compute polychoric correlations, one can calculate and import the matrix from R using the polycor package (Fox, 2016). Pearson, Spearman, and Polychoric correlation matrices for the 5 items are shown in Table 1.

With IBM SPSS Statistics, exploratory factor analysis (EFA) can now be conducted on the polychoric correlation matrix, along with other desired correlation matrices. Table 2 provides a summary of the relevant results from the three EFAs conducted on the Pearson, Spearman, and Polychoric correlation matrices. Factors were extracted using Kaiser’s criterion (eigenvalue greater than 1).

The variance extracted by each of the two retained factors is higher in the analysis conducted on the polychoric correlation matrix compared to the analyses on the Pearson or Spearman correlation matrices (59.9% vs. 53.3% vs. 51.7%, respectively, for the first factor F1; the differences for the second factor F2 are even smaller). Similarly, the factor loadings for each factor are higher when using the polychoric correlations, followed by Pearson’s and then Spearman’s correlations. These differences, often to the second decimal place, are visually depicted in Figure 2, highlighting their significance.

3.2. Case Study 2. Confirmatory Factor Analysis of Ordinal Items with 5 Points

Case study two involves confirmatory factor analysis (CFA) of the USEI, a psychometric inventory assessing academic engagement. Comprising 15 ordinal items, the USEI’s response categories range from “1-Never” to “5-Always” (Maroco, Maroco, Bonini Campos, & Fredricks, 2016). It conceptualizes academic involvement as a second-order hierarchical structure, with three first-order factors: Behavioral, Emotional, and Cognitive engagement, each consisting of five ordinal items (Maroco et al., 2016). Although most items exhibit slight leftward skewness, observed skewness (sk) and kurtosis (ku) values do not indicate severe departures from normal distribution (Marôco, 2021b).

Using the R software and the lavaan package (Rosseel, 2012), CFA can be conducted employing various correlation matrices and appropriate estimators based on the measurement scales of observed variables. The recommended estimators include maximum likelihood for normally distributed, at least interval-scale items, or WLSMV/DWLS for ordinal items with underlying bivariate normal latents (Finney et al., 2016). Traditional CFA, utilizing the covariance matrix and maximum likelihood estimation of the model described in equation (11), can be executed with the lavaan library in R using the Pearson correlation matrix (see supplementary material for R code). To perform CFA on the polychoric correlation matrix, ordinal item specification with the ordered keyword is necessary (by default, lavaan utilizes the DWLS estimator). Table 3 summarizes the primary outcomes of these two analyses. Both models exhibit good fit, whether using the traditional method (maximum likelihood on the covariance matrix) or the more recent method (DWLS on the polychoric correlation matrix). All items display standardized factor loadings exceeding the typical cutoff value in factor analysis (λ ≥ 0.5). Differences in factor loading estimates are minimal, mostly favoring the DWLS method on the polychoric correlation matrix.

3.3. Case Study 3. Confirmatory Factor Analysis of Ordinal Items - Artificial Results?

The case studies thus far suggest an advantage in conducting either exploratory or confirmatory factor analysis of ordinal items using the polychoric correlation matrix over Pearson’s correlation. However, it’s crucial to note that the observed differences in standardized factor loadings between the two methods are typically minimal, only affecting the second decimal place of estimates, and do not significantly alter the researcher’s conclusions regarding the underlying factor structure. Additionally, as highlighted by Finney et al. (2016), while the results with ML and DWLS estimators are promising, there remains uncertainty about their performance with ordinal items exhibiting varying numbers of categories and distributions.

To further illustrate this uncertainty, a case study of the Postpartum Bonding Questionnaire (PBQ) with 25 ordinal items and 6 categories each is presented (Saur et al., 2016). Descriptive statistics for the PBQ items are provided (see supplementary material). Employing the lavaan package, CFA was conducted on both the Pearson correlation matrix with ML and the Polychoric correlation matrix with DWLS estimation. The results, summarized in Table 4, indicate notable differences favoring the analysis of the polychoric correlation matrix with DWLS.

In this case study, contrary to the previous ones, the differences in standardized factor loadings, now extending to the first decimal place, significantly favor the analysis of the polychoric correlation matrix with DWLS. While ML analysis on the covariance matrix yields standardized weights mostly below 0.5, DWLS analysis on the polychoric correlation matrix results in higher factor loadings. Consequently, the goodness of fit for the tetra-factorial model obtained with DWLS on the polychoric correlation matrix is substantially better than that obtained with ML on the Pearson covariance matrix. Specifically, using DWLS on the polychoric correlation matrix, the model fit is deemed good, whereas ML on the covariance matrix yields a relatively poor fit.

The apparent superiority of CFA estimates derived from polychoric correlations compared to those from Pearson’s covariance/correlation matrix is evident. However, the reliability of these polychoric correlation-based estimates relative to those from the Pearson correlation matrix warrants scrutiny. Particularly in Case Study 3, where items exhibit significant asymmetry, it raises questions about the accuracy of these estimates.

To address this, a simulation study using various levels of skewness and kurtosis for 5-point ordinal variables with known factor structures can provide valuable insights. By comparing the performance of CFA estimates derived from polychoric correlations with those from Pearson’s correlation matrix under different data conditions, we can evaluate the robustness and reliability of the former approach.

4. Simulation Studies

Let’s consider a two-factor population model with three items each (see supplementary material), where the standardized factor loadings are 0.5, 0.6, and 0.7 for each item, and a correlation of 0.5 between the two factors. Using the lavaan package, we can generate one million standardized data observations with skewness and kurtosis values of 0. The descriptive statistics and correlations of the generated items are shown in Figure 3A, indicating normal distribution with small to moderate observed Pearson correlations between the items.

Assuming the ordinal items result from discretizing latent continuous variables into a given number of categories (ncat = 5, 7, and 9), the continuous variables generated previously can be discretized accordingly. The categorized data is depicted in Figure 3B.

Subsequent simulations were conducted under two scenarios: (i) Symmetrical (normal) ordinal variables with 5, 7, and 9 categories; (ii) Biased (non-symmetrical) ordinal data generated from biased continuous variables; and (iii) Biased (non-symmetrical) variables obtained by biased sampling from ordinal symmetrical variables.

4.1. Simulation 1. Symmetrical Ordinal Data Generated from Normal Continuous Data

A model was created with R-lavaan to fit the two-factor population data. It was then applied using DWLS or ML to data sampled with replacement from the original ordinal generated data with sample sizes ranging from 30 to 500 and 500 bootstrap replicates. Results in Figure 4 show that ML on Pearson correlations for ordinal items consistently underestimated true population factor loadings. This underestimation increased with fewer categories (from 9 to 5). Conversely, DWLS estimation quickly converged to true population factor loadings, unaffected by the number of categories, even for small sample sizes.

4.2. Simulation 2. Asymmetrical Ordinal Data Generated from Non-Normal Continuous Data

Simulation two mirrored simulation one, except the original data was skewed and leptokurtic (sk=-3, ku=7). Results from 500 bootstrap replicates are shown in Figure 5. ML and DWLS yielded biased estimates, with ML notably poor for 5-category ordinal variables. DWLS showed resilience to category count, converging even with small sample sizes in SEM (n=100-150), albeit with noticeable bias.

Figure 5. Median bias (population – sample) estimates for factor loadings and nonparametric 95% confidence intervals (shaded) derived from 500 bootstrap samples of the original skewed ordinal data with varying category counts (ncat: 5, 7, or 9) and sample sizes (30 to 500).

Figure 5. Median bias (population – sample) estimates for factor loadings and nonparametric 95% confidence intervals (shaded) derived from 500 bootstrap samples of the original skewed ordinal data with varying category counts (ncat: 5, 7, or 9) and sample sizes (30 to 500).

4.3. Simulation 3. Asymmetrical Ordinal Data Generated from Biased Sampling of Normal Ordinal Data

Simulation three explores biased responses stemming from non-random sampling practices prevalent in psychological and social sciences research. Despite a normative population exhibiting a normal distribution, non-random sampling and respondent bias lead to highly asymmetrical data. Only biased responses from the original population (N=10⁶) were sampled. As shown in Figure 6, with data from biased items featuring 7 categories, ML consistently underestimates the true population factor loadings. Surprisingly, DWLS exhibits an opposite trend, overestimating the true population factor loadings. This contrasts with the asymmetrical data generated in simulation 2.

5. Discussion

The simulation and case studies demonstrate that polychoric correlation coefficients tend to be higher than Pearson and Spearman coefficients, indicating their robustness in estimating latent variables. Polychoric correlations are derived from the normal distribution of latent variables, making them less susceptible to measurement errors. Consequently, factor loadings extracted from the polychoric correlation matrix are higher, supporting consistent theoretical interpretations, even with small sample sizes. Despite this, relying solely on polychoric correlations and weighted least squares methods for factor analysis with ordinal variables and few categories may not be justifiable.

In case study 1 and 2, both Pearson’s ML estimation and polychoric DWLS estimation produced similar results, with negligible differences in standardized factor loadings, supporting consistent interpretations of factor validity. However, case study 3 revealed substantial disparities between estimates obtained from Pearson correlations with ML estimation and polychoric correlations with DWLS estimation. The latter yielded substantially higher factor loadings, leading to conflicting conclusions regarding factorial validity.

Moreover, simulation study 3 highlighted the unreliability of polychoric correlations with DWLS under biased data conditions. ML estimation on Pearson correlations consistently underestimated true population parameters, while DWLS on polychoric correlations overestimated them, particularly evident in biased data scenarios.

The inadequacy of least squares weighted methods with polychoric correlation matrices has also been reported, both under conditions of moderate (Li, 2015; 2016) and severe non-normality (Foldnes & Grønneberg, 2020). Recommendations for method selection based on category counts exist, but our findings suggest that methodological differences are minor under sensible data conditions but significant under strongly biased scenarios. Blind reliance on polychoric correlations with weighted least squares estimates may yield misleading conclusions about factor structure validity and reliability.

6. Concluding Remarks

Method selection should consider data characteristics, with sensitivity analyses imperative, especially under biased conditions. Blind reliance on polychoric correlations with weighted least squares estimates may yield misleading conclusions about factor structure validity and reliability. Further research and sensitivity analyses are essential to ensure robust psychometric evaluations. These findings contribute to ongoing discussions about the appropriate use of correlation matrices and estimation methods in factor analysis.