Circumplex Models With Ordinal Data

Abstract Circumplex models arrange personality or mood variables around a circle such that the relative locations of these variables reflect their correlations. Although circumplex models were originally developed for continuous variables, their applications often involve ordinal variables. To accommodate ordinal variables, we propose estimating Browne’s circumplex correlation structure with polychoric correlations instead of Pearson correlations. We consider ordinary least squares estimation due to its computational robustness with polychoric correlation matrices. We illustrate the newly developed method with an empirical study that involves 18 binary variables and 12,108 participants. We also conduct a simulation study to explore its statistical property and compare it with maximum likelihood estimation with polychoric correlations, ordinary least squares estimation with Pearson correlations, and maximum likelihood estimation with Pearson correlations. The new method provides essentially unbiased point estimates and more satisfactory confidence intervals than other methods.


Introduction
The circumplex model is an alternative to the common factor model for understanding correlations among multiple measurements (Guttman, 1954). The common factor model allows researchers to explain such correlations using a smaller number of factors and to interpret the factors using a factor loading matrix. The circumplex model allows researchers to explain such correlations using a circular representation of those measurements. The larger the correlation between two measurements, the closer their representations are on the circle.
The circumplex model has been instrumental in the development of several psychological theories. Watson et al. (1999) described a circumplex model for human affects. Wiggins (2003) described a circumplex model for interpersonal behaviors. Hofstee et al. (1992) integrated circumplex models with the Big Five model. Locke (2019) described a circumplex model for the interpersonal culture in work teams and organizations. Fabrigar et al. (1997) surveyed applications of circumplex models in social and personality psychology, and they found that the most popular method of analysis is explanatory factor analysis (principal component analysis) and multidimensional scaling. Several researchers contributed to direct specifications of the circumplex model as a covariance/correlation structure (Anderson, 1960;Browne, 1992;Cudeck, 1986;Guttman, 1954). In particular, Browne (1992)  All currently available methods for estimating circumplex models require continuous data, but many applications of circumplex models involve ordinal variables. Participants need to endorse one response from several alternatives (e.g., "disagree," "neutral," and "agree"). A common practice is to estimate circumplex models with ordinal data as if they were normal, which is convenient but not statistically justified. In the context of factor analysis, treating ordinal data as if they were continuous can lead to biased point estimates, standard error estimates, and test statistics (Olsson, 1979b;Rhemtulla et al., 2012).
We suggest estimating circumplex models with polychoric correlations (Olsson, 1979a) to accommodate ordinal data. Such an accommodation is a common practice in structural equation models with ordinal data (Flora & Curran, 2004;Forero et al., 2009;Muth en, 1984;Rhemtulla et al., 2012), but it has never been used for circumplex models. The development is nontrivial because the corresponding standard errors and test statistics involve deriving derivatives of discrepancy functions with respect to model parameters and polychoric correlations.
The rest of the article is organized as follows. We first describe how to use polychoric correlations to accommodate ordinal variables. We then adapt Browne's circumplex model (1992) such that we can estimate it with polychoric correlations. We next demonstrate the proposed method with an empirical data set. We also assess its statistical properties with simulated data. We conclude the paper with several additional comments.

Ordinal Variables and Polychoric Correlations
A popular method to accommodate structural equation models with ordinal variables is the latent response variable approach, which assumes that a continuous variable x Ã underlies an ordinal variable x. The correspondence between a c-category variable x and its latent response variable x Ã is Here s 0 ¼ À1, s 1 <s 2 < . . . <s cÀ1 , s c ¼ 1 are threshold values and "()" denotes "is equivalent to." Although we do not directly observe such latent response variables, we can estimate their correlations from the observed ordinal variables. These correlations are called polychoric correlations (tetrachoric correlations for binary variables). We estimate polychoric correlations using Olsson's (1979b) two-stage method. This method first estimates threshold values for each of the p ordinal variables. It then estimates a polychoric correlation for each pair of variables with the threshold values fixed as those obtained at the first stage. Thus, a p Â p polychoric correlation matrix contains pðpÀ1Þ=2 unique polychoric correlations.
Estimating a structural equation model with polychoric correlations provides point estimates, but standard errors and test statistics involve the pðpÀ1Þ=2 Â pðpÀ1Þ=2 asymptotic covariance matrix (ACM) of the polychoric correlations. Muth en (1984) derived the ACM using a Hessian matrix of log-likelihood functions, and J€ oreskog (1994) derived the ACM using four-way contingency tables of ordinal variables. Christoffersson and Gunsj€ o (1996) presented a conceptually simpler derivation of the ACM. The derivation uses estimation equations that involve polychoric correlations, thresholds, and two-way contingency tables of ordinal variables. However, they did not include details like the derivatives of estimation equations with respect to polychoric correlations and thresholds, the derivatives of estimation equations with respect to two-way contingency tables. Zhang et al. (2021) recently presented such details and also described how to improve computing efficiency using parallel computing.
Although these three methods are asymptotically equivalent, their performance may not be satisfactory for small to medium-sized samples. Monroe (2018) suggests a Monte Carlo method to estimate the ACM in these situations. The Monte Carlo ACM differs from the asymptotic methods in how to estimate the four-way contingency tables of ordinal variables. It generates a very large sample with the parameter values chosen to be the sample estimates. Therefore, it avoids many issues associated with estimating four-way contingency tables with small to medium-sized samples. For example, many entries in the four-way contingency tables are empty when we construct them using sample data, but the number of empty entries is much smaller with large simulated data. The Monte Carlo ACMs are most useful for small to medium-sized samples because four-way contingency tables made directly from these samples are often poor estimates of their population counterparts. On the other hand, the four-way contingency tables are accurate estimates in large samples, which render Monte Carlo ACMs unnecessary.

Model Specification and Model Estimation
Browne (1992) derived a circumplex model for continuous variables and estimated the model with Pearson correlations. We now adapt his approach for ordinal variables. We propose estimating the circumplex correlation structure with polychoric correlations instead of Pearson correlations.

A Circumplex Correlation Structure
Let P x Ã be a p Â p polychoric correlation matrix of latent continuous response variables x Ã 1 , x Ã 2 , . . . , x Ã p : Note that Equation (1) specifies the correspondence between these latent response variables and observed ordinal variables x 1 , x 2 , . . . , x p : We now specify the circumplex correlation structure (Browne, 1992, Eq. (3)) for the polychoric correlation matrix Here D Ã f is a diagonal scaling matrix, P c is the common score correlation matrix, and D m is a diagonal matrix of measurement error variances. Let f Ã ii be the ith diagonal element D Ã f and ii be the ith diagonal element of D m : The element f Ã ii is the correlation between the continuous latent response variable x Ã i and its common part c i , and we refer to it as the communality index for the ith variable. The communality index f Ã ii is a function of the corresponding measurement error variance ii Let qðc i , c j Þ be a typical element of the common score correlation matrix P c and c i and c j be the common parts of x Ã i and x Ã j , respectively. Browne (1992), Eq. (34) expressed such common score correlations as weighted sums of cosine functions.
The parameter h i is the angle between c i and a reference variable in a circular placement. The angle parameter is between 0 degrees and 360 degrees. Note that both an angle of 1 degree and an angle of 359 degrees indicate that c i is close to the reference variable. An angle of 180 degrees indicates that c i is the farthest from the reference variable. Equation (4) involves m cosine functions. One can improve the model fit by considering a larger m, but it also increases the model complexity at the same time. Choosing an appropriate m for a circumplex model is similar to choosing an appropriate number of factors in factor analysis. The goal is to find a balance between a good model fit and model parsimony.
The coefficients b 1 , b 2 , Á Á Á , b m are all nonnegative and indicate the relative contributions of the m cosine functions. The intercept b 0 is also nonnegative and is the simple function of other coefficients b 0 ¼ 1Àðb 1 þ Á Á Á þ b m Þ: These coefficients determine the minimal common score correlation Here ½m=2 is the largest integer that is less than m=2: To summarize, the model involves p À 1 angle parameters, m coefficients of cosine functions, and p communality indices (measurement error variances).

Point Estimates
Let c be a vector that contains all the model parameters. When manifest variables are normally distributed, one can estimate c by minimizing the maximum likelihood (ML) discrepancy function Here R is a Pearson correlation matrix of normal manifest variables and PðcÞ is the correlation structure implied by Equations (2) and (4).
When the manifest variables are ordinal, one may consider replacing the Pearson correlation matrix with a polychoric correlation matrix in Equation (6). This method will have two issues. First, polychoric and Pearson correlations have different distribution properties. Second, Equation (6) requires R to be positive definite, but polychoric correlation matrices are often not. We propose estimating c by minimizing the ordinary least squares (OLS) discrepancy function f OLS ¼ ðrÀqðcÞÞ 0 ðrÀqðcÞÞ: Here r is a column vector that collects nonduplicated elements of a polychoric correlation R. Similarly, qðcÞ is the corresponding vector for the model implied correlation matrix PðcÞ 1 .

Standard Error and a Test Statistic
Letĉ be the vector of parameter estimates obtained by minimizing the OLS discrepancy function of Equation (7). We can compute its asymptotic covariance matrix (acov) using a sandwich method: Here A is a matrix of second-order derivatives o 2 f OLS ococ 0 : The matrix B is Here D is a pðpÀ1Þ=2 by 2pÀ1 þ m matrix: The elements of D are derivatives of the circumplex correlation structures with respect to the model parameters. At the middle of Equation (9) is the ACM of the ploychoric correlations, Appendices A and B contain the technical details of the matrices D and A, respectively.
The large sample covariance matrix forĉ of Equation (8) is similar to Proposition 2 in Browne (1984), but our result is specifically for the circumplex model with ordinal data. The square roots of diagonal elements of the large sample covariance matrix provide standard error estimates forĉ: We can use them to construct confidence intervals for c: We adapt Proposition 4 in Browne (1984) to assess how a circumplex model fits the ordinal data. Let e be a vector of residuals e ¼ rÀqðcÞ: The test statistic is Here D c is a matrix of the null space of D such that D 0 D c ¼ 0: Let p Ã be the number of nonduplicated elements of R and let q be the number of parameters of the circumplex model. The order of D c is p Ã Â ðp Ã ÀqÞ: The large sample distribution of the test statistic is a central v 2 distribution with p Ã Àq degrees of freedom if the model perfectly fits the population correlation matrix. When the model closely fits the population correlation matrix, the distribution becomes a noncentral v 2 distribution with a noncentral parameter Tðc, q).
Both the sandwich standard error estimator and the test statistic involve the ACM of the polychoric correlations. We estimate the ACM using two methods. The first method involves estimation equations (Christoffersson & Gunsj€ o, 1996;Zhang et al., 2021) and we call the regular ACM the "sample-based" ACM because the estimation uses only sample data. In contrast, the second method (Monroe, 2018) involves simulating an additional large sample, and we call the resulting ACM a "Monte Carlo" ACM. Pairing these two types of ACMs with the standard error estimator of Equation (8) and the test statistic of Equation (12) will produce two types of standard error estimates and test statistics.

An Empirical Illustration
We now illustrate the circumplex model with ordinal data using a real data set. Trampe et al. (2015) reported a study on the emotion structure of 12,108 participants. In this study, participants rated themselves on 18 emotions (nine positive emotions: "pride," "love," "hope," "gratitude," "joy," "satisfaction," "awe," "amusement," and "alertness"; nine negative emotions: "anxiety," "disdain," "offense," "guilt," 1 Adding a weight matrix to the center of Equation (7) changes it to a weighted least squares discrepancy function. Although weighted least squares estimation is theoretically superior to OLS estimation, several researchers (Forero et al., 2009;Yang-Wallentin et al., 2010) reported that OLS estimation performs more satisfactorily in practice for factor analysis with ordinal variables. "disgust," "fear," "embarrassment," "sadness," and "anger"). The rating scales were binary (0 ¼ not experiencing a certain emotion and 1 ¼ experiencing a certain emotion). Figure 1 displays these 18 variables. The proportion of experiencing is lower than the proportion of not experiencing for each emotion. Additionally, participants tend to experience positive emotions more frequently than negative emotions. The frequency of positive emotions ranged from 4.49% ("awe") to 32.63% ("joy"), with a median of 21.90% ("hope"). The frequency of negative emotions ranged from 1.06% ("disdain") to 29.55% ("anxiety"), with a median of 5.67% ("guilty").
We first estimate the tetrachoric correlations of the 18 binary variables and their ACM. We then estimate the circumplex model with the tetrachoric correlation matrix using OLS. We refer to the corresponding OLS estimates as "PC OLS-A " (with sample-based standard errors) and "PC OLS-B " (with Monte Carlo method standard errors), separately. In addition, we conducted three more analyses for comparison purposes. They are ML with polychoric correlations, OLS with Pearson correlations, and ML with Pearson correlations. We refer to these three combinations of correlation types and estimation methods as "PC ML ," "PS OLS ," and "PS ML ," respectively. We compute the tetrachoric correlation matrix and its ACM using "Turbofuns" (Zhang et al., 2021) and implement the OLS estimation method using R 2 . We obtain ML estimates using the R package "CircE" (Grassi et al., 2010). Appendix C contains the corresponding R code. Table 1 reports point estimates and their standard error estimates for all 36 parameters. 3 They include 17 angle parameters (the angle of the reference emotion "pride" is always zero), one coefficient of a cosine function, and 18 communality indices. To facilitate the comparison among the four combinations of estimation-correlation, we display their results in Figure 2.
Note that the point estimates of PC OLS-A and PC OLS-B are always the same, but their standard error estimates can be different, particularly for small and medium-sized samples. The standard error estimates of PC OLS-A and PC OLS-B are very close in the illustration example due to the large sample size. Therefore, we only present the results of PC OLS-A in Table 1.
Both the type of correlation and the estimation method affect the point estimates, but their influences are different. The same estimation method produced similar estimates for angle parameters when combined with different types of correlation. The estimates for the angle parameters of PC OLS-A and PS OLS are similar. The estimates for the angle parameters of PC ML and PS ML are also similar. The point estimates for the communality indices show a completely different pattern. The type of correlation had more substantial influences on them than the estimation method. The estimates for the communality indices of PC OLS-A and PC ML are similar, and the corresponding estimates of PS OLS and PS ML are also similar. The estimates for the communality indices were much higher when the two estimation methods combined with the polychoric correlations. It is interesting to note that the standard error estimates of the PC OLS-A approach were the largest of the four approaches. Standard errors of the other three methods require normally distributed manifest variables, which is clearly violated for highly skewed binary data.
The v 2 statistics for the circumplex model are 1071.51, 12658.43, 1671.32, and 1670.00, for PC OLS-A 4 , PC ML , PS OLS , and PS ML , respectively. We compare them against the The development is available as a user-friendly R package "CircumO." 3 We also constructed confidence intervals. We include the confidence intervals in the online support file.  reference v 2 distribution with 117 degrees of freedom. We thus reject the test of perfect fit with all four statistics. We computed the root mean square error of approximation (RMSEA) (Browne & Cudeck, 1993) to further assess the model fit. They are 0.026, 0.094, 0.033, and 0.033, for PC OLS-A , PC ML , PS OLS , and PS ML, respectively. According to a guideline suggested by Browne and Cudeck (1993), the RMSEA value of 0.05 indicates that an SEM model closely fits the data. Thus, the only PC ML produced a poor model fit.
We need to be cautious in generalizing the results of the empirical illustration. In particular, the results were obtained with highly skewed binary data. Next, we conduct a simulation study to assess the statistical properties of the methods under a variety of conditions.

A Simulation Study
We manipulated four factors in the simulation study: the population model (Model I and Model II), the number of categories of ordinal variables (two-category and three-category), the shape of ordinal variables (symmetric distributions and four types of asymmetric distributions), and the sample size (100, 200, 500, 1000, and 12,108). Crossing these four factors results in 100 conditions. We simulate 1,000 samples under each condition.

Population Parameters and Data Generation
The simulation study includes two population models, which represent an ideal case and a more realistic case, respectively. Model I is a theoretical illustration of the circumplex model where eight variables are evenly distributed around a circle (Russell, 1980). Figure 3 displays the model graphically. 5 Model II is the same as the one considered in the empirical example. The upper left panel of Figure 2 displays the model, and Table 1 presents the corresponding parameter values.
We consider two-category ordinal variables and threecategory ordinal variables because they deviate from normal distributions more substantially than ordinal variables with more categories. Therefore, accommodation of ordinal variables is most needed in these situations. Note that the data set in the empirical illustration involves only the binary variables. The expansion to both two-category and three-category variables will make the results of the simulation study relevant in more situations. We choose four sample size conditions (N ¼ 100, 200, 500, and 1000) to reflect the typical sample size in practice. The last sample size condition (N ¼ 12,108) corresponds to the sample size in the empirical example. Such a large sample size is relevant to internetbased research and offers an opportunity to assess asymptotic properties. Figure 4 displays the five shapes of distributions for the ordinal variables. We manipulate the distribution shapes for two reasons. First, some ordinal variables in the empirical study are highly skewed, and we want to assess how these methods perform when the ordinal variables are of a variety of shapes. Second, many studies show that the distribution of ordinal variables affects different statistical methods in the context of structural equation modeling (SEM) and factor analysis (Olsson, 1979b;Parry & McArdle, 1991;Garrido et al., 2013). We choose the five distributions of the ordinal variables as similar to those considered in Parry and McArdle (1991) and Garrido et al. (2013). We denote them using "S0," "S1," "MS1," "S2," and "MS2." Under the "S0" condition, all ordinal variables are evenly distributed in different categories, resulting in zero skewness. Under the "S1" condition, all ordinal variables are positively skewed. Under the "MS1" condition, half of the variables are positively skewed and the other half are negatively skewed. However, the skewness levels of "S1" and "MS1" are the same. Under the "S2" condition, all ordinal variables are positively skewed, but their skewness level is higher than that of "S1." Under the "MS2" condition, half of the variables are positively skewed and the other half are negatively skewed. The skewness levels of "S2" and "MS2" are the same. Note that the letter "M" in "MS1" and "MS2" indicates that the skewed variables are a mixture of different directions.
We generate continuous latent response variables from a multivariate normal distribution MVNð0, PðcÞ). The population covariance matrix PðcÞ is the model implied correlation matrix of Model I and Model II. We then convert the continuous latent response variables to observed ordinal variables according to Figure 4. For example, the middle right panel of Figure 4 displays "MS1" condition for threecategory variables. The odd-numbered variables (x 1 , x 3 , x 5 , . . . , x 17 ) are positively skewed and the proportions for "1," "2," and "3" are 50%, 30%, and 20%, respectively. These proportions correspond to the threshold s 1 ¼ 0.000 and s 2 ¼ 0.842. Therefore, latent response variable scores À0.4, 0.4, and 1.4 will correspond to the observed scores x ¼ 1, x ¼ 2, and x ¼ 3, respectively. The even-numbered variables (x 2 , x 4 , x 6 , . . . , x 18 ) are negatively skewed, and the three proportions are 20%, 30%, and 50%, respectively. These proportions correspond to the thresholds s 1 ¼ À0:842 and s 2 ¼ 0:000: The same latent response variable scores, À0.4, 0.4, and 1.4, will correspond to the observed scores x ¼ 2, x ¼ 3, and x ¼ 3, respectively. Therefore, we convert the continuous response variable scores x Ã 1 ¼ 0.1 and x Ã 2 ¼ 0.2 to the ordinal variable scores x 1 ¼ 2 and x 2 ¼ 3. We can use a similar method to convert continuous response variables to ordinal variables under other conditions.

Model Estimation
We estimate the circumplex models using both OLS and ML with both polychoric correlations and Pearson correlations. The corresponding standard errors and test statistics require ACM estimates for polychoric correlations.
We use the same notations ("PC OLS-A ," "PC OLS-B ," "PC ML, " "PS OLS ," and "PS ML ") introduced in the empirical illustration to denote these five methods. The two methods PC OLS-A and PC OLS-B have the same point estimates but different standard error estimates and test statistics.

Clean Convergence Rates
We first examine the computational robustness of the five methods. We regard an estimation as successful if it produces a clean convergence. In a clean convergence, none of its parameter estimates are on the boundary. Therefore, we exclude convergences with boundary estimates. Table 2 reports the clean convergence rates of four estimation Figure 4. The proportion of each category per condition. "S0," "S1," "MS1," "S2," and "MS2" stand for five types of distributions for ordinal variables. The numbers above the bars represent the percentage of each category. methods. 6 A quick comparison shows that PS OLS , PC OLS-A , PS ML , and PC ML have the highest, second highest, third highest, and lowest successful convergence rates under most conditions, respectively. We also investigated why some methods failed to produce clean convergences. The failures of the ML estimation were mainly due to nonconvergence, but those of the OLS estimation were mainly due to boundary estimates (zero error variance estimates). 7 The clean convergence rates of Model II are mostly higher than their counterparts in Model I, with one exception. The exception is PC ML , and it performed worse in Model II. As expected, the clean convergence rates are lower for smaller samples, more skewed distributions, and two-category ordinal variables. The sample sizes of 100 and 200 seem too small for Model I, but they seem large enough for Model II when the data distributions do not substantially deviate from the normal distribution except for the estimation of PC ML . In addition, the advantages of PC OLS-A over PS ML are most substantial in Model II with the three distributions "S0," "S1," and "MS1." On the other hand, when the data distributions substantially deviate from the normal distributions ("S2" and "MS2") and the sample size is small, PS ML has higher convergence rates than PC OLS-A . Therefore, we should be cautious when interpreting the results from small samples. Table 2 presents clean convergence rates where we successfully obtained point estimates, but we also need standard error estimates and test statistics to make statistical inferences about the circumplex model. As described earlier (Equation (8)), standard errors and test statistics of PC OLS-A and PC OLS-B require ACM estimates. In particular, computing test statistics involves inverting a large matrix (D 0 cĈD c in Equation (12)). Matrix inversion is more challenging for PC OLS-A with larger models and small samples because computing sample-based ACMs will be more difficult in these situations. For example, the matrix inversion of PC OLS-A failed in a large proportion of the simulated samples for Model II at N ¼ 100. We only include conditions for further analysis if we can successfully obtain point estimates, standard error estimates, and test statistics in more than 50% of the simulated samples.
samples. PS OLS and PS ML produced essentially unbiased estimates for b 1 under most conditions. The only exception is that their point estimates are negatively biased when all variables are substantially skewed in the same direction ("S2"). We observed similar patterns with the other parameters in Model I. Next, we examine whether the point estimates are biased in Model II. Similarly, we select an angle parameter (h 16 ), a coefficient of a cosine function (b 1 ), and a communality index (f 16 ). Figure 6 compares their point estimates with their population values. The results of Model I are largely replicated in Model II. In particular, PC OLS-A and PC ML produced essentially unbiased estimates when the sample size is 500 or larger. The point estimates of PS OLS and PS ML are even more biased than those of Model I. We observed similar patterns with the other parameters in Model II.

Confidence Intervals for Parameters
We now examine the empirical coverage rates of 95% confidence intervals of the parameters. The empirical coverage rate of a parameter is the proportion of simulated samples in which the confidence interval covers its population value. Confidence intervals combine point estimates and standard error estimates. We consider two ways of estimating standard errors: PC OLS-A (with sample-based ACMs) and PC OLS-B Figure 5. Means of parameter estimates (h 5 , b 1 , and f 5 ) in Model I. Rows correspond to different sample sizes and columns correspond to different parameters. The horizontal line in each plot represents the population value. "PC OLS-A -2C," "PC OLS-A -3C,""PC ML -2C," "PC ML -3C," "PS OLS -2C," "PS OLS -3C," "PS ML -2C," and "PS ML -3C" represent combinations of an estimation method (PC OLS-A , PC ML , PS OLS , and PS ML ) and the number of categories ("2C": two categories; "3C": three categories). "S0," "S1," "MS1," "S2," and "MS2" stand for five types of distributions for ordinal variables; "h": an angle parameter; "b": a cosine function coefficient parameter; and "f": a communality index parameter.
(with Monte Carlo ACMs). For comparison, we also examined the empirical coverage rates of three other types of confidence intervals (PC ML , PS OLS , and PS ML ). Table 3 presents the empirical coverage rates in Model I. Each entry is a mean empirical coverage rate computed across all parameters (seven angle parameters, one cosine function parameter, and eight communality index parameters) in the model. The mean coverage rates for PC OLS-A and PC OLS-B are close to 95% under all conditions where these two types of confidence intervals are assessed. Under some conditions, we encountered difficulties in constructing these two types of confidence intervals in many (larger than 50%) simulated samples. These difficult conditions tend to be combinations of small samples, binary data, and skewed distributions. In contrast, the mean empirical coverage rates of PC ML , PS OLS , and PS ML are lower than 95% under most conditions. In particular, the mean empirical coverage rates of PS OLS and PS ML become even lower for larger samples. 9 Table 4 presents the empirical coverage rates in Model II. Each entry is a mean empirical coverage rate computed Figure 6. Means of parameter estimates (h 16 , b 1 , and f 16 ) in Model II. Rows correspond to different sample sizes and columns correspond to different parameters. The horizontal line in each plot represents the population value. "PC OLS-A -2C," "PC OLS-A -3C," "PC ML -2C," "PC ML -3C," "PS OLS -2C," "PS OLS -3C," "PS ML -2C," and "PS ML -3C" represent combinations of an estimation method (PC OLS-A , PC ML , PS OLS , and PS ML ) and the number of categories ("2C": two categories; "3C": three categories). "S0," "S1," "MS1," "S2," and "MS2" stand for five types of distributions for ordinal variables; "h": an angle parameter; "b": a cosine function coefficient parameter; and "f": a communality index parameter. 9 The unsatisfactory empirical coverage rates of PS OLS and PS ML are mainly due to low empirical coverage rates for the communality indices. The unsatisfactory empirical coverage rates of PC ML are due to low empirical coverage rates for all parameters. We report the empirical coverage rates of all parameters in an online support file. across all parameters (17 angle parameters, one cosine function parameter, and 18 communality index parameters) in the model. Similarly to the results of Model I (reported in  are five methods to construct confidence intervals. "S0," "S1," "MS1," "S2," and "MS2" are five types of distributions for ordinal variables. ÃÃÃ Indicates conditions under which the construction of confidence intervals encountered difficulty. The shaded entries indicate the conditions under which the empirical coverage rates deviate from 95% (less than 92.5% or greater than 97.5%). PC OLS-A , PC OLS-B , PC ML , PS OLS , and PS ML are five methods to construct confidence intervals. "S0," "S1," "MS1," "S2," and "MS2" are five types of distributions for ordinal variables. ÃÃÃ Indicates conditions in which the construction of confidence intervals encountered difficulty.
The shaded entries indicate the conditions under which the empirical coverage rates deviate from 95% (less than 92.5% or greater than 97.5%).
(e.g., N ¼ 100). More specifically, PC OLS-B confidence intervals involve standard errors computed with Monte Carlo ACM estimates, which are less likely to encounter computation problems than sample-based ACM estimates for small samples.

Test Statistics and Type I Error Rates
We explore the statistical properties of test statistics by comparing them with their reference distributions. The reference distributions are v 2 distributions (12 degrees of freedom in Model I and 117 degrees of freedom in Model II). We denote the test statistics obtained with PC OLS-A , PC OLS-B , PC ML , PS OLS , and PS ML using "T1," "T2," "T3," "T4," and "T5," respectively. Figure 7 displays the empirical cumulative distribution functions (CDFs) of these five types of test statistics with binary variables in Model I. "T1" and "T2" are close to the reference distribution under most conditions and become closer to the reference distribution with larger samples. Therefore, PC OLS-A and PC OLS-B produce proper test statistics asymptotically even with skewed data. In contrast, "T3" substantially deviates from the reference distribution under all conditions. The patterns of "T4" and "T5" are surprising. They are close to the reference distribution at N ¼ 500 and 1,000, but they deviate from the reference distribution for skewed data at N ¼ 12,108. Such deviations from their reference distributions are more substantial when the ordinal variables are more skewed ("S2" and "MS2" versus "S1" and "MS1"). Therefore, we should be cautious in estimating the model with PS OLS and PS ML due to their biased parameter estimates and the surprising phenomenon of their test statistics. 10 Finally, PS OLS is computationally the most robust method among the five methods for small samples. For example, it successfully produced estimates for both symmetric and skewed data at N ¼ 100, but the other four methods failed under these conditions. Figure 8 compares the CDFs of "T1," "T2," "T3," "T4," and "T5" against the reference distribution for the three-category variables under different conditions in Model I. Although the overall patterns are similar to those of the two-category variables displayed in Figure 7, the CDFs of the three-category variables show two new phenomena when the sample size is 100 or 200. First, they involve fewer computational issues than two-category variables, particularly for symmetric distributions or mildly skewed distributions. The reason for this phenomenon is that threecategory variables contain more information than two-category variables. The computational issues disappear for both two-category variables and three-category variables eventually for large samples due to information increase. The second phenomenon is that "T2" CDFs are closer to the . Conditions marked with "I" indicate computing issues of test statistics (e.g., I 1,2,3 means computing issues of T1, T2, and T3). 10 We could understand this phenomenon by decomposing deviations of sample estimates from population values into two parts: random sampling error and "model error" introduced by nonsymmetric distributions. Random sampling error decreases as the sample size increases but such "model error" stays the same at different sample sizes. At small samples, random sampling error could dominate such "model error." Consequently, the statistics of "T4" and "T5" are close to their reference distributions. At larger samples, such "model error" will eventually dominate random sampling error and the test statistics will be larger than v 2 distributions. reference distribution than "T1" CDFs. The reason for this phenomenon is that the Monte Carlo ACMs of PC OLS-B are more stable than the sample-based ACMs of PC OLS-A under these conditions. The five statistics ("T1," "T2," "T3," "T4," and "T5") allow us to test circumplex models with ordinal variables. We can compute a p value of the test by comparing the sample test statistic against the v 2 distribution with 12 degrees of freedom. We reject the circumplex model if the p value is smaller than a prespecified a level (e.g., .05). Thus, the empirical Type I error rate under a certain condition is the proportion of simulated samples in which p < .05. Table 5 presents the empirical Type I error rates of the five test statistics for Model I. These Type I error rates are consistent with the CDFs shown in Figures 7 and 8. We interpret an empirical Type I error rate as "too low," "acceptable," and "too high" if it is less than 0.025, between 0.025 and 0.075, and greater than 0.075, respectively (Bradley, 1978). Overall, the best-performing test statistic is "T2" because its empirical Type I error rate is acceptable under most conditions. "T1" performs similarly to "T2" under most conditions, with several exceptions. Its empirical Type I error rates are slightly higher than those of "T2" in conditions that combine three-category variables and small samples. "T4" and "T5" perform similarly in most conditions. Their empirical Type I error rates are acceptable for symmetric variables and slightly skewed variables, but their empirical Type I error rates deviate from 0.05 with more skewed data and larger samples. In particular, their empirical Type I error rates deviate from 0.05 at N ¼ 12,108, indicating that these two statistics are not asymptotically appropriate. Finally, "T3" does not provide an acceptable empirical Type I error rate under any condition. Table 6 presents the empirical Type I error rates of the five test statistics for Model II. 11 The empirical Type I error rates of the five test statistics deviate from 0.05 more substantially in Model II than in Model I. Nevertheless, the comparisons of the five test statistics across different conditions in Model II are similar to those in Model I. The test statistic "T2" emerges as the best-performing one in both models. Its empirical Type I error rates are either acceptable or slightly higher than 0.075 for symmetric or slightly skewed data. Even for substantially skewed data, its empirical Type I error rates become closer to 0.05 as the sample size increases. Although the empirical Type I error rates of "T1" are also acceptable or slightly higher than 0.075 at N ¼ 12,108, its empirical Type I error rates are much less satisfactory than those of "T2" at other sample sizes. Such advantages of "T2" over "T1" are due to the more accurate Monte Carlo ACM estimates of PC OLS-B over the samplebased ACM estimates of PC OLS-A under these conditions. Finally, the empirical Type I error rates for "T3," "T4," and "T5" are too high under all conditions.

Discussion
We propose to estimate circumplex models with polychoric correlations to accommodate ordinal data. We consider OLS estimation due to its computational robustness with polychoric correlation matrices. We derive the corresponding standard errors and test statistics using two methods, and combining them with the OLS point estimates produces 1.00 "T1," "T2," "T3," "T4," and "T5" are five different test statistics. "S0," "S1," "MS1," "S2," and "MS2" are five types of distributions for ordinal variables. ÃÃÃ Indicates the conditions under which a test encountered computing problems in most simulated samples.
The shaded entries indicate that the empirical Type I error rates deviate from 0.05 (less than 0.025 or greater than 0.075). 11 We include the figures comparing the CDFs of the five test statistics against the reference distribution in online support file.
PC OLS-A and PC OLS-B . We compare these two methods with three additional methods (PC ML , PS OLS , and PS ML ) with both empirical data and simulated data. The two new methods produced more satisfactory confidence intervals than the other three methods under most conditions. Overall, the test statistic of PC OLS-B performed best under most conditions. The statistical properties of the two new methods improved with larger samples, but such a phenomenon did not always occur for the other three methods.
The influences of the correlation type on point estimates are more substantial for communality indices than for angle parameters. The point estimates for the communality indices obtained with PC OLS-A , PC OLS-B , and PC ML were larger than the corresponding estimates obtained with PS OLS and PS ML , but we did not observe this phenomenon for the angle parameters. These two sets of parameters play different roles in the circumplex model. Angle parameters account for correlations among the common components of variables, and communality indices indicate the contributions of these common components to their respective variables. Pearson correlations tend to be smaller in absolute value than the corresponding polychoric correlations, and the reductions are often larger for larger correlations. The circumplex model combines communality indices and angle parameters in a multiplicative way to produce variable correlations. More specifically, a variable correlation matrix is a rescaled version of the correlation matrix of their common components. If polychoric correlations are perfectly correlated with Pearson correlations, fitting the model with the two types of correlations will produce the same angle parameters but different communality indices. In practice, polychoric correlations are often highly (but not perfectly) correlated with Pearson correlations, and the magnitude reductions of Pearson correlations lead to negatively biased estimates for communality indices.
A related phenomenon is that the methods with PS OLS and PS ML produced satisfactory empirical Type I error rates for symmetric or slightly skewed ordinal variables. The null hypothesis of the test is that we can arrange the common parts of the variables in a circular form. When the magnitude reductions from the polychoric correlations to the Pearson correlation are proportional to the absolute values of the polychoric correlations, the common parts obtained with both polychoric correlations and Pearson correlations satisfy the circumplex model. Such proportional reductions are more likely to occur for symmetric ordinal variables. On the other hand, when the ordinal variables are more skewed, the estimates with Pearson correlations are more biased, and the corresponding empirical Type I error rates deviate from the nominal level.
We can relate the circumplex model to a two-factor factor analysis model, but the "factor loadings" demonstrate a circular pattern. In contrast, factor analysts prefer the factor loading matrix to be similar to Thurstone's simple structure. It is well known that factor analyzing ordinal variables as if they were continuous is inappropriate unless the number of categories is five or more (Rhemtulla et al., 2012). When a factor analysis model is appropriate for latent response variables, factor analyzing a Pearson correlation matrix of ordinal variables often leads to negatively biased estimates for factor loadings and positively biased estimates for error variances (Olsson, 1979b). We found that point estimates for communality indices are negatively biased (corresponding to positively biased estimates for error variances) when we obtain them with Pearson correlations. In addition, conducting factor analysis with Pearson correlations produces 1.00 1.00 "T1," "T2," "T3," "T4," and "T5" are five different test statistics. "S0," "S1," "MS1," "S2," and "MS2" are five types of distributions for ordinal variables. ÃÃÃ Indicates the conditions under which a test encountered computing problems in most simulated samples.
The shaded entries indicate that the empirical Type I error rates deviate from 0.05 (less than 0.025 or greater than 0.075). even more biased point estimates and inflated test statistics when the ordinal variables become more skewed and particularly skewed in different directions. In comparison, conducting factor analysis with polychoric correlations produces unbiased point estimates and proper test statistics for both symmetric and skewed ordinal variables (Forero et al., 2009;Rhemtulla et al., 2012). Similarly, PC OLS-A and PC OLS-B produced unbiased point estimates for all parameters and appropriate test statistics for both symmetric and skewed ordinal variables under most conditions. Although PC ML produced unbiased point estimates, the corresponding test statistics are greatly inflated. Dolan (1994) reported a similar result to PC ML in the context of factor analysis. Factor loading estimates were essentially unbiased, but test statistics were inflated.
Satisfactory confidence intervals and test statistics of the circumplex model with ordinal data require a sufficiently large sample. A sufficiently large sample size depends on several characteristics of the study, for example, the number of categories, the skewness of the ordinal data, and population parameter values. The simulation study shows that a sample size of 100 or 200 is large enough for symmetric ordinal variables, but a sample size of 500 or 1,000 is needed under more skewed conditions. Furthermore, PC OLS-A and PC OLS-B performed more satisfactorily than the other three methods under almost all conditions. We considered only two-category and three-category ordinal variables in the simulation study due to their substantial deviations from normal distributions. Therefore, accommodations for ordinal variables are much needed under these conditions for these data. Rhemtulla et al. (2012) conducted a simulation study to assess the influence of the number of categories on confirmatory factor analysis models and found that accommodations are needed for four or fewer categories. They also found that PS ML performed well for five or more categories. Therefore, we expect that the advantage of PC OLS-A and PC OLS-B over PS OLS and PS ML will become less substantial with more categories in the context of the circumplex model.
Standard errors and test statistics of PC OLS-A and PC OLS-B involve estimating the ACM of polychoric correlations. Although sample-based ACM estimates and Monte Carlo ACM estimates produced similar standard error estimates, the corresponding test statistics can be substantially different in small to moderately large samples. For example, the test statistics with Monte Carlo ACM estimates are much closer to the reference distribution than their counterparts with sample-based ACM estimates in Model II when the sample size is 200 and 500. The sample-based ACM estimate is the default setting in many SEM packages (e.g., lavaan [Rosseel et al., 2017], MPlus [Muth en & Muth en, 1998, and LISREL [J€ oreskog & S€ orbom, 1996]) and is recommended for large samples. A sample size of 200 or 500 seems too small for accurately estimating the ACM with the samplebased method in Model II. The Monte Carlo method seems to perform better under these conditions. The Monte Carlo method assumes that the latent response variables are normally distributed. In the simulation study, we generated ordinal variables by categorizing normal latent response variables. We introduced the asymmetry by manipulating the threshold values. The satisfactory performance of Monte Carlo ACM estimates may no longer hold when the latent response variables are nonnormal (Foldnes & Grønneberg, 2020;Monroe, 2018). Although we do not explicitly consider model error, we interpret a nonsignificant test statistic as that the amount of model error is comparable to the random sampling error. Model misspecification is unavoidable in practice. The sandwich standard error estimates of Equation (11) and the test statistics of Equation (12) are still valid with model error. The distribution of test statistics becomes a noncentral v 2 distribution. Therefore, one could conduct a test of close fit to allow some amount of model error. The model is less useful when it is severely misspecified, and researchers should modify it rather than interpreting its parameters.
Many applications of the circumplex model involve Likert variables. The proposed method allows researchers to pursue their questions in a more statistically justifiable way. In addition, the corresponding R package makes the new development readily available to researchers.
A typical element of the term o 2 Px omom 0 is Here The derivative dKl 1 l 1 dml 2 l 2