Fits efa_fit() to each of several imputed datasets, aligns the
factor solutions to a common factor space, and pools the resulting estimates
and selected fit quantities across imputations.
efa_mi(
data_list,
p = 0.05,
target_method = c("first_target", "consensus"),
align_unrotated = c("signed_tucker_congruence", "none", "procrustes"),
fit_pool_method = c("D2"),
consensus_args = list(),
procrustes_args = list(),
rmsea_ci_level = 0.9,
rmsr_upper = TRUE,
...
)A list of class c("efa_mi", "EFA_POOLED", "efa", "EFA") containing
pooled estimates, residuals, fit indices, the individual fits, and MI
diagnostics. The trailing legacy classes keep inherits(x, "EFA_POOLED") and
the single-fit EFA accessors and S3 dispatch working. In
addition to the slots inherited from efa_fit() (including SE, CI, and,
on the bootstrap path, replicates), the object carries:
Multiple-imputation diagnostics for each pooled parameter family.
On the bootstrap path: unrot_loadings, h2, residuals, optionally
rot_loadings, Phi, Structure, and fit_indices_descriptive, plus
integer vectors bootstrap_source_failures (replicates the component efa_fit()
could not fit), bootstrap_rotation_failures (replicates whose Procrustes
alignment to the target was invalid), and bootstrap_rotation_valid (those
that entered the pool, B - source - rotation failures). Both paths use the
plain Rubin (1987) df. On the analytic path (se = "information"):
unrot_loadings and uniquenesses, plus, when a rotation was requested,
rot_loadings, h2, and (oblique) Phi and Structure. Each per-family
entry is a list with RIV (relative increase in variance), FMI (the
fraction of missing information, reported as Rubin's asymptotic
\(\lambda = RIV / (1 + RIV)\), equal to lavaan.mi's fmi), and df; the
rotated families on the analytic path additionally carry a method string
recording the gauge alignment used ("gauge_invariant" for communalities and
"signed_permutation_approx" for rotated loadings and, for oblique rotations,
factor correlations and structure coefficients). fit_indices_descriptive, on
the bootstrap path, pools every per-imputation fit index, so the structural
constants among them (df, df_null) appear with a standard error of 0.
On the se = "sandwich" (MI2S) path only: the single efa_fit()
fit on the pooled correlation matrix \(\bar r\) and pooled asymptotic
covariance \(\tilde\Gamma\). Its orig_R is \(\bar r\) and its Gamma
is \(\tilde\Gamma\); the pooled SE, CI, and fit_indices are taken
from it. MI is NULL on this path because the imputation uncertainty is
carried by \(\tilde\Gamma\) rather than by per-parameter Rubin pooling.
A list of length \(m\), where \(m\) is the number of
imputations. Each list element is a data frame or matrix of raw data, or a
correlation matrix. See argument x in efa_fit().
Numeric in \((0, 1)\). One minus the confidence level used for
pooled Wald-type bootstrap/MI confidence intervals when bootstrap replicates
are available. For example, p = .05 gives 95% intervals.
Character. How rotated solutions are aligned across imputations
before pooling: "first_target" (the default) aligns every imputation to the first
imputation's rotated solution, while "consensus" refines a centroid target by
Generalized Procrustes Analysis (orthogonal rotations only). See Aligning solutions
across imputations in Details.
Character. How unrotated loadings are aligned before pooling:
"signed_tucker_congruence" (the default; sign/permutation via Tucker congruence,
anchored on the medoid imputation and returned in the extraction's canonical
gauge), "procrustes" (orthogonal Procrustes to the first imputation), or
"none". See Aligning solutions across imputations in Details.
Character. Currently only "D2" is implemented
for chi-square-type fit. If no chi-square is available, only residual-based
fit and descriptive quantities are returned. See Pooling the model chi-square and
fit indices in Details.
List of additional arguments controlling the
GPA-consensus iteration when target_method = "consensus". Recognised tuning
parameters include the convergence tolerances tol and loss_tol, the
iteration bounds min_iter and max_iter, the target-update damping alpha,
and the multi-start controls multi_start and starts.
List of additional arguments passed to efa_procrustes()
for fixed-target alignment.
Numeric. Confidence level for the RMSEA CI.
Logical. If TRUE, compute RMSR from the unique
off-diagonal residual correlations. If FALSE, use the full off-diagonal
matrix.
Additional arguments passed to efa_fit() (e.g. estimator, rotation, se,
n_factors, N). These select the estimator, rotation, standard-error method, and
fit indices used for every imputation; see efa_fit() for the available options, their
properties, and which combinations are valid.
All conditions are classed (prefix efa_pooled_, or efa_consensus_ for the
consensus target; the dots validation shared with efa_fit() signals
efa_flat_knob_in_dots and efa_renamed_arg) so they can be caught by
class. The ones most likely to be encountered:
Inputs. efa_pooled_min_fits (at least two fits are required);
efa_pooled_mixed_se (every imputation must use the same se).
Alignment. efa_consensus_oblique_unsupported (target_method = "consensus" is orthogonal-only).
Standard errors. efa_pooled_se_unavailable (a warning: pooled SEs could
not be produced, so only point estimates are returned); efa_pooled_no_vcov
and efa_pooled_unreliable_vcov (the analytic "procrustes" path needs a
reliable vcov_unrot_loadings on every fit).
Two-stage (se = "sandwich"). efa_pooled_mi2s_inputs_inconsistent
(every imputation must use se = "sandwich" with the same cor_method);
efa_pooled_mi2s_n_too_small (a warning below 20 imputations);
efa_pooled_mi2s_acov_not_psd (the pooled covariance is indefinite -- use
more imputations); efa_pooled_mi2s_alignment_ignored (a warning that the
alignment arguments do not apply here).
The remaining conditions concern partial or insufficient bootstrap replicates and unequal sample sizes across imputations.
Andreas Soteriades, Markus Steiner
efa_mi() is the multiple-imputation route to handling missing data:
several imputed datasets are each fitted with efa_fit() and the solutions pooled.
A single-fit alternative is full-information maximum likelihood, available
directly in efa_fit() as cor_method = "fiml", which EM-estimates a two-stage
correlation from one raw dataset with missing values. Both feed the same
correlation-scale EFA core and differ only in how the missingness is handled;
FIML is intentionally not routed through efa_mi(), which is a multi-fit
pooler by construction.
The pooling pathway is selected automatically from the se method recorded on
the component efa_fit() fits, which must be identical across imputations:
se = "none": no standard errors are pooled.
se = "information": the per-imputation expected-information standard errors
are pooled with Rubin's (1987) rules (Wald intervals).
se = "sandwich": the two-stage pooled-inputs (MI2S) approach fits a single
model on the Rubin-pooled correlation matrix and asymptotic covariance.
se = "np-boot": the non-parametric bootstrap replicates are re-aligned to
the multiple-imputation target and Rubin-pooled.
On the information and np-boot routes, if pooled standard errors cannot be
produced (for example an unreliable analytic covariance or too few bootstrap
replicates) the pool falls back to point-estimate-only pooling and downgrades
settings$se to "none". The MI2S route is the exception: its single fit
fuses the point estimates and standard errors through the pooled asymptotic
covariance, so a structural failure aborts directly rather than falling back.
The same efa_fit() model is fitted to each imputed dataset and the solutions are
put into a common factor space before averaging. For oblique solutions the
factor intercorrelations are aligned together with the loadings so the model
stays internally consistent.
target_method controls how rotated solutions are aligned. "first_target"
(the default) aligns every imputation to the first imputation's rotated
solution by one Procrustes rotation each. "consensus" instead refines a
centroid target by Generalized Procrustes Analysis (Gower 1975; van Ginkel &
Kroonenberg 2014; Lorenzo-Seva & Van Ginkel 2016). The two give the same pooled
estimate for orthogonal rotations (consensus is just more expensive), and
"consensus" is only supported there. Anchoring on the first imputation can
understate the between-imputation variability when the imputations disagree
substantially, whereas "consensus" is more robust to an atypical first
imputation (van Ginkel & Kroonenberg 2014).
align_unrotated controls how unrotated loadings are aligned before pooling:
"signed_tucker_congruence" (the default) matches them up to factor reordering
and sign changes, "procrustes" aligns them to the first imputation by
orthogonal Procrustes rotation, and "none" averages them as returned by
efa_fit().
The default anchors that matching on the medoid imputation -- the one closest
in aligned squared distance to all the others -- rather than on whichever
imputation happens to come first, so the pooled unrotated solution does not
depend on the order of data_list. The rotated solution is aligned separately,
against a reference chosen by target_method, and still depends on that
reference. The pooled matrix is then returned in the same gauge every
component fit uses, by restoring the constraint that identifies the unrotated
solution: diagonal \(L'L\) for a principal-axis extraction, diagonal
\(L' \Psi^{-1} L\) for maximum likelihood (Anderson & Rubin 1956; Lawley &
Maxwell 1971). Which one
applies is read off the component fits themselves, and a solution in neither
gauge -- an improper one, say -- is left as aligned. Without this step the
average of several aligned solutions sits in a gauge no single fit uses and
cannot be compared element-by-element with an efa_fit() solution. The
rotation is orthogonal and
common to all imputations, so communalities, the total variance accounted for,
the model-implied correlation matrix, the residuals and RMSR are unchanged by
it; only the split of variance across factors moves. "procrustes" and
"none" keep their first-imputation anchor and are returned as aligned.
Point estimates are pooled by arithmetic averaging after alignment. For oblique rotations the structure matrix is recomputed from the pooled pattern matrix and pooled factor correlations, \(Structure = \Lambda \Phi\), and communalities are the diagonal of the reproduced correlation matrix, \(diag(\Lambda \Phi \Lambda')\) for oblique rotations and \(diag(\Lambda \Lambda')\) otherwise. Residuals are not averaged across imputations; they are the pooled observed correlation matrix minus the model-implied correlation of the pooled solution, so RMSR/SRMR are based on these pooled residuals. Both are returned, though the print and summary methods show SRMR only.
The model chi-square and the indices derived from it (RMSEA, ECVI, and the
descriptive AIC/BIC) are pooled with the D2 rule (Li, Meng, Raghunathan &
Rubin, 1991), not arithmetically averaged. Because D2 shrinks the pooled
chi-square in proportion to the between-imputation variability, the pooled RMSEA
can fall below the mean of the per-imputation RMSEAs (as it does in
lavaan.mi); read it together with the per-imputation fit. The incremental
indices CFI (Bentler, 1990) and TLI (Tucker & Lewis, 1973) are instead the
average of the per-imputation indices, which keeps them consistent with the
component fits and avoids the out-of-range values that separately pooling the
model and baseline noncentralities (as lavaan.mi/semTools do) can produce;
those separately pooled noncentralities remain available in mi_diagnostics.
AIC and BIC, if returned, are chi-square-derived descriptive quantities and are
not likelihood-based MI information criteria. They are reported only where the
component fits report them: whenever a component withholds them -- any
cor_method = "fiml" fit, and any fit whose chi-square is a scaled statistic,
such as se = "sandwich" -- the pooled AIC, BIC, and ECVI are NA too, matching
what efa_fit() returns for a single such fit. On the sandwich/MI2S route the
chi-square is the single fit's scaled statistic rather than a D2 pool.
If each component efa_fit() call was run with se = "np-boot", pooled bootstrap SEs
and Wald-type MI confidence intervals are computed for loadings, communalities,
residuals, and, when applicable, factor correlations and structure
coefficients. The unrotated bootstrap replicates are re-aligned to the final MI
target before the within-imputation covariance is estimated, and Rubin pooling
is applied with \(T = Ubar + (1 + 1/m) B\). The confidence level of the pooled
intervals is set by p, not by the component efa_fit() calls' ci.
With se = "information", the analytic unrotated-loading and uniqueness SEs
returned by each fit are pooled element-wise with Rubin's rules
(\(T = Ubar + (1 + 1/m) B\)), with Wald intervals on the plain Rubin (1987)
degrees of freedom (the analytic loadings are asymptotically normal, so the
Barnard-Rubin (1999) adjustment reduces to this form, matching lavaan.mi). NA
propagation is fail-closed: if any imputation is NA at an element, all pooled
outputs for that element are NA. When a rotation was requested, the rotated
loadings, communalities, and (for oblique rotations) factor correlations and
structure coefficients are pooled as well; residual SE pooling is available only
on the bootstrap path. Under align_unrotated = "procrustes" the full unrotated
covariance vcov_unrot_loadings (populated by se = "information") is
propagated through the alignment, so it must be present and reliable on every
fit. The default alignment also mixes loading columns, through the common
canonical-gauge rotation, and so propagates the same covariance; where a fit
does not carry it, the unrotated standard errors are returned as NA rather
than aborting, and the remaining families still pool.
A rotated-loading standard error is conditional on the rotation criterion
(Archer & Jennrich 1973; Jennrich 1973, 1974; Zhang, Preacher, & Jennrich
2012; Zhang & Preacher 2015). For both orthogonal
and oblique rotations the within-imputation variance is therefore each fit's own
criterion-aware delta-method rotated SE (the quantity efa_fit() returns), reused
after a signed-permutation alignment to the MI target, and the
between-imputation variance is the sample variance of the aligned rotated
loadings. This is a deliberate approximation -- each SE is conditional on its
own fit's rotation optimum rather than on a common gauge -- and is flagged by
MI$<param>$method = "signed_permutation_approx". Communalities are
rotation-invariant and pool element-wise. For a fully gauge-consistent rotated
uncertainty, cross-check with se = "np-boot".
With se = "sandwich" (robust SEs from a polychoric/tetrachoric or
continuous-Pearson asymptotic covariance), pooling follows the two-stage,
pooled-inputs approach (Chung & Cai 2019; Sriutaisuk, Liu, Chung, Kim & Gu
2025): the correlation matrix and the asymptotic covariance of its off-diagonal
entries are Rubin-pooled across imputations,
$$\tilde\Gamma = \Gamma_W + \left(1 + \frac{1}{m}\right)\Gamma_B,$$
and a single EFA model is fitted to the pooled correlation with
\(\tilde\Gamma\) as the robust meat (its diagonal as the weights for
estimator = "DWLS"). Because there is only one fit and one rotational gauge, this
route bypasses the per-imputation alignment: target_method and
align_unrotated do not apply. The fitted object carries native scaled-shifted
chi-square statistics and sandwich SEs that already reflect the
multiple-imputation uncertainty, so the chi-square is not D2-pooled and the
likelihood-ratio-based AIC/BIC/ECVI are NA; it is returned in the mi_fit
slot, with the per-imputation fits retained for diagnostics. The pooled fit
uses the same estimate_control() and rotate_control() tuning (including
any rotation-engine extras) as the per-imputation fits. At least 20
imputations are recommended for the scaled-shifted statistic, and more (around
100) at higher rates of missingness (Sriutaisuk et al. 2025). The
polychoric/tetrachoric (ordinal) case is the primary, best-evaluated target; the
continuous-Pearson case uses the same recipe but is less benchmarked.
Anderson, T. W., & Rubin, H. (1956). Statistical inference in factor analysis. In Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability (Vol. 5, pp. 111-150). University of California Press.
Archer, C. O., & Jennrich, R. I. (1973). Standard errors for rotated factor loadings. Psychometrika, 38(4), 581-592.
Barnard, J., & Rubin, D. B. (1999). Small-sample degrees of freedom with multiple imputation. Biometrika, 86(4), 948-955.
Bentler, P. M. (1990). Comparative fit indexes in structural models. Psychological Bulletin, 107(2), 238-246.
Chan, K. W., & Meng, X.-L. (2022). Multiple improvements of multiple imputation likelihood ratio tests. Statistica Sinica, 32, 1489-1514.
Chung, S., & Cai, L. (2019). Alternative multiple imputation inference for categorical structural equation modeling. Multivariate Behavioral Research, 54(3), 323-337.
Gower, J. C. (1975). Generalized Procrustes analysis. Psychometrika, 40(1), 33-51.
Li, K. H., Meng, X.-L., Raghunathan, T. E., & Rubin, D. B. (1991). Significance levels from repeated p-values with multiply-imputed data. Statistica Sinica, 1(1), 65-92.
Jennrich, R. I. (1973). Standard errors for obliquely rotated factor loadings. Psychometrika, 38(4), 593-604.
Jennrich, R. I. (1974). Simplified formulae for standard errors in maximum-likelihood factor analysis. British Journal of Mathematical and Statistical Psychology, 27(1), 122-131.
Lorenzo-Seva, U., & Van Ginkel, J. R. (2016). Multiple imputation of missing values in exploratory factor analysis of multidimensional scales. Anales de Psicologia, 32(2), 596-608.
Rubin, D. B. (1987). Multiple imputation for nonresponse in surveys. Wiley.
Schoenemann, P. H. (1966). A generalized solution of the orthogonal Procrustes problem. Psychometrika, 31(1), 1-10.
Sriutaisuk, S., Liu, Y., Chung, S., Kim, H., & Gu, F. (2025). Evaluating imputation-based fit statistics in structural equation modeling with ordinal data: The MI2S approach. Educational and Psychological Measurement, 85(1), 82-113.
Tucker, L. R., & Lewis, C. (1973). A reliability coefficient for maximum likelihood factor analysis. Psychometrika, 38(1), 1-10.
van Ginkel, J. R., & Kroonenberg, P. M. (2014). Using generalized Procrustes analysis for multiple imputation in principal component analysis. Journal of Classification, 31(2), 242-269.
Zhang, G., & Preacher, K. J. (2015). Factor rotation and standard errors in exploratory factor analysis. Journal of Educational and Behavioral Statistics, 40(6), 579-603.
Zhang, G., Preacher, K. J., & Jennrich, R. I. (2012). The infinitesimal jackknife with exploratory factor analysis. Psychometrika, 77(4), 634-648.
Other factor analysis:
efa_average(),
efa_fit(),
efa_group(),
plot.efa_group(),
print.efa_group()
# create a list of three datasets, mimicking a list you would obtain from
# e.g. mice.
dat_list <- lapply(1:3, function(x) GRiPS_raw[sample(1:nrow(GRiPS_raw), replace = TRUE),])
mod <- efa_mi(dat_list, n_factors = 1, estimator = "ML")
mod
# \donttest{
# add computation of standard errors and CIs
mod <- efa_mi(dat_list, n_factors = 1, estimator = "ML", se = "np-boot")
mod
# }
Run the code above in your browser using DataLab