Should Regularization Replace Simple Structure Rotation in Exploratory Factor Analysis?

Exploratory factor analysis (EFA) is an important tool when the measurement structure of psychological constructs is uncertain. Typically, factor rotation is applied to obtain interpretable results resembling a simple structure. However, an overwhelming multitude of rotation techniques is available of which none is unequivocally superior. Recently, regularization has been suggested as an alternative to factor rotation. In two simulation studies, we addressed the question if regularized EFA is a suitable alternative for rotated EFA. We compared their performance in recovering predefined factor loading patterns with varying amounts of cross-loadings. Elastic net regularized EFA yielded estimates comparable to rotated EFA. For complex loading patterns, both rotated and regularized EFA tended to underestimate cross-loadings and inflate factor correlations, but regularized EFA was able to recover loading patterns as long as a subset of items followed a simple structure. We conclude that regularization is a suitable alternative to factor rotation for psychometric applications.

Exploratory factor analysis (EFA) is one of the most commonly used statistical methods in psychological research. EFA allows researchers to summarize the observed data (e.g., item responses) as a function of a few latent variables (e.g., traits), typically called factors (see, e.g., Mulaik, 2010, for a general introduction). After an initial solution has been estimated, factor rotation is used in order to obtain a more interpretable solution. A variety of rotation techniques have been proposed for this purpose including, for example, Varimax, Geomin, or Quartimax (see Browne, 2001;Mulaik, 2010, for overviews). However, the multitude of rotation techniques makes it difficult for applied researchers to choose an appropriate technique for their application scenario. Furthermore, research found that rotation techniques differ in their ability to uncover a known population factor structure and that the performance of a specific rotation technique depends on the population pattern itself Sass & Schmitt, 2010;Schmitt & Sass, 2011). As the population factor structure is unknown in practical EFA applications, the choice of the rotation technique is a purely subjective step.
Regularized (or sparse) EFA has been suggested as an alternative to the factor rotation step (e.g., Trendafilov, 2014;Yamamoto, Hirose, & Nagata, 2017). Instead of rotating factor loadings, regularized EFA tries to achieve a more interpretable solution such as a simple structure by penalizing factor loadings and/or factor correlations directly in the estimation stepshrinking non-substantial parameters toward zero. Regularized EFA addresses the subjectivity of the rotation approach to some extent, because the tuning parameters that are used for penalization can be determined in a more objective way, for example, by drawing on information criteria such as the Bayesian information criterion (BIC; e.g., Hastie, Tibshirani, & Friedman, 2009;Jacobucci, Grimm, & McArdle, 2016;James, Witten, Hastie, & Tibshirani, 2013).
Despite the growing body of methodological literature on regularized EFA, including illustrative examples of the potential of regularized EFA as a substitute for factor rotation (Trendafilov, 2014;Yamamoto et al., 2017), applications of regularized EFA are rare. One reason for this might be that applied researchers are lacking sufficient information to judge the usefulness of regularized EFA for their specific use case. Previous research focused on direct comparisons of either rotation techniques (e.g., Schmitt & Sass, 2011) or regularization methods (e.g., Hirose & Konishi, 2012), but extensive direct comparisons of factor rotation and regularization across many data situations are not yet available. The present article aims to fill this gap by comparing the performance of different rotation techniques with different regularization techniques with regard to parameter estimation. In addition, we investigated whether only the factor loadings or both factor loadings and factor correlations should be treated as to-be-regularized parameters in regularized EFA.
In the following, we first describe how the parameters of the EFA model are conventionally estimated by factor rotation techniques. Then, we explain how regularized EFA parameters can be obtained using penalized maximum likelihood (ML) estimation with a least absolute shrinkage and selection operator (lasso), ridge, or elastic net (enet) penalty (e.g., Hastie et al., 2009). Thereafter, we present the results of two simulation studies comparing the performance of factor rotation and regularization in the estimation of factor loadings and factor correlations. In addition, we investigated the influence of the sparsity of the factor loading pattern on the recovery of the population parameters.

FACTOR ROTATION IN EFA
The EFA model describes the p observed variables as a weighted linear combination of m factors (e.g., Mulaik, 2010): where Y is the p Â n matrix of observed variables measured from n observations, η is the m Â n matrix of factor scores, Λ is a p Â m matrix of factor loadings, and is a p Â n matrix of error terms. The parameters of the EFA model are often estimated with a ML approach. Here, estimates are obtained that minimize the discrepancy between the modelimplied covariance matrix of the variables Σ and their observed covariance matrix S (Jöreskog, 1969): F ML ðΣ; SÞ ¼ ln Σ j j þ trðΣ À1 SÞ À ln S j j À n It is well known that the EFA model is rotationally indeterminate. That is, an infinite set of equally well fitting factor solutions exists for a given data set that may be transformed into each other by a rotation matrix H (e.g., Mulaik, 2010, p. 276): Here, ΛH are the rotated factor loadings and H À1 η are the rotated factor scores. Consequently, researchers have to 'choose' a factor solution that describes the observed data in a convenient (i.e., interpretable) way. This is achieved by rotating the initial model, that is, computing an optimal rotation matrix H according to a prespecified criterion. Most of the common rotation techniques optimize a simple structure criterion of the factor loadings requiring that each variable should load highly onto one factor and should have low cross-loadings onto the other factors (Browne, 2001;Thurstone, 1935Thurstone, , 1954. About 50 rotation techniques have been proposed in the methodological literature (Trendafilov, 2014). The most apparent distinction can be made between orthogonal rotation techniques (e.g., Varimax; Kaiser, 1958Kaiser, , 1959, which constrain the factor correlation to zero, and oblique rotation techniques (e.g., Geomin; Yates, 1987), which allow the factors to be correlated. Apart from that, rotation techniques differ in the exact criterion that is used to operationalize a simple factor loading pattern (e.g., Browne, 2001), and hence in their tolerance for cross-loadings (Schmitt & Sass, 2011). In the following, we will briefly contrast three common rotation techniques: Varimax rotation, Geomin rotation, and Facparsim as a member of the general Crawson-Ferguson (CF) rotation family.
Varimax (Kaiser, 1958(Kaiser, , 1959 is one of the most widely applied rotation techniques (Fabrigar, Wegener, MacCallum, & Strahan, 1999). It assumes uncorrelated factors and maximizes the variance of the squared loadings: An initial Varimax rotation is an essential part of the popular oblique Promax criterion (Hendrickson & White, 1964). The Geomin rotation was explicitly developed in order to represent more complex loading patterns. It minimizes the (variable-wise) geometric mean of the squared factor loadings (Browne, 2001): Here, is an additional parameter to ensure that the rotation criterion is generally differentiable even if one of the factor loadings is exactly zero for each variable. Traditionally, ¼ 0:01 is used but a slightly higher value of ¼ 0:5 has been suggested in the literature to better represent more complex factor structures (Marsh et al., 2010(Marsh et al., , 2009Morin, Marsh, & Nagengast, 2013). The CF rotation criterion considers both variable-wise and factor-wise complexity in the rotation: The first and second term in this sum quantify the row (i.e., variable) and column (i.e., factor) complexity, respectively, in the factor loading matrix. The weight k ¼ ½0; 1 determines which complexity receives more emphasis during factor rotation with higher values indicating more emphasis on factor complexity. Many rotation criteria may be described as special cases of the CF criterion with different values of k (e.g., Sass & Schmitt, 2010). Due to the different simple structure criteria, rotation techniques differ in their ability to recover population patterns that differ in their amount of cross-loadings. More specifically, simulation studies showed that orthogonal rotations of oblique factor patterns yield spurious cross-loadings even if the population factor loading pattern is a simple structure (Schmitt & Sass, 2011). Furthermore, oblique rotations tend to yield inflated factor correlations and underestimated cross-loadings in the presence of high cross-loadings in the population factor loading pattern Sass & Schmitt, 2010;Schmitt & Sass, 2011). Thus, factor rotation techniques can achieve a more unique assignment of variables to factors at the cost of less distinct factors (or vice versa), and it is desirable to achieve a reasonable trade-off between factor correlation and cross-loadings.
To summarize, factor rotation is utilized to find a factor solution that is as interpretable as possible according to some simplicity criterion. However, applying rotation techniques has two major drawbacks: First, researchers have to choose among a multitude of rotation techniques that is still growing (e.g., Beauducel, 2018;Ertel, 2011;Jennrich, 2004Jennrich, , 2006Yamamoto & Jennrich, 2013). In addition, some rotation techniques have tuning parameters (e.g., k and mentioned above) that also have to be chosen (and that have profound consequences for parameter estimation). Second, the performance of a rotation technique in terms of the suitability of the parameter estimates depends on the true data generating mechanism in the population, especially on the amount of cross-loadings, which is unknown in practice. Therefore, alternative approaches that perform consistently well across a wide range of factor patterns would be preferableobviating the need for such a subjective choice.

REGULARIZED EFA
The rotation problem in EFA can be reconceptualized as a variable (or model) selection problem in which a set of indicator variables for each factor needs to be chosen among all indicators (e.g., Hirose & Konishi, 2012;Hirose & Yamamoto, 2014). A lot of work in the SEM literature has been done on efficient model selection via heuristic search algorithms in the context of model modification (e.g., Glymour, Madigan, Pregibon, & Smyth, 1997;Marcoulides & Drezner, 2003;Marcoulides, Drezner, & Schumacker, 1998;Marcoulides & Falk, 2018;Marcoulides & Ing, 2013). Essentially, these approaches are best subset selection methods that try to go through the space of possible models as efficiently as possible but still apply conventional estimators (such as ML) to estimate the model parameters. The best model is then selected based, for instance, on the BIC (Schwarz, 1978). In contrast, regularization strives for a solution in which as many parameters as possible are (close to) zero directly during the estimation of the model with the goal that only few but substantial variables remain in the final model (Hastie et al., 2009).
A vector of parameters in which most of the entries are zero is called sparse. With respect to EFA, where the parameters to be estimated are the factor loadings and factor correlations, a perfect simple structure of the factor loadings can be seen as a special case of sparsity, and this is the reason why regularized EFA has been suggested as an alternative to factor rotation (Trendafilov, 2014). However, unlike the simple structure criterion, the sparsity condition does not refer to a specific pattern of the non-zero estimates within the set of parameters (e.g., that cross-loadings are zero). In that sense, a conventional rotated EFA and a regularized EFA both aim for a simple factor loading pattern but the latter does not distinguish between factor and variable complexity (see Hirose & Yamamoto, 2015a;Yamamoto et al., 2017, for more formal treatments of this notion). This property may enable regularized EFA to flexibly recover a larger variety of factor loading patterns than rotated EFAchallenging the predominance of rotation in typical psychometric applications.
Several variations of regularized EFA have been proposed in the literature (e.g., Arruda & Bentler, 2017;Hirose & Konishi, 2012;Hirose & Yamamoto, 2014;Huang, Chen, & Weng, 2017a, 2017bJacobucci et al., 2016;Jung & Lee, 2011;Jung & Takane, 2007;Trendafilov & Adachi, 2015;Trendafilov, Fontanella, & Adachi, 2017). Here, we focus on approaches that directly add a penalty term at the ML-estimation stage with the aim of finding a sparse measurement model (Hirose & Yamamoto, 2014;Jacobucci et al., 2016). That is, the parameters of the model are estimated by minimizing a penalized version of the ML fit function (Eq. 2): 578 SCHARF AND NESTLER F regEFA ðΣ; SÞ ¼ F ML ðΣ; SÞ þ α Á PðθÞ Here, PðθÞ is a penalty function of a vector of model parameters θ that may, in principle, contain any parameter of the model, that is, factor loadings or factor correlations. The tuning parameter α determines the amount of penalty applied during estimation and needs to be determined in a separate step. In general, the penalty term will increase as a function of the number of nonzero parameter estimates so that the estimation procedure prefers models with many low or zero parameter estimates. In that respect, rotation and regularization have similar objectives (see Trendafilov, 2014;Yamamoto et al., 2017, for illustrative examples). The vector of penalized parameters θ may contain the factor loadings (i.e., θ ¼ vec ðΛÞ) or both the factor loadings and the factor correlations (i.e., θ ¼ ½vec ðΛÞ φ, where φ denotes a vector containing all factor correlations). 1 Furthermore, a variety of different penalty functions have been proposed in the literature including ridge, lasso, and enet penalties. Both ridge (Hoerl & Kennard, 1970) and lasso (Tibshirani, 1996) penalties are based on vector norms of the parameter vectors. Specifically, ridge uses the sum of the squared parameter estimates as penalty term, while lasso penalizes the sum of the absolute values of the parameter estimates: Here, k Á k denotes the respective norm operator and the sum is taken across all parameters contained in θ. Both penalties result in a shrinkage of the parameters toward zero but only lasso can shrink the parameter estimates to exactly zero (i.e., for α ! 1), allowing for variable selection (Hastie et al., 2009). Importantly, the variable selection property of regularization also removes the rotational indeterminacy so that regularized EFA solutions are unique (except for reordering of factors and sign switches; Choi, Oehlert, & Zou, 2010). Another consequence of the penalty term is that regularized EFA solutions tend to fit the data slightly worse than rotated EFA solutions (e.g., Jin, Moustaki, & Yang-Wallentin, 2018;Trendafilov et al., 2017).
The ability to conduct variable selection is an advantage of lasso over ridge penalization. However, lasso regression performs worse than ridge if the number of parameters exceeds the number of observations by far (Zou & Hastie, 2005). For these situations, the enet penalty has been proposed which applies both a lasso and a ridge penalty. That is, enet considers both the sum of the absolute values of the parameter estimates and the sum of the squared parameter estimates: Enet can be seen as a generalization of lasso and ridge. It includes an additional weight parameter β that determines which of the penalties receives more weight. Notably, when β ¼ 1 or β ¼ 0, enet is equivalent to ridge and lasso, respectively. When 0 < β < 1, continuously less parameters are shrunken to exactly zero the more β approaches 1 (Hastie et al., 2009). For the sake of completeness, it should be mentioned that further penalty functions have been proposed with the explicit goal of achieving sparser solutions than factor rotation (Fan & Li, 2001;Hirose, 2016;Hirose & Yamamoto, 2015b;Zhang, 2010). This is especially important for data sets with a very large number of variables (relative to the sample size) such as genome data (e.g., Carvalho et al., 2008), or fMRI data (e.g., Hirose, 2016). However, as outlined in the context of factor rotation, simpler (or sparser) solutions are typically accompanied by inflated factor correlations for psychometric data sets. Considering even sparser solutions would rather aggravate the outlined problems. Therefore, in the present paper, we focus on penalties that do not specifically aim for sparser solutions than factor rotation (i.e., ridge, lasso, and enet).
Apart from the choice of the penalty function, the tuning parameter α heavily influences the estimates of regularized EFA. In general, the parameter estimates and the proportion of nonzero parameters decrease as α increases. In that sense, regularization may be seen as a continuous and objective approach to achieve a simple (or sparse) solution as the tuning parameter is determined utilizing an objective criterion. Typically, the model is estimated over a range of possible values for α, and the set of parameters is chosen that yields the best cross-validated fit as indicated by the root mean squared error or some information criterion (Hastie et al., 2009) such as the BIC. For regularized factor analysis models, the BIC performs well in finding a penalty weight that results in reasonable parameter estimates (Hirose & Yamamoto, 2014;Jacobucci et al., 2016).
In summary, regularized EFA aims for a sparse factor loading matrix with as many zero-elements as possible but without assuming a specific form of simple structure (e.g., that each variable should have at least one zero loading). This makes regularization a potential alternative to factor rotation also for psychometric applications (see also Trendafilov, 2014). Importantly, it should not be expected that regularized EFA always performs better than rotated EFA -although there may be conditions under which regularized EFA is generally superior in terms of parameter estimation. Rather, regularized EFA may provide a better compromise between the ability to uncover simple structure when it exists with the ability to offer reasonably interpretable results when this is not the case. In that sense, regularized EFA could already be considered a suitable alternative for factor rotation if its estimates are not substantially worse than the best rotated EFA for a given factor loading pattern. Previous research convincingly demonstrated that regularized EFA achieves interpretable solutions for some specific applications (e.g., genome data or popular standard examples for EFA). Extending these studies, we investigated if this observation can be generalized over a wider range of populations inspired by typical psychometric data sets in which we systematically varied the size of the cross-loadings.

THE PRESENT STUDY
Extensive empirical comparisons of factor rotation techniques for typical psychometric data sets are available Sass & Schmitt, 2010;Schmitt & Sass, 2011) and a number of studies have compared different penalty functions in the context of both exploratory factor analysis and regression (Fan & Li, 2001;Hirose, 2016, Hirose & Yamamoto, 2014, 2015bHuang et al., 2017a;Zhang, 2010). Occasionally, rotation techniques have been included in simulation studies on regularization but these comparisons were limited to one specific rotation technique per study and considered only a small variety of factor loading patterns (Hirose & Yamamoto, 2015b;Ning & Georgiou, 2011;Trendafilov & Adachi, 2015, 2015.
Hence, more direct and extensive comparisons of factor rotation and regularization on the same data set are necessary to judge the usefulness of regularization for typical psychometric applications. Closing this gap, we report the results of two simulation studies comparing the performance of the described regularization methods with factor rotation. In the first simulation, we compared the asymptotic properties of factor rotation and regularization for large samples, and, in the second study, we investigated whether the results from the large samples also apply to samples sizes that are more realistic in psychological research.

STUDY 1: ASYMPTOTIC PERFORMANCE
In this study, we compared the performance of factor rotation and regularization for large samples. We based our comparison on factor loading patterns for which the performance of common factor rotation techniques is known. In addition, we investigated the performance of regularization and factor rotation on extended factor loading patterns with more items and varying degrees of sparsity. We focused on oblique CF-Varimax rotated EFA due to its popularity, and on Geomin and Facparsim rotated EFAs as comparison techniques because they performed best for complex loading patterns in previous simulations Schmitt & Sass, 2011). We estimated regularized EFAs with ridge, lasso, and enet penalties. We either penalized the factor loadings only or both the factor loadings and the inter-factor correlations.
As outlined earlier, the success of factor rotation depends on the degree to which the population factor loadings follow a simple structure. Regarding the performance of regularization, it is an important property, especially of lasso penalties, that the performance of the regularization depends on the degree to which the true model is actually sparse, (Donoho, 2006;Donoho & Stodden, 2006). Only if a sufficient degree of sparsity holds in the population, regularization is likely to recover the correct factor loading pattern. In the present context, the factor loading patterns are sparser the closer they approximate a perfect simple structure. Therefore, we expected that regularization behaves similarly as factor rotation in case of a population pattern that conforms to a simple structure. For more complex factor loading patterns, the cross-loading estimates are shrunken toward zero at the cost of inflated factor correlations. We expected that this effect is reduced if the penalty also considers the inter-factor correlations. In this case, inflated factor correlations would result in a higher penalty term, which in turn should allow the estimation procedure to aim for a better compromise between minimal cross-loadings and minimal factor correlations.

Simulation model
We simulated factor loading patterns of varying complexity with 3 factors and 18 (basic conditions) or 36 variables (extended conditions). For the sake of comparability, we adapted the simulation setup that was used in Schmitt and Sass (2011, Table 1) comprising 18 variables that follow a perfect simple pattern, an approximate simple pattern with small cross-loadings (< 0.20) or a complex pattern with substantial cross-loadings of up to 0.40 (standardized loadings). The standardized main loadings varied between 0.63 and 0.75 and the factors were standardized (variance of 1) and substantially correlated (0.40). The residual covariance matrix was a diagonal matrix, that is, the residuals were uncorrelated.
In addition to these basic conditions, we also explored the behavior of factor rotation and regularization for 5 extended factor loading patterns with 36 variables. These patterns were constructed by concatenating different combinations of the basic patterns (Table 2). These extended conditions enabled us to investigate the influence of the number of items and the degree of sparsity on the performance of factor rotation an regularization. For instance, in the condition 'Extended Simple 1' composed of a Perfect Simple pattern and an Approximate Simple pattern, that is, the variables 1 to 18 loaded on the 3 factors with 6 main loadings (0.75) per factor and otherwise zero-loadings (Table 1, left-most pattern), and the variables 19 to 36 loaded on the 3 factors with small cross-loadings (Table  1, middle pattern). In this pattern, both factor and variable complexity are higher than in the basic Perfect Simple pattern but the sparsity is preserved to some extent due to the zero cross-loadings of the first 18 variables. Hence, this condition is optimally suited to investigate differential behavior of rotated and regularized EFA.
Finally, in order to exclude that differences between basic and extended conditions may be attributed to the increased number of variables, we included two conditions (Extended Simple 2, Extended Complex 3) that differed to the respective basic condition only in the number of variables.

Procedure
All simulations and analyses were conducted in R (Version 3.4.4, R Core Team, 2018). All scripts necessary to reproduce the simulations and analyses are available from the Open Science Framework. For each factor loading pattern, we derived the implied covariance matrix of the variables from the common factor model (e.g., Mulaik, 2010, p. 138). In order to investigate the asymptotic behavior of the discussed approaches, we drew large random samples of N ¼ 10000 participants from a continuous multivariate normal distribution using the package mvtnorm (Genz et al., 2017).
For each condition, the respective data set was subjected to rotated ML-EFAs and regularized EFAs, extracting three factors. ML-EFA was conducted as implemented in the package psych (Version 1.7.5; Revelle, 2016) and rotated using oblique CF-Varimax, oblique Geomin ( ¼ 0:01 or Note. The factor correlations were .40 among all factors. These simulation parameters were adapted from Schmitt and Sass (2011; where the same patterns are used but unstandardized loadings are presented). The factor loadings were standardized with respect to the total variable variances (Muthén, 2004, Appendix 3). F Á = Factor.  (Bernaards & Jennrich, 2005). Regularized EFA was conducted using the package regSEM, a general package that estimates regularized structural equation models (regSEMs) allowing the user to select which parameters of the model should be penalized (Jacobucci, 2017). In order to estimate a regularized EFA, we specified three measurement models but no structural model. Each variable was allowed to load on each factor, the factor variances were fixed to 1, and all factor correlations were freely estimated. Both rotated and regularized EFA were conducted on z-standardized data.
We compared the performance of ridge, lasso, and enet penalty on either the factor loadings alone (Ridge Λ , Lasso Λ , Enet Λ ) or the factor loadings and factor correlations (Ridge Λ;Φ , Lasso Λ;Φ , Enet Λ;Φ ). The tuning parameters α and β (for enet, Eq. 11) were automatically chosen so that they minimized the BIC over the respective sample (Jacobucci et al., 2016). For α, a grid of 100 values starting from α ¼ 0:001 with a step size of 10 À5 was used. For all models, we ensured that the final parameter estimate was not at the boundary of the grid (indicating that the parameter space should be enhanced). For β, we tested values between 0.05 and 0.95 with a step size of 0.05.

Dependent measures
The standardized root mean residual (SRMR) was calculated as measure of fit between the model-implied and the observed correlation matrices following the procedure described, for instance, by Asparouhov and Muthén (2018, Section 2.2). Before calculating the dependent measures describing the recovery of the population parameters, factors were inverted if the sum of their loadings was negative (e.g., Asparouhov & Muthén, 2009, Appendix D) and factor alignment was ensured by reordering the factors according to their highest Tucker congruency with the respective factors in the population loading pattern (e.g., Lorenzo-Seva & Ten Berge, 2006). The average congruency across all three factors also served as measure of similarity between the estimated and the population loadings. In addition, the average bias (across all factors and variables) was calculated separately for main loadings, cross-loadings, and the inter-factor correlations.
To avoid misunderstandings, it should be reiterated that due to the rotational indeterminacy of the EFA model an infinite set of factor loadings and factor correlations has the same fit for a given correlation matrix. Hence, the concept of parameter bias does not apply in the same manner to the EFA model as to other models (e.g., linear regression). Following conventions of previous studies Sass & Schmitt, 2010;Schmitt & Sass, 2011), we define bias as the deviation of the estimated parameters from the data generating parameters. The main purpose of the reported bias measures is to summarize the estimated factor loading patterns in an efficient way, and they should not be understood as deviations from a ground truth.

Results
The main results of this simulation are summarized in Table 3. The model fit as indicated by the SRMR was nearly perfect across all methods but slightly worse for lasso and enet regularized EFA than for rotated or ridge regularized EFA. The estimated factor loadings and factor correlations for each condition are available from the OSF. Except for ridge penalized regularized EFAs, the congruencies and biases of the main loadings indicated that all investigated methods recovered the general factor pattern sufficiently. However, we observed strong differences between the methods with respect to the biases in the cross-loadings and factor correlations. Notably, across all conditions, methods, and factors, the bias in the cross-loadings and the bias in the factor correlations were strongly correlated, r spearman ¼ À 0:98, indicating that the bias in the factor correlation reliably increased the more the cross-loadings were underestimated.

Basic conditions
Geomin ð ¼ 0:01Þ rotation performed very well for conditions with simple structure in the population, yielding unbiased estimates of main loadings, cross-loadings and factor correlations. However, in the presence of substantial cross-loadings, Geomin ð ¼ 0:01Þ underestimated the cross-loadings and overestimated both main-loadings and factor correlations. This pattern was more pronounced the more complex the factor pattern was (Approximate Simple vs. Complex condition). The alternative tuning parameter in Geomin ð ¼ 0:5Þ had profound influences on the performance of the Geomin rotation, resulting in much weaker biases for Approximate Simple and Complex patterns but performing worse for Simple patterns where it introduced spurious cross-loadings and underestimated factor correlations. A similar trend as for Geomin ð ¼ 0:5Þ was observed for Facparsim and CF-Varimax rotations, but both were slightly less biased for Approximate Simple and Complex factor patterns and more biased for Perfect Simple patterns.
Lasso and enet regularization estimates (see Table A1 for an overview of the enet β weights) were very similar to Geomin estimates. They perfectly recovered the Simple factor loading pattern, and resulted in underestimated crossloadings and inflated factor correlations for Approximate Simple and Complex patterns. Notably, lasso and enet were less biased compared to Geomin ð ¼ 0:01Þ but more biased than Geomin ð ¼ 0:5Þ or Facparsim in the presence of cross-loadings. Ridge penalized EFAs behaved notably different yielding more biased estimates than all other tested methods for Perfect and Approximate Simple 582 SCHARF AND NESTLER patterns. Complex factor patterns, however, were almost perfectly recovered by ridge regularization, especially when the factor correlations were included in the penalty term (Ridge Λ;Φ ). Overall, the factor correlations tended to be less inflated for Complex factor patterns when they were included in the penalty term. However, the differences between regularization of only the factor loadings or factor loadings and factor correlations were rather smallexcept for the Ridge Λ;Φ .

Extended conditions
In general, we observed similar trends in the extended conditions, that is, the higher the cross-loadings of a respective pattern, the more difficult it was for the majority of the investigated methods to recover the population pattern. Especially, Geomin ð ¼ 0:01Þ rotation yielded severely biased estimates in all extended conditions. As in the basic conditions, crossloadings were underestimated and factor correlations overestimated, and this pattern was more pronounced the higher the cross-loadings. Notably, we observed that even a set of additional items with perfect simple structure (Extended Simple 1, Extended Complex 1) did not reduce the biases sufficiently for Geomin ð ¼ 0:01Þ. Geomin ð ¼ 0:5Þ and Facparsim were clearly less biased in the presence of cross-loadings butas in the basic conditionssuffered from spurious cross-loadings in the simpler conditions (Extended Simple 1).
Unlike all investigated rotation techniques, lasso and enet penalties recovered the factor loading patterns in the extended conditions perfectly if the additional variables followed a perfect simple structure (Extended Simple 1, Extended Complex 1). In all other extended conditions, lasso and enet yielded similar estimates as in the basic Complex condition, that is, cross-loadings and factor correlations were biased to a similar degree. In contrast to enet and lasso, ridge penalties resulted in fairly distorted estimates in the extended conditions. Remarkably, the ridge penalty additionally distorted the main loadings by a substantial amount. 2

Discussion
In the first simulation study, we compared the asymptotic performances of regularized EFAs and traditional rotated EFAs. For factor rotation, we replicated previous simulation results indicating that the performance of factor rotation depends on the combination of factor rotation technique, rotation parameter (here: ) and population factor loading pattern. In line with previous notions (Morin et al., 2013), a modified Geomin ð ¼ 0:5Þ criterion and Facparsim rotation performed especially well in conditions with moderate to high variable complexity but performed poorly in conditions with simple structure. This was also the case for oblique CF-Varimax rotation. The performance of regularized EFA depended largely on the choice of the penalty function: Lasso and enet recovered the population pattern if it contained a sufficient amount of zero-loadings and otherwise yielded similar estimates as factor rotation. Overall ridge penalties were less successful than lasso and enet. While ridge penalties were superior in some selected conditions, they resulted in severe distortions in some other conditions.
In sum, the advantages of ridge for complex loading patterns wereby faroutweighed by the distortions in other conditions and the inability to recover simple structure in the population. Enet and lasso, however, recovered simple structure where it existed and were able to handle the extended conditions in which additional simple structure items were appended. These results are in line with previous studies showing that lasso (but not ridge) asymptotically selects the correct subset of variables (here, items as indicators of the factors) if the sparsity assumption holds, that is, if a sufficient proportion of the parameters to be estimated is zero in the population (Donoho, 2006;Donoho & Stodden, 2006). Put simply, lasso and enet strive for the sparsest parameter matrix and do not distinguish between factor and variable complexitywhich enables them to recover the extended factor loading patterns in which there was high variability with respect to variable complexity. Across all conditions, an enet penalty on factor loadings and factor correlations showed slightly better performance than a simple lasso penalty, indicating that it combined the strengths of lasso and ridge to some extent.
With respect to the question if regularization is a suitable alternative to factor rotation, the present results are promising because, even in the conditions where the sparsity assumption was violated, enet was able to match up with factor rotation techniquesoutperforming traditional Geomin ð ¼ 0:01Þ rotation in every single of the tested conditions. Compared to the modified Geomin ð ¼ 0:5Þ and Facparsim rotations, enet performed worse for patterns with moderate cross-loadings (Approximate Simple) and comparably for patterns with high cross-loadings (Complex). However, it should be considered that both rotations achieved their partial superiority (similarly to ridge regularization) at the cost of an inability to recover the loadings of simple structure factor patterns. In practice (i.e., without a known ground truth), an informed choice of the rotation criterion is not possible. From that perspective, enet regularized EFA has reasonable all-round properties without the necessity of additional subjective choices.
Although these results were convincing in favor of regularization, it should be acknowledged that they were obtained from unrealistically large samples. This is especially relevant for regularized EFA because the penalty is considered directly in the estimation step, and estimation performance largely depends on sample size. In order to conclude that regularization is a suitable alternative to factor rotation, it needs to be established whether these results also hold for more practically common samples sizes. Apart from that, an open question is how factor rotation and regularization compare with respect to the stability (i.e., standard errors) of the parameter estimates. We addressed these questions in the second simulation study.

STUDY 2: SMALL-SAMPLE PERFORMANCE
In this simulation, we compared the performances of factor rotation and regularization for more realistic sample sizes (N ¼ 100, N ¼ 200). In order to achieve a feasible simulation time, we focused on the basic conditions from Schmitt and Sass (2011) and the extended conditions with partial simple structure. These conditions were chosen in order to investigate if the performance advantages of regularized EFA from Study 1 are preserved in realistic samples. In addition, we only used the two most successful regularization methods (Enet Λ & Enet Λ;Φ ) from Study 1. These were compared with both versions of the Geomin rotation and CF-Varimax rotation. We evaluated the average recovery of the factor loading pattern and their empirical standard errors across samples.

Simulation model
The same simulation model was used as in the first simulation, but we reduced the number of conditions to keep the simulation feasible. Specifically, we included only the basic conditions and the extended conditions with partial simple structure (i.e., Extended Simple 1 and Extended Complex 1; Table 2).

Procedure
The procedure in Study 2 differed from Study 1 in two aspects: First, instead of simulating one large sample, we drew N rep ¼ 500 random samples from a multivariate normal distribution with either N ¼ 100 or N ¼ 200 observations per sample. Second, we only applied Geomin ( ¼ 0:01 or 0:5) and CF-Varimax rotated ML-EFA and enet penalized regularized EFA. For the penalty weight α, a grid of 10 values starting from α ¼ 0:001 with a step size of 10 À4 was used. For the enet weight β, we tested values between 0 and 1 with a step size of 0.1, that is, enet was allowed to result in pure lasso or ridge penalties if this optimized the sample BIC.

Dependent measures
For each sample, sign and order indeterminacies of EFA estimates were taken into account and the same dependent measures as in Study 1 were calculated. We report means and standard deviations across all samples for all measures.

Results
The main results of the simulation are summarized in Table 4. The average estimated factor loadings and factor correlations for each condition are available from the OSF. ML-EFA converged normally in all samples. Regularized EFA converged normally in all but 2 samples with N ¼ 100 in which Enet Λ penalized solutions were improper (factor correlations > 1). For these samples, the solution with the tuning parameters (α; β) entered the results that yielded the next best BIC. Overall, the model-implied covariance matrices fit the observed covariance matrices of the data very well as indicated by the SRMR. The fit was slightly worse for smaller samples and for regularized EFA compared to rotated EFA. The congruencies and biases of the main loadings indicated that all investigated methods recovered the general factor pattern sufficiently. Factor recovery and stability were marginally better for N ¼ 200 than for N ¼ 100. For the factor loadings, the stability of the estimates did not differ substantially between rotated and regularized solutions. For the factor correlations, Geomin ð ¼ 0:5Þ and CF-Varimax yielded smaller standard errors than the other methods. Notably, unlike in Study 1, none of the tested methods was able to perfectly recover the Perfect Simple pattern but rather yielded underestimated factor correlations and spurious cross-loadings. Across all conditions and methods, the biases of factor correlations and cross-loadings were almost perfectly correlated, r spearman ¼ À0:98 .
Both Geomin and CF-Varimax rotations resulted in underestimated cross-loadings and inflated cross-loadings for complex factor loading patterns. As in Study 1, Geomin ð ¼ 0:5Þ and CF-Varimax perfectly recovered the Approximate Simple pattern whereas Geomin ð ¼ 0:01Þ yielded similar (but slightly weaker) distortions as in the Complex condition. The enet regularization behaved very similar to the Geomin rotations. The inclusion of the factor correlation in the penalty term was more influential than in Study 1. Specifically, Enet Λ;Φ like Geomin ð ¼ 0:5Þ yielded less distorted estimates than Enet Λ in the presence of cross-loadings (Approximate Simple and Complex) but performed worse for Perfect Simple patterns. In the extended conditions, regularized EFA was superior to rotated EFA as indicated by (nearly) unbiased estimates (especially for the factor correlations). Altogether, the performance of factor rotation and regularization were very similar with respect to both accuracy and stability of the factor solutions.

Discussion
In this simulation, we investigated whether regularization yields comparable results as factor rotation for realistic sample sizes. Overall, the results of Study 1 generalized quite well to the small sample case and all investigated methods recovered the factor loading pattern with reasonable accuracy (Tucker's congruencies > 0.95; Lorenzo-Seva & Ten Berge, 2006) and yielded reasonable fit (SRMR < 0.05; Bentler & Yuan, 1999). In contrast to Study 1, none of the investigated methods was able to recover Perfect Simple patterns without biases. This was almost negligible for Geomin ð ¼ 0:01Þ but all other methods underestimated the factor correlations to  Note. The index symbols indicate which parameters were penalized with Λ = Factor loading matrix and Φ = Factor correlation matrix. The numbers outside the parentheses give the averages across all samples of the respective measure; the numbers in parentheses give the standard deviation across all samples of the respective measure.

586
SCHARF AND NESTLER a considerable extent. As in study 1, only regularized EFA was able to recover the parameters in the extended conditions.
Despite additional small sample biases, the present results support the notion that (enet) regularization may be used to estimate EFA without a separate rotation step. Differences between Geomin rotated and enet regularized estimates with respect to accuracy and stability were rather small. As in Study 1, both methods tended to oversimplify loading patterns with high cross-loadings, resulting in inflated factor correlations. We note that the general size of the biases was rather small, especially for the factor loadings, and would arguably not lead to fundamentally different interpretations of the factor loading pattern (e.g., a different selection of items). With respect to factor correlation bias, we note that interpretations of factor correlations should acknowledge the trade-off between the size of the cross-loadings and the size of the factor correlation that is inherent to all oblique factor analysis methods. In sum, we conclude regularization is a suitable alternative to the traditional factor rotation approach even in the case of small samples.

GENERAL DISCUSSION
In two simulation studies, we compared the performance of factor rotation and regularization in recovering predefined factor loading patterns that resembled typical psychometric data sets. Both with respect to asymptotic (Study 1) and finite sample performance (Study 2), regularization resulted in similar estimates as factor rotation. In line with previous notions (Morin et al., 2013), Geomin rotation with an increased rotation parameter ¼ 0:5 was superior for factor loading patterns with substantial cross-loadings but was unable to recover perfect simple structures. An enet penalty on both factor loadings and factor correlations showed the best overall performance among the investigated regularization methods and provided reasonable balance between the ability to recover simple structure, if it exists, and the ability to handle complex loading patterns.
We set out to investigate whether regularization is a suitable alternative to simple structure rotation in EFA. With respect to the estimation performance, our results confirm that regularization is a viable alternative to factor rotation. This is not to say that regularization is always the better approach to estimate the EFA parameters but, across all conditions in our simulations, it performed very well compared to the arguably best among the commonly used rotation techniques (cf. Schmitt & Sass, 2011). In addition, regularization was able to recover the EFA parameters when only a subset of the items followed a simple structurewhere factor rotation failed to do so (cf. Extended Simple/Complex 1 in Study 1). Thus, for typical psychometric data sets it can be expected that the results of regularized EFA match the results of simple structure rotated EFA very well.
Despite the promising performance of regularized EFA, the method is not without limitations: First, the penalization approach implemented in regSEM is very simplistic because all penalized parameters are treated equally, no matter whether they are factor loadings, structural coefficients or factor (co-) variances (or correlations). This could be a drawback because the number of factor loadings always exceeds the number of factor correlations or structural parameters; hence, the factor loadings will have a stronger relative influence on the estimation process (cf. Jacobucci et al., 2016;Jacobucci, 2017, for technical details). Second, the lasso and enet penalty (just as factor rotation) oversimplified complex loading patterns. It is well-known that lasso can result in overshrinkage of parameter estimates (here: cross-loadings), therefore, alternative methods for obtaining the final parameters may be considered (see Jacobucci et al., 2016, for a discussion). Lastly, our simulations were limited to the cases where all assumptions of ML-EFA and regSEM were fulfilled. Further research is needed to assess the sensitivity of regSEM to violations of distributional assumptions and model misspecification (e.g., correlated residuals). Especially the use of Likert scales may have profound consequences on the estimation performance (DiStefano, 2002;DiStefano & Morgan, 2014). Moreover, in the light of the increasing availability of large data sets (e.g., Kosinski, Wang, Lakkaraju, & Leskovec, 2016), future investigations should also consider conditions in which the number of variables exceeds the number of observations. Taken together, the present and previous research (e.g., Trendafilov, 2014;Yamamoto et al., 2017) suggest that, regarding estimation performance, there is not much to lose when replacing factor rotation with regularization but potential gains are also rather small (except for situations as in our extended conditions with partial simple structure). Therefore, applied readers may wonder why they should consider replacing wellestablished rotation procedures with regularization. We think that regularization has advantages both at a pragmatic and at a conceptual level: First, from a pragmatic perspective, a general use of regularization over factor rotation considerably alleviates the subjectivity of the analyses. Admittedly, researchers must still choose a penalty function just as they have to choose a rotation techniquewhich cannot be optimally done without a known ground truthbut the elastic net may provide a reasonable default choice and at least the tuning parameter can be determined in an objective and (nearly) automatized way.
At a conceptual level, regularization offers a generalization that subsumes EFA and confirmatory factor analysis (CFA). These methods are often considered separate methods that have different purposes. Using the penalized likelihood approach, however, both are simply two extremes in the space of possible models -differing only in which parameters enter the penalty term PðθÞ in Equation 8. EFA, on the one hand, has measurement models where all paths are allowed and all paths are considered to-be-regularized. CFA, on the other hand, has measurement models where only the theoretically motivated paths are estimated, all other paths are constrained to zero, and no parameter enters the penalty term. As a generalized method that contains both EFA and CFA as special cases, regularization allows all possible variations in between these two extremes. In particular, it has been considered to specify the theoretically motivated main loading paths but to not include them into the penalty term and only regularize all cross-loadings. Such a semi-confirmatory approach allows researchers to specify a model that considers there prior beliefs about the factor structure but that does not completely depend on the validity of these beliefs (Huang et al., 2017b).
Beyond basic factor analysis applications, the present findings have implications for the more generalized Exploratory Structural Equation modeling (ESEM) approach in which factor rotation is an important step as well ). ESEMs extent EFA by a structural model of the latent variables. As we operationalized regularized EFA from a regSEM perspective, regularized EFA can be easily extended to a regularized ESEM by (enet) penalizing the factor loadings and correlations of the exploratory factors. Such a regularized ESEM may be used to address similar research questions as ESEM (see also Huang et al., 2017a). Our results on factor rotation and regularization in the context of EFA should in principle apply for a comparison of ESEM and regularized ESEM as well. That is, the more complex the exploratory measurement model, the more rotation and regularization will underestimate the cross-loadings and overestimate the factor correlations. Consequently, the structural parameters may be distorted as well (Mai, Zhang, & Wen, 2018). Nevertheless, future research should investigate this notion empirically and investigate how additional regularization of the structural parameters affects the solutions.
Considering this connection, some of the central limitations of ESEMs may be solved in the regSEM framework. For instance, in ESEM, it is mandatory that the relationships between all or none of the exploratory factors with a given predictor are specified (otherwise, the rotation procedure is not valid, see . RegSEM is much more flexible in that respect, allowing researchers to choose which parameters of the model should be penalized and also add restrictions such as equality constraints to the model. Consequently, regSEM obviates the need for workarounds to ESEM limitations such as the ESEM within Confirmatory Factor Analysis approach (Morin et al., 2013). In this context, it would, for instance, be interesting to extend regSEM to the multi-level framework (e.g., Asparouhov & Muthén, 2012;Rabe-Hesketh, Skrondal, & Zheng, 2007). All in all, regSEM offers a consistent translation of the rotation problem into an estimation problemallowing for a unified framework of both confirmatory and exploratory techniques.
Future research should directly compare the performance of regSEM with other methods of semi-automatic model specification such as specification search in order to develop recommendations which method is appropriate for which purposes. In this context, it should be noted that regSEM and specification search are closely related to Bayesian Structural Equation Modeling (BSEM, e.g., Muthén & Asparouhov, 2011) because best subset selection and common penalty functions may be seen as a special case of BSEM with specific priors (Hastie et al., 2009;Jacobucci & Grimm, 2018). Apart from providing a unifying theoretical framework, the connection to BSEM offers a range of possibilities for improvements because different priors (i.e., penalty functions) could be placed on different parameters. For instance, one could use different priors on the factor loadings and factor correlations, respectively, in order to achieve a better balance between cross-loadings and factor correlations for complex loading patterns.

CONCLUSION
Regularization is an estimation method for complex statistical models with increasing popularity among social scientists. In two simulation studies, we compared the estimates of the traditional rotated EFA approach and regularized EFA for realistic factor loading patterns with varying complexity. Regularized EFA performed very similar to common factor rotation techniques in the majority of the considered conditions, indicating that regularization is a suitable alternative to the traditional rotation approach. Although regularized EFA was not unequivocally the best method across all conditions, the increased objectivity and the relation of the underlying regSEM to wider statistical frameworks such as ESEM and BSEM make it a valuable tool to be considered by social scientists.

ACKNOWLEDGMENTS
We thank Jana Pförtner for her assistance in preparing the simulation study and for proof-reading an earlier draft of the manuscript. We also thank Richard Rau for his valuable comments on an earlier draft of this manuscript.

SUPPLEMENTARY MATERIAL
Supplemental data for this article can be accessed here.  Note. The index symbols indicate which parameters were penalized with Λ = Factor loading matrix, and Φ = Factor correlation matrix.