A Simulation-Based Scaled Test Statistic for Assessing Model-Data Fit in Least-Squares Unrestricted Factor-Analysis Solutions

Aims of the Proposal

The present article proposes an empirical test statistic in chi-square metric forassessing goodness of model-data fit in UFA solutions based on the ULS criterion,specifically: Principal-Axis-Factoring (PAF), ULS-MINRES, and Minimum Rank Factor Analysis (MRFA). The UFA-ULS solutions, in turn can be based on the standard linear FA model in which the variables are treated as continuous or on the non-linear UVA-FA model in which they are treated as ordered-categorical.

Our proposal is fully empirical, and avoids theoretical developments at the cost of intensive simulation which, given the capabilities of modern computers, is perfectly affordable. The basic idea is to combine: (a) the rationale of previous empirical proposals based on the idea of sampling from a simulated scenario in which the model holds exactly with, (b) non-linear transformations that bring the distribution of the resulting fit statistic close to the expected chi-square distribution. We shall label the approach as LOSEFER, as an acronym of Lorenzo-Seva and Ferrando’s approach.

Description of the Procedure

Consider a set of m observed variables related to p common factors, the population standardized variance-covariance (i.e., correlation) matrix Σ (m × m) among the set of observed variables, and the corresponding estimate R (m × m) obtained in a sample of N observations. When R is obtained from a large and representative sample from the population, R is expected to be a good estimate of Σ.

The direct UFA model decomposes Σ as,

where $\Gamma$ is the loading matrix (m × p), and $\Psi$ is the diagonal matrix (m × m) with the unique variances in the main diagonal. When an UFA solution is fitted to sample data, the aim is to estimate matrices $\Gamma$ and $\Psi$ from the observed matrix R. In terms of the sample estimate, matrix R is decomposed as,

where A (m × k) and U (m × m) are the corresponding estimates of $\Gamma$ and $\Psi$ in Expression (1), and E (m × m) represents the amount of observed covariance in R that cannot be accounted for by the sample factor model. When k (the number of factors in the sample model) is chosen to equal p (the number of factors in the population model), the observed values in E will tend to be zero, and the estimated factor model is expected to attain an appropriate goodness-of-fit. However, the typical situation when fitting a UFA solution (especially when used with exploratory purposes) is that the value of p is not known, so that k is assigned a tentative value that aims to achieve an appropriate goodness-of-fit level for the sample factor model.

Goodness-of-Fit Assessment

In order to define a scalar goodness-of fit test statistic for a prescribed UFA solution obtained from the PAF, MINRES or MRFA approaches, the matrix E is derived from Expression (2) as,

and the discrepancy function we shall consider here is:

where the $e_{ij}^2$ terms are the non-diagonal elements of E. So, the discrepancy function (4) is the sum of non-duplicated squared residuals between the observed and reproduced correlation matrix. The test statistic we consider based on this discrepancy function is now:

The discrepancy function in (4) is the ordinary or unweighted least squares (ULS) function, which is the simplest in covariance structure analysis. As mentioned above, all the methods considered here: PAF, ULS-MINRES and MRFA are essentially ULS methods, and so are based on the minimization of (4).

If the ULS estimates were asymptotically efficient, the distributional assumptions mentioned above were met, and the null hypothesis of exact fit in (1) would hold, then (5) would be asymptotically distributed as a chi-square variable with degrees of freedom,

(see e.g., Lawley, 1940).

Because the ULS estimates are never asymptotically efficient, neither for continuous nor for ordered-categorical solutions, and the fulfillment of the distributional assumptions cannot be taken for granted, a naïve theoretical chi-square reference distribution cannot be claimed for (5). To address this problem, the proposal here is to non-linearly transform the test statistic (5) so that, when the fitted solution holds in the population, the transformed statistic closely approaches a central chi-square distribution with degrees of freedom (6). Next, the sample test statistic undergoes the same transformation, so that, when null hypothesis (H₀) holds, it is interpretable as a value sampled from a chi-square distribution. Overall, then, our proposal can be regarded as a correction of a chi-square type test statistics that brings its distribution closer to that expected theoretically when the null hypothesis holds. The main basis difference with existing corrections of this type (i.e., robust statistics) is that the correction is not based on asymptotic theory but is empirical. Finally, as for requirements and assumptions for undertaking the transformation, the basic requirement is that both, the observed and reproduced covariance matrices are positive definite (see Lorenzo-Seva & Ferrando, 2021). And, as for the basic assumptions, being (4) a sum of squares that adopts only positive values, it is assumed that the distribution of (5) will be positively skewed, and will approach normality as the model becomes larger.

Empirically Obtaining the Scaled Test Statistic in Chi-Square Metric

Once an UFA solution has been fitted to an observed correlation matrix R, the reproduced variance-covariance matrix is defined as,

Let’s take this R*, as if it was a true population matrix, and decompose it using Cholesky’s method,

If Z (N × m) is a random matrix with columns normally distributed N(0, 1), the product,

produces a population matrix X (N × m). From this population matrix X, random samples X_i (N_i × m) of observed scores are sampled for which the corresponding correlation matrix R_i is an estimate of the true population matrix R*. Matrix R_i is then factor analyzed in order to obtain estimates A_i, U_i, and

Finally, the test statistic c_i is obtained using expression (5) applied to matrix E_i.

The process is repeated for an arbitrary number of times K (i = 1…K), in order to obtain a vector of c that contains the K values c_i. The distribution of the elements of this vector is then the distribution of the uncorrected test statistic when the null hypothesis holds. In the studies presented in this document, we used K = 1,000 and N = 100,000, and N_i equal to the size of the sample used to obtain observed matrix R, and arrived at acceptable results. In addition, it must be pointed out that, if X is a set of ordinal variables, each X_i must be discretized using the empirical thresholds estimated from X, and the computed correlation matrices must be based on polychoric correlations.

So far, our proposal is similar to previous existing developments, and particularly to Bollen and Stine’s (1992) bootstrapped approach (see also Corrêa Ferraz et al., 2022). The basic idea, in effect, is to sample (or resample in the bootstrapping approach) from a population in which the null hypothesis exactly holds. Furthermore, this condition is obtained by using the parent sample matrix as a basis (i.e., R* is taken as if it was a Σ in Equation 1). More specifically, the original data is transformed (according to (8) and (9) in our proposal) so that is forced to satisfy the null hypothesis.

From here on, however, our proposal differs from the previous related developments. Bollen and Stine (1992) considered only the ML-based scenario for continuous outcomes, in which the test statistic was expected to truly follow a central chi-square distribution under the null hypothesis. Second, Bollen and Stine’s (1992) proposal was not intended to make the empirical bootstrap distribution closer to the theoretical chi square distribution, but rather, to obtain reference p values to which the untransformed sample test statistic could be compared.

Continuing with our proposal, once c is available, what we propose is nonlinearly transform it using a third degree polynomial,

so that the first moment, and the second, third, and fourth central moment estimates of the transformed c coincide with those of the reference chi-square distribution with degrees of freedom (df) in (6). In more detail, the coefficients of the polynomial (11) are obtained by solving the following system,

Where E() is used for expectation. Technical details on how to determine the polynomial coefficients in (11) from the system (12) can be obtained from the authors. As a summary, of the different solving procedures we tried, the most effective was two step. First, the original c values were: (a) cube-root transformed, and (b) transformed to have the first four moments of a standard normal variable (i.e., the first four moments in the system being 0, 1, 0 and 3). Second, the normal-transformed variable was transformed again using Fleishman’s (1978) procedure to obtain a chi-square distribution from a standard normal distribution, so that the final transformed variable have the first four moments as close as possible to those in (12).

Once the coefficients a, b₁, b₂, and b₃ have been obtained, we can now factor analyze the sample correlation matrix R and compute the sample-observed c statistic by using expressions (4) and (5). Next, c is transformed to y using the fitted polynomial (11) and this transformed value is interpreted with relation to a chi-square distribution with degrees of freedom (6).

A crucial point in interpreting the transformed sample c value as a proper chi-square type test statistic is, indeed, that not only the first four moments, but its distribution in general would adhere to the theoretical distribution under the null hypothesis of model-data fit. Strictly speaking, this adherence cannot be guaranteed, so we assessed this point using intensive simulation. Full results are provided below. However, it can be advanced that our proposal shows a considerable viability in this respect.

Beyond the Scaled Test Statistic: GOF Indices Test of Close Fit and Power Analysis

Provided that y approaches closely enough the corresponding reference distribution, it can be further used as a basis for computing meaningful point estimates of selected GOF statistics, so that the proposed solution can be more thoroughly assessed. We shall propose to derive CFI and RMSEA point estimates directly from the transformed y statistic. However, it should be clear that this is only an initial proposal that is expected to be updated as more information about the performance of GOF indices in UFA will become available.

Let λ₁ = y₁ - df₁ the noncentrality parameter estimate for the solution under study, and λ_o = y_o - df_o the corresponding estimate for the solution with zero common factors (i.e., the null or baseline model). In terms of y, the comparative fit index point estimate can then be obtained as:

And the RMSEA point estimate as:

Our view, however, is that meaningful information of a GOF statistic requires not only the point estimate to be reported, but also the corresponding confidence interval. In the implementation approach in which the indices proposed here are programed (see below), 90% confidence intervals are reported based on bootstrap resampling.

The choice of the RMSEA as an index directly derived from the test statistic allows two further pieces of important information to be obtained: the test of close fit and power analysis (Lee et al., 2012). In our view, this information is highly relevant for avoiding two common pitfalls when fitting UFA solutions. The first is to use too small samples in order to achieve a better fit (at the cost of a gross loss of power). The second is to over-factor with the same aim, a practice that ends up in considering as relevant trivial, uninterpretable minor factors which are devoid of any substantive interest (Ferrando & Lorenzo-Seva, 2018).

The implementation is straightforward. In our proposal we have set the null (close fit) and alternative hypothesis as,

With regards to power assessment, we have chosen the approach by Lee et al. (2012) in which the noncentrality parameter that expresses the lack of fit in the population is obtained by setting RMSEA values under H₀ and H₁. In particular, in our implementation, power is computed as the capacity for distinguishing between a close fit solution (H₀: RMSEA = 0.05) and a moderately misspecified solution (H₁: RMSEA = 0.08).

First Simulation Study

A Monte Carlo simulation study was carried out using samples drawn from a true population model. Based on Expression (1), a population loading matrix was defined in which each observed variable had a salient loadings had a values of .70, and unicity equal to 1 minus communality of the variable. The number of factors and the number of salient variables per factor were manipulated in order to produce models with different degrees of freedom:

The sample sizes were also manipulated: 200, 500, and 800. From each true population matrix, 500 samples were obtained, and a total of 6,000 were factor analyzed.

For each sample, the true number of factors was extracted using two extraction methods: ULS/Minres, and Principal Axes (PA); and the y statistic computed after each extraction. The statistic y was computed using population samples of 100,000 and the number of random samples extracted from the population were K = 1,000.

The four firsts moments of the empirical distribution of y were computed and compared to the expected ones for the theoretical distribution of chi square statistic with degrees of freedom df. Table 1 shows the outcomes related to ULS/MINRES extraction. Kurtosis values are printed as zero centered (i.e., kurtosis minus 3). Mean and variances were slightly overestimated, especially in small samples. The estimates of skewness and kurtosis in general do not differ from the values expected in the population, and only in small models (5 degrees of freedom) and large samples (N = 800) the estimates are overestimating the expected values.

Rejection rates after ULS/MINRES factor extraction are shown in Table 2. The worse rejection rates were observed when the sample size was small (N = 200). In addition, the worse rejections rates estimates were the ones expected to be .001. It means that is at the farthest tail of the distribution of statistic y is where less adherence to the chi-square distribution is observed. From a practical point of view, the rejection levels observed are reasonable.

Table 3 shows the outcomes related to PA extraction. Mean and variances were slightly underestimated, especially in large samples. Again, the estimates of skewness and kurtosis in general do not differ from the values expected in the population, and only in small models (5 degrees of freedom) and large samples (N = 800) the estimates are underestimating the expected values.

Rejection rates after PA factor extraction are shown in Table 4. The outcomes are quite similar to those obtained after ULS/MINRES extraction pattern. However, PA extraction obtained slightly better rejection rates than ULS/MINRES. Again, rejection rates improve when sample sizes are large.

Finally, RMSEA, CFI and NNFI goodness-of-fit indices were computed. As the model that was fitted to the sample data systematically corresponded to the true population model, the values of these indices should indicate in all cases that an acceptable model fit had been attained. Table 5 shows the mean of goodness-of-fit indices after ULS/MINRES factor extraction. As can be observed, a good model fit was always reported. In addition, as the sample size became larger, the goodness-of-fit values improved for all the indices.

Table 6 shows the mean of goodness-of-fit indices after PA factor extraction. Once again, a good model fit was always reported. It must be pointed out that the values obtained here suggested a better fit than those obtained after ULS/MINRES extraction. In addition, the estimates did not seem so much influenced by the sample size, as it happened after ULS/MINRES extraction.

The conclusion of the first simulation study is that, when the correct model is fitted to the sample data, the distribution of the y statistic shows a reasonable adherence to the expected chi-square distribution, and the related goodness-of-fit indices can be safely interpreted.

Second Simulation Study

The second study is mainly concerned with assessing the power and sensitivity of our proposed statistic. To assess so, we replicated the first simulation study with two variations. First, the number of factors retained in the sample data was one less than the true number in the population. Second, only three population models were considered: the models with 64, 207 and 492 degrees-of-freedom in the population.

ULS/MINRES and PA extraction methods were used to extract the incorrect number of factors. As the sample model was not the one that existed in the population, goodness-of-fit indices derived from y should suggest an improper model adjustment. In addition, the values should worsen when the discrepancy between the sample and the population models increases. Thus, if the population model had four factors, and the sample model is adjusted for three factors, the misspecification is not so strong as if the population model had two factors, and the sample model is adjusted for one single factor. The values of the goodness-of-fit indices should be sensitive to the different degrees of misspecification.

RMSEA, CFI and NNFI goodness-of-fit indices were again computed. Table 7 shows the mean of the goodness-of-fit indices after ULS/MINRES factor extraction. As can be observed, an unacceptable model fit was always reported for all the indices, and the worse goodness-of-fit values were related to the solutions with the lowest degrees of freedom. It must be noted that, when the sample size was large, the indices reported better goodness-of-fit values, but these values were still farther away from the usual threshold values used in applied research for judging the fit as acceptable.

Table 8 shows the mean of goodness-of-fit indices after PA factor extraction. As can be observed, the values of the goodness-of-fit are quite similar to the ones obtained after ULS/MINRES extraction.

The conclusion of the second simulation study is that, when the sample model is incorrectly specified, the goodness of fit indices derived from the proposed y statistic are expected to detect that the proposed model is wrong under most of the conditions expected to occur in practice.