esem: Exploratory Structural Equation Modeling

In bifactory: (Bifactor) ESEM with Continuous (MLR) or Ordinal (WLSMV) Data

Exploratory Structural Equation Modeling

Description

Fits ESEM using lavaan's native efa() block. For continuous data (default), estimation and rotation are integrated in lavaan's ML optimiser—the same single-step approach as Mplus (*1) syntax. For ordered indicators, supply ordered and the call is routed to esem_ordered (set-ESEM / WLSMV; not identical to continuous esem()).

Usage

esem(
  data,
  nfactors,
  indicators = NULL,
  rotation = "geomin",
  target = NULL,
  estimator = "MLR",
  std.lv = TRUE,
  ordered = NULL,
  group = NULL,
  group_equal = NULL,
  missing = "listwise",
  factor_names = NULL,
  rotation_args = list(),
  heywood_fix = TRUE,
  ...
)

Arguments

data	A data.frame containing the observed indicators (and optionally a grouping variable).
nfactors	Integer. Number of latent factors.
indicators	Character vector of item names to include. If NULL (default), all columns of data except group are used.
rotation	Character. Rotation criterion. Passed to lavaan::cfa(rotation = ...). Common choices: "geomin" (default) – oblique geomin, the Mplus default "target" – requires target matrix "varimax" – orthogonal varimax "oblimin" – oblique oblimin "none" – no rotation (confirmatory EFA) The full list of lavaan-supported rotations is in ?lavaan::efaRotate.
target	A numeric matrix (items x factors) for target rotation. Create with make_target. Required when rotation = "target".
estimator	Character. lavaan estimator. Default "MLR" (ML with Huber-White robust SEs and Satorra-Bentler scaled chi-square). Use "WLSMV" for ordered indicators.
std.lv	Logical. Fix factor variances to 1 for identification? Default TRUE (recommended for ESEM).
ordered	Character vector of ordinal item names. When non-NULL, routes to esem_ordered (WLSMV). Default method there is "lavaan"; see ?esem_ordered for the custom "rotation" path.
group	Character. Name of a grouping variable in data for multi-group ESEM (configural by default).
group_equal	Character vector of lavaan parameter labels to constrain equal across groups, e.g. "loadings" (metric) or c("loadings", "intercepts") (scalar).
missing	Character. Missing data handling for continuous fits. Default "listwise"; use "fiml" for full-information ML. When ordered is set, passed to esem_ordered (default "pairwise" there unless you override).
factor_names	Optional character vector of length nfactors. Defaults to F1, F2, ....
rotation_args	Named list of extra arguments passed to the rotation function (e.g. list(geomin.epsilon = 0.001)). When rotation = "geomin" (or any geomin variant) and geomin.epsilon is not supplied, it defaults to 0.0001 / 0.001 / 0.01 for 2 / 3 / 4+ factors respectively, matching Mplus's defaults.
heywood_fix	Logical. Retry rotation with Cholesky unrotation if a standardised loading exceeds 1 (single-group continuous fits only). Default TRUE. If correction runs, std_loadings may differ from lavaan::standardizedsolution(lavaan_fit(x)); see Heywood fix and loadings.
...	Additional arguments passed to lavaan::cfa().

Details

Mplus equivalence

Mplus syntax	bifactory equivalent
⁠F1-F3 BY y1-y15 (*1);⁠	esem(data, nfactors = 3)
⁠(*1)⁠ with geomin (default)	rotation = "geomin"
⁠(*1)⁠ with target	⁠rotation = "target", target = tgt⁠
GROUPING = g;	group = "g"
Metric invariance	group_equal = "loadings"
Scalar invariance	group_equal = c("loadings", "intercepts")

Why not automate cross-loadings from a CFA?

Adding cross-loadings stepwise from modification indices is a different (exploratory CFA) approach and is not recommended: it capitalises on chance, inflates Type I error, and produces a different model on every dataset. ESEM instead estimates all cross-loadings simultaneously, with rotation acting as a mathematical penalty for complexity – giving a reproducible, theory-neutral solution.

Value

An object of class "esem_fit" containing:

lavaan_fit

The lavaan fit object for continuous ESEM. For ordered / custom WLSMV paths see esem_ordered and The lavaan_fit slot on custom WLSMV fits.

syntax

The lavaan model string that was estimated.

nfactors

Number of factors.

rotation

Rotation method used.

factor_names

Factor names.

indicators

Item names included in the model.

call

The matched call.

How this differs from a two-stage EFA + CFA workaround

Many R implementations run EFA first, extract the loading matrix, then paste it into a CFA as starting values. That is an approximation. bifactory instead uses lavaan's efa() block syntax (available since lavaan 0.6-12), which estimates the rotation and the SEM parameters simultaneously – exactly as Mplus does with the (*1) syntax.

Estimator paths

Continuous ESEM / B-ESEM

esem and besem use lavaan::cfa() with a native efa() block (integrated MLR estimation and rotation).

Ordered ESEM

esem(ordered = ...) calls esem_ordered. Default method = "lavaan" (lavaan WLSMV, efa() block, post-hoc rotation). method = "rotation" uses a custom DWLS pipeline (polychoric correlations + GPArotation).

Ordered B-ESEM

besem(ordered = ...) calls besem_ordered. Default method = "rotation" (custom DWLS/WLSMV + orthogonal targetT; this is what run_comparison uses for ordered data). method = "set-esem" fits a lavaan WLSMV bifactor CFA with non-primary specific loadings fixed at zero; it does not match Mplus B-ESEM loadings.

Multi-group invariance

esem_invariance uses lavaan multi-group efa() models plus explicit syntax patches for ordered B-ESEM. See that help page for scope limits.

Missing-data defaults

Defaults differ by entry point; pass missing explicitly when fitting ESEM and B-ESEM separately on the same ordered dataset.

• specify_model / run_comparison: "pairwise" for ordered data, "listwise" for continuous.

• esem_ordered: "pairwise".

• besem_ordered: "listwise" (pipeline still passes missing from the spec when used via run_comparison).

• esem / besem (continuous): "listwise".

Heywood fix and loadings

When heywood_fix = TRUE (default on single-group esem and esem_ordered), std_loadings may reflect a post-hoc rotation correction while lavaan_fit still holds the pre-correction lavaan solution. Prefer std_loadings(x) for reported loadings in that case.

See Also

make_target for target matrices, esem_ordered for ordinal data, esem_compare for ESEM vs CFA comparison, std_loadings for the standardised loading matrix.

Examples

data("HolzingerSwineford1939", package = "lavaan")
d <- HolzingerSwineford1939[, paste0("x", 1:9)]

# Basic ESEM: 3 factors, geomin rotation (Mplus default)
fit <- esem(d, nfactors = 3)
round(std_loadings(fit), 2)
lavaan::fitMeasures(lavaan_fit(fit), c("cfi", "tli", "rmsea", "srmr"))

# Named factors
fit2 <- esem(d, nfactors = 3,
             factor_names = c("Visual", "Textual", "Speed"))


# Target rotation
tgt <- make_target(list(Vis = 1:3, Txt = 4:6, Spd = 7:9), nitems = 9)
fit3 <- esem(d, nfactors = 3, rotation = "target", target = tgt)

# Multi-group configural ESEM
fit_mg <- esem(HolzingerSwineford1939, nfactors = 3,
               indicators = paste0("x", 1:9), group = "sex")

bifactory documentation built on Oct. 28, 2026, 1:06 a.m.