besem: Bifactor Exploratory Structural Equation Modeling (B-ESEM)

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

Bifactor Exploratory Structural Equation Modeling (B-ESEM)

Description

Fits a Bifactor ESEM model – a general factor (G) loading on all indicators plus domain-specific factors each targeting a subset of indicators. All factors are orthogonal (uncorrelated), matching the Mplus B-ESEM specification on continuous data via lavaan's efa() block. With ordered indicators, the call routes to besem_ordered (default custom DWLS/WLSMV path; see ?besem_ordered).

Usage

besem(
  data,
  specific_factors,
  indicators = NULL,
  g_name = "G",
  estimator = "MLR",
  std.lv = TRUE,
  ordered = NULL,
  group = NULL,
  group_equal = NULL,
  missing = "listwise",
  n_starts = 30L,
  ...
)

Arguments

data	A data.frame of observed indicators.
specific_factors	A named list mapping specific factor names to their primary indicator names. The general factor G is added automatically. Example: list(EX = c("y1","y2"), MD = c("y3","y4"), CI = c("y5","y6")).
indicators	Optional character vector of all indicator names. If NULL (default), all items from specific_factors are used.
g_name	Character. Name for the general factor. Default "G".
estimator	Character. Default "MLR" (ML with Huber-White robust SEs and Satorra-Bentler scaled chi-square).
std.lv	Logical. Fix factor variances to 1. Default TRUE.
ordered	Character vector of ordinal item names. When non-NULL, routes to besem_ordered. For Mplus-aligned ordered B-ESEM use method = "rotation" there (default in run_comparison).
group	Character. Grouping variable for multi-group B-ESEM.
group_equal	Character vector of lavaan equality constraints.
missing	Character. Missing data handling for continuous B-ESEM. Default "listwise". When ordered is set, passed to besem_ordered (default "listwise" there; the pipeline passes "pairwise" from specify_model when used via run_comparison).
n_starts	Integer. Number of random orthogonal starting matrices for the target rotation (forwarded as rstarts to lavaan's rotation engine). Multiple random starts help escape local optima of the rotation criterion, particularly for bifactor models with many specific factors. Default 30L (matches Mplus).
...	Additional arguments passed to lavaan::cfa().

Details

What makes B-ESEM different from ESEM

	ESEM	B-ESEM
Factor structure	k oblique specific factors	1 general + k orthogonal specific
Factor correlations	Freely estimated	All fixed to zero (orthogonal)
Cross-loadings	Estimated via rotation	Estimated via orthogonal target rotation
G factor	None	Loads freely on all items
Rotation	Oblique target/geomin	Orthogonal target

Orthogonality

B-ESEM uses orthogonal target rotation ("targetT" in lavaan), which constrains all factors to be uncorrelated. This means:

• G is uncorrelated with EX, MD, CI

• EX, MD, CI are uncorrelated with each other

This matches Mplus ROTATION = TARGET (orthogonal).

Interpreting results

• G loadings: variance shared across all items regardless of domain

• Specific loadings: domain-specific variance after accounting for G

• omega_h (omega hierarchical): reliability of G (use psych::omega())

• omega_s (omega specific): reliability of each specific factor

Target matrix structure

       G   EX   MD   CI
y1     1    1    0    0   <- G free on all; EX primary; MD/CI targeted to 0
y2     1    1    0    0
y6     1    0    1    0   <- MD item
y9     1    0    0    1   <- CI item

Value

An object of class c("besem_fit", "esem_fit") with the same structure as esem, plus:

g_name

Name of the general factor.

specific_factors

Named list of specific factor assignments.

Mplus target syntax

Mplus continuous / ordered target syntax:

ROTATION = TARGET (orthogonal);
MODEL:
  G BY batEX1-batCI5 (*1);
  EX BY batEX1~1 ... batMD1~0 ... (*1);
  MD BY batEX1~0 ... batMD1~1 ... (*1);
  CI BY batEX1~0 ... batCI1~1 ... (*1);

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".

The lavaan_fit slot on custom WLSMV fits

For besem_ordered(method = "rotation") and esem_ordered(method = "rotation"), $lavaan_fit is often an auxiliary one-factor WLSMV CFA used only to extract DWLS weight matrices—not the fitted ESEM/B-ESEM model. Use std_loadings, parameters, and fitMeasures(x) on the esem_fit wrapper; fitMeasures(x) reads wlsmv_stats when present. Do not interpret summary(x$lavaan_fit), modindices(x), or coef(x) as the rotated solution unless you know the fit used the lavaan efa() path (method = "lavaan" or method = "set-esem" for B-ESEM).

See Also

esem for standard oblique ESEM, besem_ordered for ordinal B-ESEM, make_bifactor_target for the target matrix, generate_mplus_besem_syntax for Mplus comparison.

Examples

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


fit_b <- besem(
  data = d,
  specific_factors = list(
    Visual  = c("x1", "x2", "x3"),
    Textual = c("x4", "x5", "x6"),
    Speed   = c("x7", "x8", "x9")
  ),
  n_starts = 5L
)

summary(fit_b, fit.measures = TRUE, standardized = TRUE)
std_loadings(fit_b)        # rows = items, cols = G + specific factors
factor_correlations(fit_b) # should all be ~0 (orthogonal)

# Compare B-ESEM vs standard ESEM
fit_esem <- esem(d, nfactors = 3)
lavaan::fitMeasures(lavaan_fit(fit_b),    c("cfi", "rmsea", "aic"))
lavaan::fitMeasures(lavaan_fit(fit_esem), c("cfi", "rmsea", "aic"))

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