This vignette provides an overview for the functionalities of the EFAtools package. The general aim of the package is to provide flexible implementations of different algorithms for an exploratory factor analyses (EFA) procedure, including factor retention methods, factor extraction and rotation methods, as well as the computation of a Schmid-Leiman solution and McDonaldβs omega coefficients.

The package was first designed to enable a comparison of EFA (specifically, principal axis factoring with subsequent promax rotation) performed in R using the **psych** package and EFA performed in SPSS. That is why some functions allow the specification of a type, including `"psych"`

and `"SPSS"`

, such that the respective procedure will be executed to match the output of these implementations (which do not always lead to the same results; see separate vignette **Replicate_SPSS_psych** for a demonstration of the replication of original results). This vignette will go through a complete example, that is, we will first show how to determine the number of factors to retain, then perform different factor extraction methods, run a Schmid-Leiman transformation and compute omegas.

The package can be installed from CRAN using `install.packages("EFAtools")`

, or from GitHub using `devtools::install_github("mdsteiner/EFAtools")`

, and then loaded using:

`library(EFAtools)`

In this vignette, we will use the `DOSPERT_raw`

dataset, which contains responses to the Domain Specific Risk Taking Scale (DOSPERT) of 3123 participants. The dataset is contained in the `EFAtools`

package, for details, see `?DOSPERT_raw`

. Note that this vignette is to provide a general overview and it is beyond its scope to explain all methods and functions in detail. If you want to learn more on the details and methods, please see the respective help functions for explanations and literature references. However, the dataset is rather large, so, just to save time when building the vignette, we will only use the first 500 observations. When you normally do your analyses, you use the full dataset.

```
# only use a subset to make analyses faster
<- DOSPERT_raw[1:500,] DOSPERT_sub
```

The first step in an EFA procedure is to test whether your data is suitable for factor analysis. To this end, the `EFAtools`

package provides the `BARTLETT()`

and the `KMO()`

functions. The Bartlettβs test of sphericity tests whether a correlation matrix is significantly different from an identity matrix (a correlation matrix with zero correlations between all variables). This test should thus be significant. The Kaiser-Meyer-Olkin criterion (KMO) represents the degree to which each observed variable is predicted by the other variables in the dataset and thus is another indicator for how correlated the different variables are.

We can test whether our `DOSPERT_sub`

dataset is suitable for factor analysis as follows.

```
# Bartlett's test of sphericity
BARTLETT(DOSPERT_sub)
#> βΉ 'x' was not a correlation matrix. Correlations are found from entered raw data.
#>
#> β The Bartlett's test of sphericity was significant at an alpha level of .05.
#> These data are probably suitable for factor analysis.
#>
#> πΒ²(435) = 5843.05, p < .001
# KMO criterion
KMO(DOSPERT_sub)
#> βΉ 'x' was not a correlation matrix. Correlations are found from entered raw data.
#>
#> ββ Kaiser-Meyer-Olkin criterion (KMO) ββββββββββββββββββββββββββββββββββββββββββ
#>
#> β The overall KMO value for your data is meritorious.
#> These data are probably suitable for factor analysis.
#>
#> Overall: 0.87
#>
#> For each variable:
#> ethR_1 ethR_2 ethR_3 ethR_4 ethR_5 ethR_6 finR_1 finR_2 finR_3 finR_4 finR_5
#> 0.908 0.927 0.906 0.818 0.926 0.837 0.880 0.816 0.853 0.857 0.853
#> finR_6 heaR_1 heaR_2 heaR_3 heaR_4 heaR_5 heaR_6 recR_1 recR_2 recR_3 recR_4
#> 0.861 0.912 0.898 0.892 0.889 0.882 0.935 0.888 0.927 0.905 0.849
#> recR_5 recR_6 socR_1 socR_2 socR_3 socR_4 socR_5 socR_6
#> 0.840 0.920 0.714 0.784 0.735 0.791 0.840 0.786
```

Note that these tests can also be run in the `N_FACTORS()`

function.

As the goal of EFA is to determine the underlying factors from a set of multiple variables, one of the most important decisions is how many factors can or should be extracted. There exists a plethora of factor retention methods to use for this decision. The problem is that there is no method that consistently outperforms all other methods. Rather, which factor retention method to use depends on the structure of the data: are there few or many indicators, are factors strong or weak, are the factor intercorrelations weak or strong. For rules on which methods to use, see, for example, Auerswald and Moshagen, (2019).

There are multiple factor retention methods implemented in the `EFAtools`

package. They can either be called with separate functions, or all (or a selection) of them using the `N_FACTORS()`

function.

Letβs first look at how to determine the number of factors to retain by calling separate functions. For example, if you would like to perform a parallel analysis based on squared multiple correlations (SMC; sometimes also called a parallel analysis with principal factors), you can do the following:

```
# determine the number of factors to retain using parallel analysis
PARALLEL(DOSPERT_sub, eigen_type = "SMC")
#> βΉ 'x' was not a correlation matrix. Correlations are found from entered raw data.
#> Parallel Analysis performed using 1000 simulated random data sets
#> Eigenvalues were found using SMC
#>
#> Decision rule used: means
#>
#> ββ Number of factors to retain according to ββββββββββββββββββββββββββββββββββββ
#>
#> β SMC-determined eigenvalues: 10
```

Generating the plot can also be suppressed if the output is printed explicitly:

```
# determine the number of factors to retain using parallel analysis
print(PARALLEL(DOSPERT_sub, eigen_type = "SMC"), plot = FALSE)
#> βΉ 'x' was not a correlation matrix. Correlations are found from entered raw data.
#> Parallel Analysis performed using 1000 simulated random data sets
#> Eigenvalues were found using SMC
#>
#> Decision rule used: means
#>
#> ββ Number of factors to retain according to ββββββββββββββββββββββββββββββββββββ
#>
#> β SMC-determined eigenvalues: 10
```

Other factor retention methods can be used accordingly. For example, to use the empirical Kaiser criterion, use the `EKC`

function:

```
# determine the number of factors to retain using parallel analysis
print(EKC(DOSPERT_sub), plot = FALSE)
#> βΉ 'x' was not a correlation matrix. Correlations are found from entered raw data.
#>
#> Empirical Kaiser criterion suggests 4 factors.
```

The following factor retention methods are currently implemented: comparison data (`CD()`

), empirical Kaiser criterion (`EKC()`

), the hull method (`HULL()`

), the Kaiser-Guttman criterion (`KGC()`

), parallel analysis (`PARALLEL()`

), scree test (`SCREE()`

), and sequential model tests (`SMT()`

). Many of these functions have multiple versions of the respective factor retention method implemented, for example, the parallel analysis can be done based on eigenvalues found using unity (principal components) or SMCs, or on an EFA procedure. Another example is the hull method, which can be used with different fitting methods (principal axis factoring [PAF], maximum likelihood [ML], or unweighted least squares [ULS]), and different goodness of fit indices. Please see the respective function documentations for details.

`N_FACTORS()`

If you want to use multiple factor retention methods, for example, to compare whether different methods suggest the same number of factors, it is easier to use the `N_FACTORS()`

function. This is a wrapper around all the implemented factor retention methods. Moreover, it also enables to run the Bartlettβs test of sphericity and compute the KMO criterion.

For example, to test the suitability of the data for factor analysis and to determine the number of factors to retain based on parallel analysis (but only using eigen values based on SMCs and PCA), the EKC, and the sequential model test, we can run the following code:

```
N_FACTORS(DOSPERT_sub, criteria = c("PARALLEL", "EKC", "SMT"),
eigen_type_other = c("SMC", "PCA"))
#>
π β― β― Running EKC
β π β― Running PARALLEL
β β π Running SMT
!
β β β Done#>
#> ββ Tests for the suitability of the data for factor analysis βββββββββββββββββββ
#>
#> Bartlett's test of sphericity
#>
#> β The Bartlett's test of sphericity was significant at an alpha level of .05.
#> These data are probably suitable for factor analysis.
#>
#> πΒ²(435) = 5843.05, p < .001
#>
#> Kaiser-Meyer-Olkin criterion (KMO)
#>
#> β The overall KMO value for your data is meritorious with 0.87.
#> These data are probably suitable for factor analysis.
#>
#> ββ Number of factors suggested by the different factor retention criteria ββββββ
#>
#> β Empirical Kaiser criterion: 4
#> β Parallel analysis with PCA: 5
#> β Parallel analysis with SMC: 10
#> β Sequential πΒ² model tests: 14
#> β Lower bound of RMSEA 90% confidence interval: 6
#> β Akaike Information Criterion: 13
```

If all possible factor retention methods should be used, it is sufficient to provide the data object (note that this takes a while, as the comparison data is computationally expensive and therefore relatively slow method, especially if larger datasets are used). We additionally specify the method argument to use unweighted least squares (ULS) estimation. This is a bit faster than using principle axis factoring (PAF) and it enables the computation of more goodness of fit indices:

```
N_FACTORS(DOSPERT_sub, method = "ULS")
#>
π β― β― β― β― β― β― Running CD
β π β― β― β― β― β― Running EKC
β β π β― β― β― β― Running HULL
β β β π β― β― β― Running KGC
β β β β π β― β― Running PARALLEL
β β β β β π β― Running SCREE
β β β β β β π Running SMT
!
β β β β β β β Done#>
#> ββ Tests for the suitability of the data for factor analysis βββββββββββββββββββ
#>
#> Bartlett's test of sphericity
#>
#> β The Bartlett's test of sphericity was significant at an alpha level of .05.
#> These data are probably suitable for factor analysis.
#>
#> πΒ²(435) = 5843.05, p < .001
#>
#> Kaiser-Meyer-Olkin criterion (KMO)
#>
#> β The overall KMO value for your data is meritorious with 0.87.
#> These data are probably suitable for factor analysis.
#>
#> ββ Number of factors suggested by the different factor retention criteria ββββββ
#>
#> β Comparison data: 5
#> β Empirical Kaiser criterion: 4
#> β Hull method with CAF: 6
#> β Hull method with CFI: 6
#> β Hull method with RMSEA: 6
#> β Kaiser-Guttman criterion with PCA: 8
#> β Kaiser-Guttman criterion with SMC: 4
#> β Kaiser-Guttman criterion with EFA: 3
#> β Parallel analysis with PCA: 5
#> β Parallel analysis with SMC: 10
#> β Parallel analysis with EFA: 7
#> β Sequential πΒ² model tests: 14
#> β Lower bound of RMSEA 90% confidence interval: 6
#> β Akaike Information Criterion: 13
```

Now, this is not the scenario one is happy about, but it still does happen: There is no obvious convergence between the methods and thus the choice of the number of factors to retain becomes rather difficult (and to some extend arbitrary). We will proceed with 6 factors, as it is what is typically used with DOSPERT data, but this does not mean that other number of factors are not just as plausible.

Note that all factor retention methods, except comparison data (CD), can also be used with correlation matrices. We use `method = "ULS"`

and `eigen_type_other = c("SMC", "PCA")`

to skip the slower criteria. In this case, the sample size has to be specified:

```
N_FACTORS(test_models$baseline$cormat, N = 500,
method = "ULS", eigen_type_other = c("SMC", "PCA"))
#> Warning in N_FACTORS(test_models$baseline$cormat, N = 500, method = "ULS", : ! 'x' was a correlation matrix but CD needs raw data. Skipping CD.
#>
β π β― β― β― β― β― Running EKC
β β π β― β― β― β― Running HULL
β β β π β― β― β― Running KGC
β β β β π β― β― Running PARALLEL
β β β β β π β― Running SCREE
β β β β β β π Running SMT
!
β β β β β β β Done#>
#> ββ Tests for the suitability of the data for factor analysis βββββββββββββββββββ
#>
#> Bartlett's test of sphericity
#>
#> β The Bartlett's test of sphericity was significant at an alpha level of .05.
#> These data are probably suitable for factor analysis.
#>
#> πΒ²(153) = 2173.28, p < .001
#>
#> Kaiser-Meyer-Olkin criterion (KMO)
#>
#> β The overall KMO value for your data is marvellous with 0.916.
#> These data are probably suitable for factor analysis.
#>
#> ββ Number of factors suggested by the different factor retention criteria ββββββ
#>
#> β Comparison data: NA
#> β Empirical Kaiser criterion: 2
#> β Hull method with CAF: 3
#> β Hull method with CFI: 3
#> β Hull method with RMSEA: 3
#> β Kaiser-Guttman criterion with PCA: 3
#> β Kaiser-Guttman criterion with SMC: 1
#> β Parallel analysis with PCA: 3
#> β Parallel analysis with SMC: 3
#> β Sequential πΒ² model tests: 3
#> β Lower bound of RMSEA 90% confidence interval: 2
#> β Akaike Information Criterion: 3
```

Multiple algorithms to perform an EFA and to rotate the found solutions are implemented in the `EFAtools`

package. All of them can be used using the `EFA()`

function. To perform the EFA, you can use one of principal axis factoring (PAF), maximum likelihood estimation (ML), and unweighted least squares (ULS; also sometimes referred to as MINRES). To rotate the solutions, the `EFAtools`

package offers varimax and promax rotations, as well as the orthogonal and oblique rotations provided by the `GPArotation`

package (i.e., the `GPArotation`

functions are called in the `EFA()`

function in this case).

You can run an EFA with PAF and no rotation like this:

```
EFA(DOSPERT_sub, n_factors = 6)
#> βΉ 'x' was not a correlation matrix. Correlations are found from entered raw data.
#>
#> EFA performed with type = 'EFAtools', method = 'PAF', and rotation = 'none'.
#>
#> ββ Unrotated Loadings ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F2 F3 F4 F5 F6
#> ethR_1 .558 -.193 .141 .232 -.129 .049
#> ethR_2 .472 -.157 .073 .242 .037 .059
#> ethR_3 .507 -.402 .164 .112 -.087 .062
#> ethR_4 .308 -.278 .189 .294 -.103 .224
#> ethR_5 .459 -.173 .202 .042 -.082 .045
#> ethR_6 .412 -.308 .230 .201 -.066 .207
#> finR_1 .602 -.304 -.461 -.045 .244 .028
#> finR_2 .294 .281 -.255 -.066 -.341 .086
#> finR_3 .626 -.285 -.471 -.031 .264 .114
#> finR_4 .557 -.012 -.357 -.023 -.437 -.063
#> finR_5 .657 -.342 -.491 -.064 .218 .028
#> finR_6 .551 .157 -.341 -.007 -.429 -.082
#> heaR_1 .464 -.079 .047 .133 .010 .046
#> heaR_2 .409 -.020 .163 .229 .040 .068
#> heaR_3 .497 -.101 .172 .127 .007 -.325
#> heaR_4 .547 -.110 .183 .038 .070 -.322
#> heaR_5 .382 .002 .175 .250 .028 .010
#> heaR_6 .504 .063 .102 .015 .006 -.196
#> recR_1 .443 .337 .057 -.064 .007 -.063
#> recR_2 .617 -.012 .147 -.141 .126 -.232
#> recR_3 .612 .178 .103 -.247 .080 -.153
#> recR_4 .634 .228 .254 -.447 .033 .235
#> recR_5 .640 .135 .267 -.463 .051 .261
#> recR_6 .629 .238 .135 -.244 -.060 -.014
#> socR_1 .075 .528 -.024 .233 .213 .124
#> socR_2 .358 .428 -.083 .247 .147 .013
#> socR_3 .118 .528 -.160 .122 .044 .116
#> socR_4 .235 .481 -.057 .214 .166 -.093
#> socR_5 .317 .330 .014 .154 -.030 .002
#> socR_6 .260 .462 -.058 .267 -.023 .082
#>
#> ββ Variances Accounted for βββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F2 F3 F4 F5 F6
#> SS loadings 7.013 2.417 1.529 1.246 0.850 0.641
#> Prop Tot Var 0.234 0.081 0.051 0.042 0.028 0.021
#> Cum Prop Tot Var 0.234 0.314 0.365 0.407 0.435 0.457
#> Prop Comm Var 0.512 0.177 0.112 0.091 0.062 0.047
#> Cum Prop Comm Var 0.512 0.689 0.800 0.891 0.953 1.000
#>
#> ββ Model Fit βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> CAF: .48
#> df: 270
```

To rotate the loadings (e.g., using a promax rotation) adapt the `rotation`

argument:

```
EFA(DOSPERT_sub, n_factors = 6, rotation = "promax")
#> βΉ 'x' was not a correlation matrix. Correlations are found from entered raw data.
#>
#> EFA performed with type = 'EFAtools', method = 'PAF', and rotation = 'promax'.
#>
#> ββ Rotated Loadings ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F3 F6 F1 F2 F4 F5
#> ethR_1 -.020 -.025 .547 .009 .130 .122
#> ethR_2 .132 -.055 .447 .115 .098 -.043
#> ethR_3 .068 .045 .532 -.212 .117 .030
#> ethR_4 -.039 -.034 .695 .000 -.162 -.002
#> ethR_5 -.051 .159 .382 -.100 .130 .036
#> ethR_6 -.003 .086 .664 -.063 -.089 -.044
#> finR_1 .852 -.019 -.015 .003 .022 .008
#> finR_2 -.044 .121 -.050 .120 -.211 .599
#> finR_3 .895 .035 .046 .061 -.094 -.011
#> finR_4 .098 -.072 .039 -.102 .035 .768
#> finR_5 .886 -.010 -.006 -.038 .023 .057
#> finR_6 .023 -.051 -.022 .031 .061 .774
#> heaR_1 .111 .042 .317 .092 .091 .019
#> heaR_2 -.003 .020 .416 .191 .093 -.081
#> heaR_3 -.071 -.141 .133 -.036 .674 -.009
#> heaR_4 -.006 -.032 .074 -.067 .697 -.074
#> heaR_5 -.053 -.037 .385 .196 .170 -.072
#> heaR_6 -.026 .061 .044 .059 .456 .048
#> recR_1 -.060 .253 -.067 .247 .219 .094
#> recR_2 .089 .230 -.046 -.046 .575 -.090
#> recR_3 .053 .401 -.155 .035 .424 .011
#> recR_4 -.010 .925 .042 .002 -.088 -.017
#> recR_5 .034 .945 .088 -.065 -.114 -.058
#> recR_6 -.067 .510 -.020 .053 .213 .154
#> socR_1 .025 .029 .016 .679 -.137 -.162
#> socR_2 .103 -.035 .057 .602 .072 -.020
#> socR_3 .029 .061 -.084 .559 -.199 .110
#> socR_4 .023 -.083 -.090 .583 .198 -.057
#> socR_5 -.086 .031 .097 .370 .077 .117
#> socR_6 -.072 -.035 .134 .568 -.079 .142
#>
#> ββ Factor Intercorrelations ββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F2 F3 F4 F5 F6
#> F1 1.000 0.348 0.455 0.046 0.502 0.463
#> F2 0.348 1.000 0.356 0.316 0.594 0.441
#> F3 0.455 0.356 1.000 0.031 0.567 0.329
#> F4 0.046 0.316 0.031 1.000 0.287 0.324
#> F5 0.502 0.594 0.567 0.287 1.000 0.487
#> F6 0.463 0.441 0.329 0.324 0.487 1.000
#>
#> ββ Variances Accounted for βββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F2 F3 F4 F5 F6
#> SS loadings 7.013 2.417 1.529 1.246 0.850 0.641
#> Prop Tot Var 0.234 0.081 0.051 0.042 0.028 0.021
#> Cum Prop Tot Var 0.234 0.314 0.365 0.407 0.435 0.457
#> Prop Comm Var 0.512 0.177 0.112 0.091 0.062 0.047
#> Cum Prop Comm Var 0.512 0.689 0.800 0.891 0.953 1.000
#>
#> ββ Model Fit βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> CAF: .48
#> df: 270
```

This now performed PAF with promax rotation with the specification, on average, we found to produce the most accurate results in a simulation analysis (see function documentation). If you want to replicate the implementation of the *psych* R package, you can set the `type`

argument to `"psych"`

:

```
EFA(DOSPERT_sub, n_factors = 6, rotation = "promax", type = "psych")
#> βΉ 'x' was not a correlation matrix. Correlations are found from entered raw data.
#>
#> EFA performed with type = 'psych', method = 'PAF', and rotation = 'promax'.
#>
#> ββ Rotated Loadings ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F2 F6 F4 F3 F5
#> ethR_1 .554 -.002 -.029 .116 -.024 .136
#> ethR_2 .441 .119 -.054 .088 .124 -.035
#> ethR_3 .556 -.227 .052 .113 .063 .048
#> ethR_4 .704 -.005 .003 -.195 -.042 -.006
#> ethR_5 .395 -.111 .156 .131 -.052 .044
#> ethR_6 .677 -.070 .117 -.109 -.007 -.047
#> finR_1 -.002 .017 -.018 .040 .827 .029
#> finR_2 -.034 .103 .116 -.229 -.041 .593
#> finR_3 .057 .079 .050 -.082 .868 -.000
#> finR_4 .078 -.134 -.100 .026 .096 .801
#> finR_5 .013 -.027 -.010 .042 .860 .081
#> finR_6 .005 .003 -.088 .051 .023 .801
#> heaR_1 .315 .093 .038 .087 .105 .024
#> heaR_2 .401 .197 .018 .082 -.007 -.083
#> heaR_3 .125 -.044 -.210 .700 -.072 .036
#> heaR_4 .068 -.074 -.103 .732 -.009 -.032
#> heaR_5 .367 .201 -.048 .162 -.056 -.068
#> heaR_6 .036 .054 .008 .477 -.028 .071
#> recR_1 -.086 .250 .215 .233 -.059 .089
#> recR_2 -.051 -.049 .170 .618 .085 -.065
#> recR_3 -.160 .032 .348 .464 .051 .019
#> recR_4 .051 -.002 .933 -.063 -.010 -.059
#> recR_5 .102 -.069 .961 -.087 .033 -.100
#> recR_6 -.018 .044 .476 .236 -.065 .147
#> socR_1 -.043 .713 .032 -.155 .022 -.200
#> socR_2 .007 .627 -.057 .062 .097 -.034
#> socR_3 -.123 .580 .061 -.217 .028 .079
#> socR_4 -.144 .608 -.121 .197 .020 -.065
#> socR_5 .070 .377 .010 .068 -.085 .109
#> socR_6 .094 .583 -.044 -.104 -.073 .122
#>
#> ββ Factor Intercorrelations ββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F2 F3 F4 F5 F6
#> F1 1.000 0.149 0.333 0.585 0.435 0.306
#> F2 0.149 1.000 0.398 0.387 0.083 0.436
#> F3 0.333 0.398 1.000 0.637 0.316 0.501
#> F4 0.585 0.387 0.637 1.000 0.463 0.483
#> F5 0.435 0.083 0.316 0.463 1.000 0.422
#> F6 0.306 0.436 0.501 0.483 0.422 1.000
#>
#> ββ Variances Accounted for βββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F2 F3 F4 F5 F6
#> SS loadings 7.013 2.417 1.529 1.246 0.850 0.641
#> Prop Tot Var 0.234 0.081 0.051 0.042 0.028 0.021
#> Cum Prop Tot Var 0.234 0.314 0.365 0.407 0.435 0.457
#> Prop Comm Var 0.512 0.177 0.112 0.091 0.062 0.047
#> Cum Prop Comm Var 0.512 0.689 0.800 0.891 0.953 1.000
#>
#> ββ Model Fit βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> CAF: .48
#> df: 270
```

If you want to use the *SPSS* implementation, you can set the `type`

argument to `"SPSS"`

:

```
EFA(DOSPERT_sub, n_factors = 6, rotation = "promax", type = "SPSS")
#> βΉ 'x' was not a correlation matrix. Correlations are found from entered raw data.
#>
#> EFA performed with type = 'SPSS', method = 'PAF', and rotation = 'promax'.
#>
#> ββ Rotated Loadings ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F4 F3 F2 F6 F5
#> ethR_1 .547 -.025 -.020 .009 .130 .122
#> ethR_2 .447 -.055 .132 .115 .098 -.043
#> ethR_3 .532 .045 .068 -.212 .117 .030
#> ethR_4 .695 -.035 -.039 .000 -.162 -.002
#> ethR_5 .382 .159 -.051 -.100 .130 .036
#> ethR_6 .664 .086 -.003 -.063 -.089 -.044
#> finR_1 -.015 -.019 .853 .003 .022 .008
#> finR_2 -.050 .121 -.044 .120 -.211 .599
#> finR_3 .046 .035 .895 .061 -.094 -.011
#> finR_4 .039 -.071 .098 -.102 .035 .768
#> finR_5 -.006 -.010 .887 -.038 .023 .057
#> finR_6 -.022 -.051 .023 .031 .061 .774
#> heaR_1 .318 .042 .111 .092 .091 .019
#> heaR_2 .416 .020 -.003 .191 .093 -.081
#> heaR_3 .133 -.140 -.071 -.036 .674 -.009
#> heaR_4 .074 -.031 -.006 -.067 .697 -.074
#> heaR_5 .385 -.037 -.053 .197 .170 -.072
#> heaR_6 .044 .062 -.026 .059 .455 .048
#> recR_1 -.067 .253 -.060 .247 .219 .094
#> recR_2 -.046 .231 .089 -.046 .574 -.090
#> recR_3 -.155 .402 .053 .035 .423 .011
#> recR_4 .042 .926 -.010 .002 -.089 -.017
#> recR_5 .088 .945 .034 -.065 -.114 -.058
#> recR_6 -.020 .510 -.067 .053 .212 .154
#> socR_1 .016 .029 .025 .679 -.137 -.162
#> socR_2 .057 -.035 .103 .602 .072 -.020
#> socR_3 -.084 .061 .029 .559 -.199 .110
#> socR_4 -.090 -.082 .023 .583 .197 -.057
#> socR_5 .097 .032 -.086 .370 .077 .117
#> socR_6 .134 -.035 -.072 .568 -.079 .142
#>
#> ββ Factor Intercorrelations ββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F2 F3 F4 F5 F6
#> F1 1.000 0.357 0.455 0.031 0.567 0.329
#> F2 0.357 1.000 0.349 0.316 0.594 0.441
#> F3 0.455 0.349 1.000 0.046 0.502 0.463
#> F4 0.031 0.316 0.046 1.000 0.287 0.324
#> F5 0.567 0.594 0.502 0.287 1.000 0.486
#> F6 0.329 0.441 0.463 0.324 0.486 1.000
#>
#> ββ Variances Accounted for βββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F2 F3 F4 F5 F6
#> SS loadings 7.013 2.417 1.529 1.246 0.850 0.641
#> Prop Tot Var 0.234 0.081 0.051 0.042 0.028 0.021
#> Cum Prop Tot Var 0.234 0.314 0.365 0.407 0.435 0.457
#> Prop Comm Var 0.512 0.177 0.112 0.091 0.062 0.047
#> Cum Prop Comm Var 0.512 0.689 0.800 0.891 0.953 1.000
#>
#> ββ Model Fit βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> CAF: .48
#> df: 270
```

This enables comparisons of different implementations. The `COMPARE()`

function provides an easy way to compare how similar two loading (pattern) matrices are:

```
COMPARE(
EFA(DOSPERT_sub, n_factors = 6, rotation = "promax", type = "psych")$rot_loadings,
EFA(DOSPERT_sub, n_factors = 6, rotation = "promax", type = "SPSS")$rot_loadings
)#> βΉ 'x' was not a correlation matrix. Correlations are found from entered raw data.
#>
#> βΉ 'x' was not a correlation matrix. Correlations are found from entered raw data.
#> Mean [min, max] absolute difference: 0.0159 [ 0.0001, 0.0712]
#> Median absolute difference: 0.0102
#> Max decimals where all numbers are equal: 0
#> Minimum number of decimals provided: 18
#>
#> F1 F2 F6 F4 F3 F5
#> ethR_1 0.0078 -0.0106 -0.0042 -0.0138 -0.0035 0.0148
#> ethR_2 -0.0059 0.0038 0.0008 -0.0100 -0.0078 0.0082
#> ethR_3 0.0245 -0.0146 0.0066 -0.0037 -0.0052 0.0186
#> ethR_4 0.0090 -0.0054 0.0375 -0.0325 -0.0028 -0.0045
#> ethR_5 0.0128 -0.0116 -0.0035 0.0007 -0.0011 0.0073
#> ethR_6 0.0129 -0.0067 0.0315 -0.0208 -0.0038 -0.0033
#> finR_1 0.0131 0.0143 0.0012 0.0183 -0.0257 0.0209
#> finR_2 0.0162 -0.0171 -0.0048 -0.0177 0.0025 -0.0055
#> finR_3 0.0110 0.0182 0.0150 0.0121 -0.0272 0.0107
#> finR_4 0.0392 -0.0319 -0.0286 -0.0093 -0.0023 0.0331
#> finR_5 0.0191 0.0113 0.0007 0.0192 -0.0266 0.0243
#> finR_6 0.0265 -0.0279 -0.0376 -0.0097 -0.0002 0.0271
#> heaR_1 -0.0024 0.0009 -0.0042 -0.0038 -0.0062 0.0053
#> heaR_2 -0.0153 0.0059 -0.0011 -0.0105 -0.0039 -0.0020
#> heaR_3 -0.0075 -0.0086 -0.0702 0.0258 -0.0012 0.0450
#> heaR_4 -0.0064 -0.0065 -0.0712 0.0356 -0.0028 0.0424
#> heaR_5 -0.0182 0.0048 -0.0112 -0.0087 -0.0024 0.0034
#> heaR_6 -0.0085 -0.0045 -0.0537 0.0222 -0.0015 0.0230
#> recR_1 -0.0191 0.0028 -0.0382 0.0146 0.0004 -0.0053
#> recR_2 -0.0045 -0.0025 -0.0615 0.0438 -0.0046 0.0251
#> recR_3 -0.0054 -0.0023 -0.0546 0.0413 -0.0024 0.0080
#> recR_4 0.0091 -0.0036 0.0078 0.0261 0.0001 -0.0421
#> recR_5 0.0148 -0.0042 0.0154 0.0270 -0.0012 -0.0417
#> recR_6 0.0021 -0.0090 -0.0342 0.0247 0.0011 -0.0068
#> socR_1 -0.0595 0.0345 0.0032 -0.0180 -0.0031 -0.0379
#> socR_2 -0.0494 0.0255 -0.0217 -0.0095 -0.0061 -0.0144
#> socR_3 -0.0390 0.0208 0.0001 -0.0185 -0.0017 -0.0311
#> socR_4 -0.0543 0.0254 -0.0381 -0.0005 -0.0032 -0.0083
#> socR_5 -0.0270 0.0072 -0.0213 -0.0092 0.0003 -0.0082
#> socR_6 -0.0401 0.0156 -0.0089 -0.0251 -0.0005 -0.0205
```

*Why would you want to do this?* One of us has had the experience that a reviewer asked whether the results can be reproduced in another statistical program than R. We therefore implemented this possibility in the package for an easy application of large scale, systematic comparisons.

Note that the `type`

argument of the `EFA()`

function only affects the implementations of principal axis factoring (PAF), varimax and promax rotations. The other procedures are not affected (except the order of the rotated factors for the other rotation methods).

As indicated previously, it is also possible to use different estimation and rotation methods. For example, to perform an EFA with ULS and an oblimin rotation, you can use the following code:

```
EFA(DOSPERT_sub, n_factors = 6, rotation = "oblimin", method = "ULS")
#> βΉ 'x' was not a correlation matrix. Correlations are found from entered raw data.
#>
#> EFA performed with type = 'EFAtools', method = 'ULS', and rotation = 'oblimin'.
#>
#> ββ Rotated Loadings ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F2 F3 F4 F5 F6
#> ethR_1 .039 .009 .522 -.007 .146 .137
#> ethR_2 .166 -.025 .419 .100 -.008 .097
#> ethR_3 .120 .064 .511 -.220 .052 .115
#> ethR_4 -.013 -.038 .644 -.031 .008 -.104
#> ethR_5 -.003 .177 .371 -.103 .059 .126
#> ethR_6 .029 .084 .620 -.086 -.025 -.048
#> finR_1 .850 -.013 -.012 .003 .028 .004
#> finR_2 -.022 .102 -.039 .106 .551 -.137
#> finR_3 .883 .028 .039 .054 .007 -.088
#> finR_4 .160 -.054 .063 -.109 .725 .068
#> finR_5 .890 -.005 -.002 -.037 .073 .006
#> finR_6 .085 -.028 .005 .023 .737 .092
#> heaR_1 .144 .064 .302 .083 .046 .091
#> heaR_2 .029 .047 .389 .176 -.043 .095
#> heaR_3 .011 -.044 .154 -.016 .051 .558
#> heaR_4 .073 .060 .101 -.041 -.008 .571
#> heaR_5 -.016 .002 .362 .185 -.031 .158
#> heaR_6 .035 .123 .064 .074 .092 .382
#> recR_1 -.024 .282 -.050 .254 .122 .193
#> recR_2 .153 .295 -.014 -.016 -.028 .467
#> recR_3 .105 .441 -.118 .062 .056 .348
#> recR_4 .010 .881 .053 .014 .001 -.060
#> recR_5 .052 .897 .095 -.053 -.039 -.084
#> recR_6 -.016 .524 .004 .068 .178 .191
#> socR_1 -.004 .033 -.010 .656 -.134 -.104
#> socR_2 .109 -.000 .043 .585 .018 .074
#> socR_3 .009 .054 -.095 .539 .110 -.147
#> socR_4 .031 -.034 -.091 .576 -.015 .169
#> socR_5 -.061 .059 .091 .359 .137 .086
#> socR_6 -.064 -.017 .114 .542 .155 -.036
#>
#> ββ Factor Intercorrelations ββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F2 F3 F4 F5 F6
#> F1 1.000 0.303 -0.379 -0.017 0.353 0.367
#> F2 0.303 1.000 -0.272 -0.245 0.353 0.412
#> F3 -0.379 -0.272 1.000 0.030 -0.186 -0.436
#> F4 -0.017 -0.245 0.030 1.000 -0.282 -0.173
#> F5 0.353 0.353 -0.186 -0.282 1.000 0.238
#> F6 0.367 0.412 -0.436 -0.173 0.238 1.000
#>
#> ββ Variances Accounted for βββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F2 F3 F4 F5 F6
#> SS loadings 7.013 2.417 1.529 1.246 0.850 0.641
#> Prop Tot Var 0.234 0.081 0.051 0.042 0.028 0.021
#> Cum Prop Tot Var 0.234 0.314 0.365 0.407 0.435 0.457
#> Prop Comm Var 0.512 0.176 0.112 0.091 0.062 0.047
#> Cum Prop Comm Var 0.512 0.689 0.800 0.891 0.953 1.000
#>
#> ββ Model Fit βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> πΒ²(270) = 213.54, p = .995
#> CFI = 1.00
#> RMSEA [90% CI] = .00 [.00; .00]
#> AIC = -326.46
#> BIC = -1464.40
#> CAF = .48
```

Of course, `COMPARE()`

can also be used to compare results from different estimation or rotation methods (in fact, to compare any two matrices), not just from different implementations:

```
COMPARE(
EFA(DOSPERT_sub, n_factors = 6, rotation = "promax")$rot_loadings,
EFA(DOSPERT_sub, n_factors = 6, rotation = "oblimin", method = "ULS")$rot_loadings,
x_labels = c("PAF and promax", "ULS and oblimin")
)#> βΉ 'x' was not a correlation matrix. Correlations are found from entered raw data.
#>
#> βΉ 'x' was not a correlation matrix. Correlations are found from entered raw data.
#> Mean [min, max] absolute difference: 0.0279 [ 0.0001, 0.1269]
#> Median absolute difference: 0.0226
#> Max decimals where all numbers are equal: 0
#> Minimum number of decimals provided: 18
#>
#> F3 F6 F1 F2 F4 F5
#> ethR_1 -0.0590 -0.0336 0.0243 0.0157 -0.0073 -0.0248
#> ethR_2 -0.0345 -0.0298 0.0274 0.0153 0.0006 -0.0348
#> ethR_3 -0.0522 -0.0194 0.0208 0.0073 0.0025 -0.0224
#> ethR_4 -0.0253 0.0038 0.0509 0.0306 -0.0581 -0.0099
#> ethR_5 -0.0482 -0.0176 0.0110 0.0028 0.0047 -0.0227
#> ethR_6 -0.0324 0.0021 0.0434 0.0224 -0.0410 -0.0190
#> finR_1 0.0021 -0.0062 -0.0028 -0.0001 0.0176 -0.0201
#> finR_2 -0.0214 0.0192 -0.0115 0.0142 -0.0742 0.0479
#> finR_3 0.0114 0.0078 0.0072 0.0069 -0.0059 -0.0179
#> finR_4 -0.0625 -0.0172 -0.0240 0.0071 -0.0334 0.0428
#> finR_5 -0.0033 -0.0055 -0.0047 -0.0006 0.0164 -0.0168
#> finR_6 -0.0624 -0.0233 -0.0267 0.0075 -0.0315 0.0379
#> heaR_1 -0.0331 -0.0224 0.0156 0.0095 0.0004 -0.0269
#> heaR_2 -0.0316 -0.0276 0.0272 0.0147 -0.0019 -0.0376
#> heaR_3 -0.0820 -0.0963 -0.0214 -0.0202 0.1159 -0.0596
#> heaR_4 -0.0792 -0.0922 -0.0270 -0.0266 0.1269 -0.0661
#> heaR_5 -0.0376 -0.0394 0.0226 0.0119 0.0124 -0.0404
#> heaR_6 -0.0610 -0.0626 -0.0203 -0.0150 0.0733 -0.0437
#> recR_1 -0.0354 -0.0295 -0.0172 -0.0074 0.0259 -0.0280
#> recR_2 -0.0637 -0.0644 -0.0328 -0.0302 0.1085 -0.0612
#> recR_3 -0.0520 -0.0402 -0.0371 -0.0277 0.0759 -0.0450
#> recR_4 -0.0201 0.0442 -0.0107 -0.0118 -0.0283 -0.0187
#> recR_5 -0.0177 0.0483 -0.0075 -0.0117 -0.0291 -0.0186
#> recR_6 -0.0510 -0.0141 -0.0249 -0.0149 0.0211 -0.0238
#> socR_1 0.0287 -0.0045 0.0260 0.0226 -0.0332 -0.0286
#> socR_2 -0.0068 -0.0353 0.0137 0.0166 -0.0018 -0.0377
#> socR_3 0.0203 0.0073 0.0111 0.0206 -0.0523 0.0005
#> socR_4 -0.0079 -0.0487 0.0005 0.0065 0.0286 -0.0417
#> socR_5 -0.0242 -0.0274 0.0057 0.0114 -0.0082 -0.0191
#> socR_6 -0.0076 -0.0179 0.0199 0.0257 -0.0427 -0.0123
```

Finally, if you are interested in factor scores from the EFA solution, these can be obtained with `FACTOR_SCORES()`

, a wrapper for `psych::factor.scores()`

to be used directly with an output from `EFA()`

:

```
<- EFA(DOSPERT_sub, n_factors = 6, rotation = "promax")
EFA_mod #> βΉ 'x' was not a correlation matrix. Correlations are found from entered raw data.
<- FACTOR_SCORES(DOSPERT_sub, f = EFA_mod) fac_scores
```

To improve performance of the iterative procedures (currently the parallel analysis, and the PAF, ML, and ULS methods) we implemented some of them in C++. For example, the following code compares the EFAtools parallel analysis with the corresponding one implemented in the psych package (the default of `PARALLEL()`

is to use 1000 datasets, but 25 is enough to show the difference):

```
::microbenchmark(
microbenchmarkPARALLEL(DOSPERT_sub, eigen_type = "SMC", n_datasets = 25),
::fa.parallel(DOSPERT_sub, SMC = TRUE, plot = FALSE, n.iter = 25)
psych
)#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Parallel analysis suggests that the number of factors = 10 and the number of components = 4
#> Unit: milliseconds
#> expr
#> PARALLEL(DOSPERT_sub, eigen_type = "SMC", n_datasets = 25)
#> psych::fa.parallel(DOSPERT_sub, SMC = TRUE, plot = FALSE, n.iter = 25)
#> min lq mean median uq max neval cld
#> 29.07872 34.86004 36.83956 36.60048 37.79203 59.74162 100 a
#> 179.96236 196.42759 209.89361 206.36946 223.96015 278.08575 100 b
```

Moreover, the following code compares the PAF implementation (of type βpsychβ) of the EFAtools package with the one from the psych package:

```
::microbenchmark(
microbenchmarkEFA(DOSPERT_raw, 6),
::fa(DOSPERT_raw, 6, rotate = "none", fm = "pa")
psych
)#> Unit: milliseconds
#> expr min lq
#> EFA(DOSPERT_raw, 6) 12.50138 12.89573
#> psych::fa(DOSPERT_raw, 6, rotate = "none", fm = "pa") 32.15971 34.42932
#> mean median uq max neval cld
#> 13.42817 13.12023 13.42090 18.44232 100 a
#> 37.88658 36.46675 38.50189 165.29582 100 b
```

While these differences are not large, they grow larger the more iterations the procedures need, which is usually the case if solutions are more tricky to find. Especially for simulations this might come in handy. For example, in one simulation analysis we ran over 10,000,000 EFAs, thus a difference of about 25 milliseconds per EFA leads to a difference in runtime of almost three days.

Instead of relying on one of the many possible implementations of, for example, PAF, and of using just one rotation (e.g., promax), it may be desirable to average different solutions to potentially arrive at a more robust, average solution. The `EFA_AVERAGE()`

function provides this possibility. In addition to the average solution it provides the variation across solutions, a matrix indicating the robustness of indicator-to-factor correspondences, and a visualisation of the average solution and the variability across solutions. For example, to average across all available factor extraction methods and across all available oblique rotations, the following code can be run:

```
# Average solution across many different EFAs with oblique rotations
<- EFA_AVERAGE(test_models$baseline$cormat, n_factors = 3, N = 500,
EFA_AV method = c("PAF", "ML", "ULS"), rotation = "oblique",
show_progress = FALSE)
# look at solution
EFA_AV#>
#> Averaging performed with averaging method mean (trim = 0) across 162 EFAs, varying the following settings: method, init_comm, criterion_type, start_method, rotation, k_promax, P_type, and varimax_type.
#>
#> The error rate is at 0%. Of the solutions that did not result in an error, 100% converged, 0% contained Heywood cases, and 100% were admissible.
#>
#>
#> ββ Indicator-to-Factor Correspondences βββββββββββββββββββββββββββββββββββββββββ
#>
#> For each cell, the proportion of solutions including the respective indicator-to-factor correspondence. A salience threshold of 0.3 was used to determine indicator-to-factor correspondences.
#>
#> F1 F2 F3
#> V1 .11 .00 1.00
#> V2 .11 .00 1.00
#> V3 .11 .00 .94
#> V4 .11 .00 1.00
#> V5 .11 .00 .94
#> V6 .11 .00 1.00
#> V7 .11 .94 .00
#> V8 .11 1.00 .00
#> V9 .11 .94 .00
#> V10 .11 1.00 .00
#> V11 .11 .89 .00
#> V12 .11 1.00 .00
#> V13 1.00 .00 .00
#> V14 1.00 .00 .00
#> V15 1.00 .00 .00
#> V16 1.00 .00 .00
#> V17 1.00 .00 .00
#> V18 1.00 .00 .00
#>
#>
#> ββ Loadings ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> ββ Mean ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F2 F3
#> V1 .025 .048 .576
#> V2 .060 .077 .451
#> V3 .115 .066 .425
#> V4 .157 .007 .518
#> V5 .198 -.002 .412
#> V6 .002 -.028 .658
#> V7 .074 .497 .102
#> V8 .056 .538 .046
#> V9 .100 .510 .018
#> V10 .048 .625 -.046
#> V11 .082 .336 .228
#> V12 .094 .606 .007
#> V13 .597 .083 -.047
#> V14 .531 -.056 .093
#> V15 .548 .122 -.049
#> V16 .540 -.041 .097
#> V17 .633 -.033 -.009
#> V18 .542 .009 .060
#>
#> ββ Range βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F2 F3
#> V1 0.513 0.086 0.239
#> V2 0.431 0.093 0.186
#> V3 0.394 0.108 0.179
#> V4 0.415 0.110 0.214
#> V5 0.315 0.122 0.177
#> V6 0.514 0.104 0.267
#> V7 0.527 0.255 0.089
#> V8 0.520 0.275 0.078
#> V9 0.470 0.276 0.080
#> V10 0.533 0.313 0.097
#> V11 0.482 0.176 0.102
#> V12 0.548 0.324 0.103
#> V13 0.081 0.289 0.114
#> V14 0.063 0.220 0.117
#> V15 0.091 0.280 0.107
#> V16 0.072 0.230 0.122
#> V17 0.108 0.270 0.124
#> V18 0.081 0.246 0.118
#>
#>
#> ββ Factor Intercorrelations from Oblique Solutions βββββββββββββββββββββββββββββ
#>
#> ββ Mean ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F2 F3
#> F1 1.000 0.431 0.518
#> F2 0.431 1.000 0.454
#> F3 0.518 0.454 1.000
#>
#> ββ Range βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F2 F3
#> F1 0.000 1.276 0.679
#> F2 1.276 0.000 1.316
#> F3 0.679 1.316 0.000
#>
#>
#> ββ Variances Accounted for βββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> ββ Mean ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F2 F3
#> SS loadings 2.443 1.929 1.904
#> Prop Tot Var 0.136 0.107 0.106
#> Prop Comm Var 0.389 0.307 0.303
#>
#> ββ Range βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F2 F3
#> SS loadings 2.831 1.356 1.291
#> Prop Tot Var 0.157 0.075 0.072
#> Prop Comm Var 0.419 0.215 0.215
#>
#>
#> ββ Model Fit βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> M (SD) [Min; Max]
#> πΒ²: 101.73 (34.62) [53.23; 125.98]
#> df: 102
#> p: .369 (.450) [.054; 1.000]
#> CFI: 1.00 (.00) [1.00; 1.00]
#> RMSEA: .01 (.01) [.00; .02]
#> AIC: -102.27 (34.62) [-150.77; -78.02]
#> BIC: -532.16 (34.62) [-580.66; -507.91]
#> CAF: .50 (.00) [.50; .50]
```

The first matrix of the output tells us that the indicators are mostly allocated to the same factors. However, that some rowsums are larger than one also tells as that there likely are some cross loadings present in some solutions. Moreover, the relatively high percentages of salient pattern coefficients all loading on the first factor may indicate that some rotation methods failed to achieve simple structure and it might be desirable to exclude these from the model averaging procedure. The rest of the output is similar to the normal `EFA()`

outputs shown above, only that in addition to the average coefficients their range is also shown. Finally, the plot shows the average pattern coefficients and their ranges.

**Important disclaimer:** While it is possible that this approach provides more robust results, we are unaware of simulation studies that have investigated and shown this. Therefore, it might make sense to for now use this approach mainly to test the robustness of the results obtained with one single EFA implementation.

For the Schmid-Leiman transformation and computation of omegas, we will use PAF and promax rotation:

```
<- EFA(DOSPERT_sub, n_factors = 6, rotation = "promax")
efa_dospert #> βΉ 'x' was not a correlation matrix. Correlations are found from entered raw data.
efa_dospert#>
#> EFA performed with type = 'EFAtools', method = 'PAF', and rotation = 'promax'.
#>
#> ββ Rotated Loadings ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F3 F6 F1 F2 F4 F5
#> ethR_1 -.020 -.025 .547 .009 .130 .122
#> ethR_2 .132 -.055 .447 .115 .098 -.043
#> ethR_3 .068 .045 .532 -.212 .117 .030
#> ethR_4 -.039 -.034 .695 .000 -.162 -.002
#> ethR_5 -.051 .159 .382 -.100 .130 .036
#> ethR_6 -.003 .086 .664 -.063 -.089 -.044
#> finR_1 .852 -.019 -.015 .003 .022 .008
#> finR_2 -.044 .121 -.050 .120 -.211 .599
#> finR_3 .895 .035 .046 .061 -.094 -.011
#> finR_4 .098 -.072 .039 -.102 .035 .768
#> finR_5 .886 -.010 -.006 -.038 .023 .057
#> finR_6 .023 -.051 -.022 .031 .061 .774
#> heaR_1 .111 .042 .317 .092 .091 .019
#> heaR_2 -.003 .020 .416 .191 .093 -.081
#> heaR_3 -.071 -.141 .133 -.036 .674 -.009
#> heaR_4 -.006 -.032 .074 -.067 .697 -.074
#> heaR_5 -.053 -.037 .385 .196 .170 -.072
#> heaR_6 -.026 .061 .044 .059 .456 .048
#> recR_1 -.060 .253 -.067 .247 .219 .094
#> recR_2 .089 .230 -.046 -.046 .575 -.090
#> recR_3 .053 .401 -.155 .035 .424 .011
#> recR_4 -.010 .925 .042 .002 -.088 -.017
#> recR_5 .034 .945 .088 -.065 -.114 -.058
#> recR_6 -.067 .510 -.020 .053 .213 .154
#> socR_1 .025 .029 .016 .679 -.137 -.162
#> socR_2 .103 -.035 .057 .602 .072 -.020
#> socR_3 .029 .061 -.084 .559 -.199 .110
#> socR_4 .023 -.083 -.090 .583 .198 -.057
#> socR_5 -.086 .031 .097 .370 .077 .117
#> socR_6 -.072 -.035 .134 .568 -.079 .142
#>
#> ββ Factor Intercorrelations ββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F2 F3 F4 F5 F6
#> F1 1.000 0.348 0.455 0.046 0.502 0.463
#> F2 0.348 1.000 0.356 0.316 0.594 0.441
#> F3 0.455 0.356 1.000 0.031 0.567 0.329
#> F4 0.046 0.316 0.031 1.000 0.287 0.324
#> F5 0.502 0.594 0.567 0.287 1.000 0.487
#> F6 0.463 0.441 0.329 0.324 0.487 1.000
#>
#> ββ Variances Accounted for βββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> F1 F2 F3 F4 F5 F6
#> SS loadings 7.013 2.417 1.529 1.246 0.850 0.641
#> Prop Tot Var 0.234 0.081 0.051 0.042 0.028 0.021
#> Cum Prop Tot Var 0.234 0.314 0.365 0.407 0.435 0.457
#> Prop Comm Var 0.512 0.177 0.112 0.091 0.062 0.047
#> Cum Prop Comm Var 0.512 0.689 0.800 0.891 0.953 1.000
#>
#> ββ Model Fit βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> CAF: .48
#> df: 270
```

The indicator names in the output (i.e., the rownames of the rotated loadings section) tell us which domain (out of ethical, financial, health, recreational, and social risks) an indicator stems from. From the pattern coefficients it can be seen that these theoretical domains are recovered relatively well in the six factor solution, that is, usually, the indicators from the same domain load onto the same factor. When we take a look at the factor intercorrelations, we can see that there are some strong and some weak correlations. It might be worthwhile to explore whether a general factor can be obtained, and which factors load more strongly on it. To this end, we will use a Schmid-Leiman (SL) transformation.

The SL transformation or orthogonalization transforms an oblique solution into a hierarchical, orthogonalized solution. To do this, the `EFAtools`

package provides the `SL()`

function.

```
<- SL(efa_dospert)
sl_dospert
sl_dospert#>
#> EFA for second-order loadings performed with type = 'EFAtools' and method = 'PAF'
#>
#> ββ Schmid-Leiman Solution ββββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> g F1 F2 F3 F4 F5 F6 h2 u2
#> ethR_1 .487 .440 .009 -.016 .067 .093 -.018 .445 .555
#> ethR_2 .400 .360 .110 .104 .050 -.033 -.041 .318 .682
#> ethR_3 .441 .428 -.202 .054 .060 .023 .034 .427 .573
#> ethR_4 .225 .559 .000 -.031 -.084 -.001 -.026 .372 .628
#> ethR_5 .406 .308 -.095 -.041 .067 .028 .119 .290 .710
#> ethR_6 .325 .534 -.060 -.003 -.046 -.034 .064 .402 .598
#> finR_1 .523 -.012 .003 .676 .011 .006 -.014 .731 .269
#> finR_2 .266 -.041 .114 -.035 -.109 .458 .090 .316 .684
#> finR_3 .527 .037 .058 .709 -.049 -.008 .026 .789 .211
#> finR_4 .529 .031 -.097 .078 .018 .587 -.053 .643 .357
#> finR_5 .574 -.005 -.036 .703 .012 .043 -.008 .827 .173
#> finR_6 .528 -.018 .029 .018 .031 .592 -.038 .633 .367
#> heaR_1 .403 .256 .088 .088 .047 .014 .031 .246 .754
#> heaR_2 .344 .335 .182 -.002 .048 -.062 .015 .270 .730
#> heaR_3 .503 .107 -.034 -.056 .347 -.007 -.105 .400 .600
#> heaR_4 .548 .060 -.064 -.005 .359 -.057 -.024 .441 .559
#> heaR_5 .331 .310 .187 -.042 .088 -.055 -.028 .254 .746
#> heaR_6 .491 .036 .056 -.021 .234 .037 .045 .304 .696
#> recR_1 .417 -.054 .235 -.047 .113 .072 .188 .288 .712
#> recR_2 .601 -.037 -.044 .071 .296 -.069 .171 .492 .508
#> recR_3 .590 -.125 .033 .042 .218 .008 .299 .503 .497
#> recR_4 .550 .034 .002 -.008 -.045 -.013 .689 .780 .220
#> recR_5 .549 .071 -.062 .027 -.058 -.044 .704 .812 .188
#> recR_6 .585 -.016 .050 -.053 .109 .118 .379 .518 .482
#> socR_1 .030 .013 .646 .020 -.071 -.124 .021 .439 .561
#> socR_2 .306 .046 .573 .081 .037 -.015 -.026 .433 .567
#> socR_3 .081 -.068 .532 .023 -.102 .084 .045 .315 .685
#> socR_4 .217 -.073 .555 .018 .102 -.044 -.062 .376 .624
#> socR_5 .282 .078 .352 -.068 .040 .090 .023 .225 .775
#> socR_6 .211 .108 .540 -.057 -.040 .109 -.026 .366 .634
#>
#>
#> ββ Variances Accounted for βββββββββββββββββββββββββββββββββββββββββββββββββββββ
#>
#> g F1 F2 F3 F4 F5 F6
#> SS loadings 5.711 1.550 1.994 1.519 0.552 1.000 1.326
#> Prop Tot Var 0.190 0.052 0.066 0.051 0.018 0.033 0.044
#> Cum Prop Tot Var 0.190 0.242 0.309 0.359 0.378 0.411 0.455
#> Prop Comm Var 0.418 0.114 0.146 0.111 0.040 0.073 0.097
#> Cum Prop Comm Var 0.418 0.532 0.678 0.789 0.830 0.903 1.000
```

From the output, it can be seen that all, except the social domain indicators substantially load on the general factor. That is, the other domains covary substantially.

Finally, we can compute omega reliability estimates based on the SL solution. To this end, we can either specify the variable-to-factor correspondences, or let them be determined automatically (in which case the highest factor loading will be taken, which might lead to a different solution than what is desired, in the presence of cross-loadings). Given that no cross-loadings are present here, it is easiest to let the function automatically determine the variable-to-factor correspondence. To this end, we will set the `type`

argument to `"psych"`

.

```
OMEGA(sl_dospert, type = "psych")
#> Omega total, omega hierarchical, and omega subscale for the general factor (top row) and the group factors:
#>
#> tot hier sub
#> g 0.922 0.718 0.177
#> F1 0.791 0.376 0.415
#> F2 0.739 0.124 0.615
#> F3 0.912 0.344 0.568
#> F4 0.709 0.532 0.177
#> F5 0.745 0.295 0.451
#> F6 0.866 0.474 0.393
```

If we wanted to specify the variable to factor correspondences explicitly (for example, according to theoretical expectations), we could do it in the following way:

```
OMEGA(sl_dospert, factor_corres = matrix(c(rep(0, 18), rep(1, 6), rep(0, 30),
rep(1, 6), rep(0, 6), 1, 0, 1, 0, 1,
rep(0, 19), rep(1, 6), rep(0, 31), 1, 0,
1, 0, 1, rep(0, 30), rep(1, 6),
rep(0, 12)), ncol = 6, byrow = FALSE))
#> Omega total, omega hierarchical, and omega subscale for the general factor (top row) and the group factors:
#>
#> tot hier sub
#> g 0.922 0.718 0.083
#> F1 0.536 0.536 0.001
#> F2 0.742 0.082 0.660
#> F3 0.912 0.344 0.568
#> F4 0.318 0.317 0.001
#> F5 0.745 0.295 0.451
#> F6 0.479 0.479 0.000
```