---
title: "EFAtools"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{EFAtools}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
message = FALSE,
fig.width = 7,
fig.align = "center"
)
```
`EFAtools` provides a complete workflow for exploratory factor analysis (EFA). Starting
from raw data or a correlation matrix, it takes you through every stage of a typical
analysis:
1. **Screening** the data for factorability, multicollinearity, non-normality, and
outliers (`efa_screen()`, with `efa_bartlett()` and `efa_kmo()` for the individual
measures).
2. **Retention**: deciding how many factors to keep, with a comprehensive suite of factor
retention criteria (`efa_retain()` and the individual criteria it wraps).
3. **Extraction and rotation** of the factor solution (`efa_fit()`), with a choice
of estimators and orthogonal or oblique rotations.
4. **Post-processing**: comparing solutions (`efa_compare()`), averaging across analytic
choices (`efa_average()`), estimating factor scores (`efa_scores()`), a Schmid-Leiman
transformation (`efa_schmid_leiman()`), and reliability coefficients including
McDonald's omegas (`efa_reliability()`).
The computationally intensive routines are implemented in C++ for speed. This vignette
walks through the whole workflow on a single data set; it is an overview rather than a
tutorial on the methods themselves, so please consult the individual help pages for the
statistical details and literature references.
The package can be installed from CRAN using `install.packages("EFAtools")`, or from GitHub
using `pak::pak("mdsteiner/EFAtools")`, and then loaded using:
```{r}
library(EFAtools)
```
We use the `DOSPERT_raw` data set throughout, which contains responses of 3123 participants
to the Domain Specific Risk Taking Scale (DOSPERT); see `?DOSPERT_raw` for details. The
data are rather large, so, purely to keep this vignette quick to build, we work with the
first 500 observations. In a real analysis you would of course use the full data set.
```{r}
# only use a subset to make analyses faster
DOSPERT_sub <- DOSPERT_raw[1:500, ]
```
## Screening the Data
Before extracting factors it is worth checking whether the data are suitable for factor
analysis, and whether they satisfy the assumptions of the estimator you plan to use.
`efa_screen()` collects the relevant diagnostics in one place. From raw data it reports the
Kaiser-Meyer-Olkin (KMO) measure of sampling adequacy (overall and per variable), Bartlett's
test of sphericity, the determinant and condition number of the correlation matrix, the
squared multiple correlation (SMC) and per-variable diagnostics for each item, tests of
multivariate normality (Mardia and Henze-Zirkler), and multivariate outliers, and it closes
with a set of recommendations.
```{r, warning = FALSE}
efa_screen(DOSPERT_sub, seed = 2)
```
KMO measures, on a 0-1 scale, how well each variable is predicted by the others, and
Bartlett's test checks that the correlations are not all zero. The overall KMO is
meritorious and Bartlett's test is clearly significant, so the data are factorable. The determinant looks small, but it is a product of 30 eigenvalues and falls as
variables are added; the condition index (6.8) is the measure that does not move with the
number of variables, and it is well below the value of 30 that flags a near linear
dependency.
More interesting are the last few
sections: three items have a sparse response category, and both the Mardia and the
Henze-Zirkler tests reject multivariate normality. A consequence is that normal-theory
standard errors and fit indices may be biased, so robust (sandwich) or bootstrapped standard
errors are preferable, and Bartlett's test should be read with caution. For a fully ordinal
treatment you would move to polychoric correlations with a diagonally weighted least-squares
(DWLS) estimator.
The outlier section reports something other than the usual case, and it says so. The
diagnostic normally measures each case against a high-breakdown location and scatter,
estimated from a minimum covariance determinant (MCD) subset found by a randomized search.
For these data the search keeps landing on a subset where one item is constant (several
DOSPERT items have a very common lowest response), so the covariance is singular and no valid
subset exists.
`efa_screen()` falls back to classical Mahalanobis distances. It names the reason and labels
the distances it reports. A classical covariance uses every observation, outliers included, so
the outliers inflate the scatter that measures them. The diagnostic is no longer
high-breakdown and it tends to under-flag. Read the 60 flagged cases as a lower bound, and as
a sign of the non-normality the tests above report, rather than as a list of cases to remove.
The MCD does not run, so its random subset search does not run either and `seed` changes
nothing here. `seed` matters only when the robust estimate is available, as in the ordinal
example in the [EFA with ordinal and missing data](Ordinal_and_missing_data.html) vignette.
Inspect the flagged cases (via `$outliers$flagged`) rather than deleting them automatically.
Bartlett's test and the KMO criterion, introduced above, are also available as standalone
functions if you only want the two classic factorability tests:
```{r}
# Bartlett's test of sphericity
efa_bartlett(DOSPERT_sub)
# KMO criterion
efa_kmo(DOSPERT_sub)
```
Both tests can also be requested directly inside `efa_retain()`, alongside the factor
retention criteria discussed next.
## Deciding How Many Factors to Retain
One of the most consequential decisions in an EFA is how many factors to extract. There is a
large number of factor retention criteria, and no single one dominates the others: which
criterion works best depends on the structure of the data — how many indicators there are,
how strong the factors are, and how strongly they intercorrelate. The current recommendation
is therefore to examine several criteria together and check whether they converge (see, for
example, [Auerswald and Moshagen, 2019](https://doi.org/10.1037/met0000200)).
`EFAtools` implements a broad set of these criteria. They can be called individually, or all
(or a selection) of them together through `efa_retain()`.
### Calling Individual Criteria
Each criterion has its own function. For a parallel analysis based on squared multiple
correlations (SMC; also called parallel analysis with principal factors), for example, you
would run:
```{r}
# parallel analysis based on SMC eigenvalues
pa <- efa_parallel(DOSPERT_sub, eigen_type = "smc")
pa
```
Printing the result reports the suggested number of factors. Each criterion also has a
`plot()` method that draws the figure behind the decision — for a parallel analysis, the
observed eigenvalues against the reference eigenvalues from the simulated data:
```{r}
plot(pa)
```
Other criteria work the same way. The empirical Kaiser criterion, for instance:
```{r}
# empirical Kaiser criterion
efa_ekc(DOSPERT_sub)
```
The following criteria are currently implemented: comparison data (`efa_cd()`), the empirical
Kaiser criterion (`efa_ekc()`), the hull method (`efa_hull()`), the Kaiser-Guttman criterion
(`efa_kgc()`), the minimum average partial test (`efa_map()`), the next eigenvalue sufficiency
test (`efa_nest()`), parallel analysis (`efa_parallel()`), the scree test (`efa_scree()`), and
sequential model tests (`efa_smt()`). Many of them offer several variants of the underlying
method. For instance, parallel analysis, can be based on eigenvalues from unity (principal
components), SMCs, or a full EFA. See the individual help pages for the options.
### Running Several Criteria at Once with `efa_retain()`
To compare criteria, it is easier to call `efa_retain()`, a wrapper around the implemented
retention criteria that also runs Bartlett's test and computes the KMO. You can pick a subset
with the `criteria` argument:
```{r}
ret <- efa_retain(DOSPERT_sub,
criteria = c("parallel", "ekc", "kgc", "smt", "map"))
ret
```
`plot()` on the result draws the figure for every criterion that has one:
```{r}
plot(ret)
```
Omitting `criteria` does not run every criterion listed above. The default is a subset of six:
comparison data, the empirical Kaiser criterion, the hull method, the minimum average partial
test, the next eigenvalue sufficiency test, and parallel analysis. The call above instead adds
the Kaiser-Guttman criterion and the sequential model tests, and drops comparison data (slow
on larger data sets), the hull method, and the next eigenvalue sufficiency test.
This is the situation one would rather avoid, but it does happen: the criteria disagree, with
suggestions ranging from four to well above ten factors. When there is no clear convergence,
the choice becomes partly a matter of judgement. We proceed with six factors, consistent with the minimum average partial test and the RMSEA-based sequential model test above (this version of the DOSPERT scale is usually
modeled with five factors; the sixth factor here splits the financial-risk items into two).
All criteria except comparison data (which needs raw data) can also be applied to a
correlation matrix, in which case the sample size must be supplied. On a cleaner data set the
criteria often agree:
```{r}
efa_retain(test_models$baseline$cormat, N = 500, estimator = "uls",
criteria = c("parallel", "ekc", "smt"),
eigen_type_other = c("smc", "pca"))
```
Here the criteria agree on three factors, with only the conservative RMSEA lower-bound
variant of the sequential model test suggesting two.
## Extracting Factors with `efa_fit()`
`efa_fit()` performs the factor extraction and rotation. For extraction you can choose
principal axis factoring (PAF), maximum likelihood (ML), or unweighted least squares (ULS,
also known as MINRES); diagonally weighted least squares (DWLS) is additionally available
for ordinal data analysed with polychoric or tetrachoric correlations (see the
[EFA with ordinal and missing data](Ordinal_and_missing_data.html) vignette). For rotation
the package offers varimax and promax as well as a range of further rotations (orthogonal
rotations keep the factors uncorrelated, oblique ones let them correlate): quartimax and
equamax are orthogonal; oblimin, quartimin, and simplimax are oblique; and geomin, bentler,
and bifactor come in both an orthogonal form (`geominT`, `bentlerT`, `bifactorT`) and an
oblique one (`geominQ`, `bentlerQ`, `bifactorQ`). The suffix is part of the name: `geomin` on
its own matches both forms and is rejected. All are computed by rotation engines built into
the package.
An EFA with PAF and no rotation is as simple as:
```{r}
efa_fit(DOSPERT_sub, n_factors = 6)
```
To rotate the loadings, set the `rotation` argument. Here we use a promax rotation. We assign the
result, as later sections reuse this solution:
```{r}
efa_dospert <- efa_fit(DOSPERT_sub, n_factors = 6, rotation = "promax")
efa_dospert
```
We can also use a different estimator and rotation, here for example ULS and an oblimin rotation:
```{r}
efa_uls <- efa_fit(DOSPERT_sub, n_factors = 6, rotation = "oblimin", estimator = "uls")
efa_uls
```
Note that `estimator = "uls"` was written in lowercase although the canonical spelling is
`"ULS"`. This holds throughout the `efa_*` interface: wherever an argument takes a fixed set
of values, capitalization does not matter, while the result stores and prints the canonical
spelling. For an argument that takes a single value an unambiguous abbreviation is enough as
well, and a value matching nothing raises an error naming the argument and listing the
choices. Arguments that take several values at once, such as the analytic choices
`efa_average()` averages over, need the full spelling of each value.
## Comparing Solutions with `efa_compare()`
`efa_compare()` provides a quick way to see how similar two loading (pattern) matrices are —
whether they come from different implementations, different estimators, or different rotations.
For instance, to compare the PAF/promax solution with the ULS/oblimin one:
```{r}
efa_compare(efa_dospert$rot_loadings, efa_uls$rot_loadings,
x_labels = c("PAF and promax", "ULS and oblimin"))
```
## Averaging Across Analytic Choices with `efa_average()`
Comparing two solutions by hand answers the question for those two. `efa_average()` asks it
across a whole grid at once: it runs an EFA for every combination of the settings you allow it
to vary, averages the loadings across the solutions that converged, and reports the range each
loading spans. That range is the useful part — it shows how much of the result is the data and
how much is the analytic choice. Here we vary only the estimator, keeping one implementation
each of PAF, ULS, and ML with a promax rotation:
```{r}
avg_dospert <- efa_average(DOSPERT_sub, n_factors = 6,
estimator = c("paf", "uls", "ml"),
type = "EFAtools", start_method = "psych")
avg_dospert
```
Left unrestricted, the defaults sweep the different implementations of principal axis
factoring and promax rotation instead, which is considerably more work. See `?efa_average`
for the settings that can be varied, and for `plot()`, which shows the minimum, maximum, and
average loading of every indicator.
## Factor Scores with `efa_scores()`
If you need factor scores, `efa_scores()` estimates them directly from an `efa_fit()` solution
and reports score-quality diagnostics. Printing the object gives the determinacy of each score
and the Guttman indeterminacy index. `summary()` adds the factor weights, the score
intercorrelations, and a score validity/univocality matrix. The off-diagonals of that matrix
are the univocality — how strongly each score correlates with the *other* factors:
```{r}
fac_scores <- efa_scores(DOSPERT_sub, f = efa_dospert)
fac_scores
summary(fac_scores)
```
## Schmid-Leiman Transformation and McDonald's Omegas
For the Schmid-Leiman transformation and the omega coefficients we return to the PAF/promax
solution (`efa_dospert`) from above. Its indicator names encode the risk domain each item
belongs to (ethical, financial, health, recreational, and social), and its pattern
coefficients recovered these theoretical domains fairly well: items from the same domain
tend to load on the same factor. The factor intercorrelations also included some sizeable
values, which raises the question of whether a general factor underlies the domains. A
Schmid-Leiman transformation lets us examine that.
### Schmid-Leiman Transformation
The Schmid-Leiman (SL) transformation orthogonalises an oblique solution into a hierarchical
one with a general factor and orthogonalised group factors. `efa_schmid_leiman()` performs it:
```{r}
sl_dospert <- efa_schmid_leiman(efa_dospert)
sl_dospert
```
The general-factor column shows that the items load substantially on the general factor
across all domains except the social one, whose items load least strongly on it on average
(mean general-factor loading around .19, against roughly .38 to .55 in the other domains) —
that is, the non-social domains covary considerably, while the social domain stands somewhat
apart. Individual social items vary considerably around that average, so this is a statement
about the domain and not about every item in it.
### McDonald's Omegas
Finally, `efa_reliability()` computes omega estimates and related indices from the SL solution.
You can either specify the indicator-to-factor correspondences via the `factor_map` argument or
let the function determine them automatically, in which case each indicator is assigned to its
highest-loading factor. With a reasonably simple structure here, it is easiest to leave
`factor_map` unspecified — a few items do land on a neighbouring domain's factor, which is
exactly what motivates the explicit map shown afterwards:
```{r}
efa_reliability(sl_dospert)
```
To set the correspondences explicitly instead — for example to match a theoretical assignment —
pass a `factor_map` matrix:
```{r}
efa_reliability(sl_dospert,
factor_map = matrix(c(rep(1, 6), rep(0, 24), # F1: ethical
rep(0, 18), rep(1, 6), rep(0, 6), # F2: recreational
rep(0, 6), 1, 0, 1, 0, 1, rep(0, 19), # F3: financial (1, 3, 5)
rep(0, 12), rep(1, 6), rep(0, 12), # F4: health
rep(0, 24), rep(1, 6), # F5: social
rep(0, 7), 1, 0, 1, 0, 1, rep(0, 18)), # F6: financial (2, 4, 6)
ncol = 6, byrow = FALSE))
```
The columns of `factor_map` are matched to the
group factors **by position**, so the assignment has to be read off the solution rather than
assumed: a column naming the wrong set of items yields a well-formed but meaningless subscale
omega instead of an error. Comparing the two tables, the health row (F4) drops more in omega
subscale than in omega total, because three of its six items actually load higher on the
ethical factor and still contribute their true-score variance to the total through that
factor.
Which coefficient to report depends on the question being asked. Omega hierarchical asks
whether a total score measures a single construct. Omega subscale asks whether a subscale
carries anything beyond the general factor. The H index asks whether a factor is well defined
by its indicators. ECV and PUC together ask whether a unidimensional model would be defensible
at all. Alpha is reported alongside them because it is still widely
expected, but it assumes essentially tau-equivalent items — equal true-score contributions —
whereas factor analysis yields congeneric solutions, for which alpha is only a lower bound. For
a multidimensional scale such as the DOSPERT it is therefore rarely the coefficient to report,
and the omegas answer the question alpha is usually asked to answer. See `?efa_reliability` for
the definition of each coefficient and the references behind it.
## Where to Next
This vignette covered the core workflow. `EFAtools` goes further in several directions, each
of which has its own vignette or article:
* [EFA with ordinal and missing data](Ordinal_and_missing_data.html): analysing ordinal
items with polychoric correlations and handling missing data.
* [Migrating to the efa_* interface](Migrating_to_efa.html): a guide for users coming from
the older uppercase function names.
* [Multigroup EFA](https://mdsteiner.github.io/EFAtools/articles/Multigroup_EFA.html)
(on the package website): comparing factor structures across groups.
* [Simulation and power](https://mdsteiner.github.io/EFAtools/articles/Simulation_and_power.html)
(on the package website): simulating data and running power analyses for EFA.
Run `browseVignettes("EFAtools")` to see the vignettes installed with the package, or visit
the [package website](https://mdsteiner.github.io/EFAtools/) for the full set of articles.