Package {bifactory}

bifactory: Bifactor ESEM, CFA, and Measurement Invariance with Mplus-Compatible Output

Description

Mplus-compatible CFA, ESEM, and bifactor ESEM (B-ESEM) in R, with multi-group measurement invariance, ESEM-within-CFA conversion, and McDonald's omega reliability. Built on lavaan and GPArotation.

Details

Estimator paths differ by data type and model; see ?esem, ?besem, ?esem_ordered, and ?besem_ordered. Validation summary: system.file("VALIDATION.md", package = "bifactory").

Author(s)

Maintainer: Leon T. De Beer leondb@gmail.com

Authors:

• Leon T. De Beer leondb@gmail.com

See Also

Useful links:

• https://github.com/leondebeer/bifactory

• Report bugs at https://github.com/leondebeer/bifactory/issues

Align ESEM/B-ESEM standardized loadings to a reference

Description

Rotates a fitted solution's standardized loadings by an orthogonal Q so they match a reference matrix (Procrustes) or a deterministic canonical orientation. Useful when comparing bifactor solutions across software, or when you want a reproducible orientation across reruns.

Usage

align_loadings(
  x,
  target = "canonical",
  level = "configural",
  se_method = c("approx", "none")
)

Arguments

Details

Multi-group fits and esem_invariance objects are aligned per group. Item communalities (\Sigma\lambda^2), reliability indices (omega, ECV), and model fit are invariant under orthogonal rotation, so alignment affects only the per-loading partition, not substantive conclusions.

Value

An object of class aligned_loadings: a list with one element per group, each containing

loadings

aligned standardized loading matrix (items x factors)

se

aligned SEs (matching shape) or NULL

Q

the orthogonal rotation matrix applied

residual_max, residual_mean

max / mean abs residual vs target (only when an external reference was supplied)

Examples

data("HolzingerSwineford1939", package = "lavaan")

spec <- specify_model(
  Visual  = c("x1", "x2", "x3"),
  Textual = c("x4", "x5", "x6"),
  Speed   = c("x7", "x8", "x9"),
  data  = HolzingerSwineford1939,
  group = "school"
)


  inv <- esem_invariance(spec)

  # Canonical orientation (no external reference, deterministic)
  aligned <- align_loadings(inv, target = "canonical", level = "configural")
  print(aligned)

  # Within-fit: align all groups to group 1's orientation
  aligned <- align_loadings(inv, target = "group1", level = "configural")


## Not run: 
  # Align to Mplus output (requires a Mplus .out file + MplusAutomation)
  mp  <- MplusAutomation::readModels("besem_inv_configural.out")
  tgt <- extract_mplus_loadings(mp)
  aligned <- align_loadings(inv, target = tgt, level = "configural")

## End(Not run)

Alignment Ratio Check for ICM-CFA Specification

Description

Tests whether a standard CFA (ICM-CFA) is appropriately specified by computing alignment ratios from the manifest correlation matrix. Based on Mehrvarz & Rouder (2026), who prove that in a correctly specified ICM-CFA, alignment ratios must be invariant across all admissible item quadruples, with the common value equal to the squared latent correlation phi^2.

Usage

alignment_check(
  data,
  clusters,
  cfa_fit = NULL,
  is_cor = FALSE,
  min_within_r = 0.05,
  log_sd_thresholds = c(slight = 0.64, moderate = 1.6)
)

Arguments

Details

## Interpreting the log-scale dispersion (sd of log Q) s = \mathrm{sd}(\log Q) is the primary diagnostic. The log scale is natural because alignment ratios are multiplicative objects. Mehrvarz & Rouder (2026, p.24) report empirical terciles of s from their bifactor misalignment simulation:

• s < 0.64 – Slight misalignment. ICM-CFA is approximately consistent with the data; latent correlations are roughly trustworthy.

• s 0.64-1.6 – Moderate misalignment. ICM-CFA is partially misspecified; interpret latent correlations with caution.

• s > 1.6 – High misalignment (top tercile of the paper's simulation, up to ~4.4). ICM-CFA is substantially misspecified and latent correlations are likely inflated. Consider ESEM.

These cutpoints are the 33rd and 66th percentiles of s observed by Mehrvarz & Rouder across a broad range of simulated loading configurations (with \mathrm{sd}(\log \kappa_j) \in [0, 3]), not hard decision boundaries; adjust log_sd_thresholds if your application calls for stricter or looser cutoffs.

## Types of misspecification detected

• Misassignment: items are assigned to the wrong cluster. Mehrvarz & Rouder prove this necessarily inflates phi. The alignment ratios split into three plateaus at phi^2, 1, and 1/phi^2.

• Misalignment: cluster assignments are correct but cross-loadings exist (or the proportionality constraint is violated in a bifactor-like DGP). The ratios are dispersed rather than invariant.

## Number of alignment ratios per pair Mehrvarz & Rouder (2026, p.14) count one alignment ratio per admissible index quadruple: n_Q = \binom{m_A}{2}\binom{m_B}{2}. Because each unordered quadruple \{\ell, \ell'\} \times \{k, k'\} admits two between-cluster matchings ((\ell,k)/(\ell',k') and (\ell,k')/(\ell',k)), both of which equal \varphi^2 under correctly specified ICM-CFA but diverge under misspecification, this function records both matchings per quadruple – yielding 2 \binom{m_A} {2}\binom{m_B}{2} ratios in all_ratios – to maximise the information available for the dispersion diagnostic.

Value

An object of class "alignment_check" (a list) containing:

pair_results

A data.frame with one row per factor pair, giving: geometric mean of alignment ratios, log-scale SD, implied phi, CFA-estimated phi (if supplied), inflation percentage, and verdict.

all_ratios

A named list of raw alignment ratio vectors, one element per factor pair.

recommendation

Character string: overall recommendation.

cor_matrix

The manifest correlation matrix used.

clusters

Resolved cluster assignments (named list).

call

The matched call.

What are alignment ratios?

For any two items \ell, \ell' in cluster A and two items k, k' in cluster B, the alignment ratio is:

Q(\ell, \ell', k, k') = \frac{r_{\ell k} \cdot r_{\ell' k'}}{r_{\ell \ell'} \cdot r_{k k'}}

Under a correctly specified ICM-CFA, all such ratios equal phi^2 and lie in [0, 1] (Mehrvarz & Rouder, 2026, Eq. 5). Dispersion in the ratios – measured by their log-scale standard deviation s = \mathrm{sd}(\log Q) – signals misspecification: either misassignment (items in the wrong cluster) or misalignment (cross-loadings exist but are fixed to zero).

References

Mehrvarz, M., & Rouder, J. N. (2026). The geometry and brittleness of latent correlations in confirmatory factor analysis.

See Also

esem for the recommended follow-up under moderate or high misalignment.

Examples

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

# Vector form of cluster assignment
clusters <- c(
  x1 = "Visual",  x2 = "Visual",  x3 = "Visual",
  x4 = "Textual", x5 = "Textual", x6 = "Textual",
  x7 = "Speed",   x8 = "Speed",   x9 = "Speed"
)


# Run alignment check (data only)
check <- alignment_check(d, clusters)
print(check)

# Run with a fitted CFA to quantify inflation
cfa_model <- "
  Visual  =~ x1 + x2 + x3
  Textual =~ x4 + x5 + x6
  Speed   =~ x7 + x8 + x9
"
cfa_fit <- lavaan::cfa(cfa_model, data = d, std.lv = TRUE)
check   <- alignment_check(d, clusters, cfa_fit = cfa_fit)
print(check)

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

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 ordered-categorical 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"))

Bifactor ESEM for Ordered-Categorical Data

Description

Fits B-ESEM on ordered-categorical indicators. Called from besem when ordered is set. For Mplus-aligned loadings and fit indices use method = "rotation" (default; also used by run_comparison).

Usage

besem_ordered(
  data,
  specific_factors,
  indicators = NULL,
  g_name = "G",
  method = c("rotation", "set-esem"),
  n_starts = 30L,
  r_obs_override = NULL,
  group = NULL,
  group_equal = NULL,
  missing = "listwise",
  std.lv = TRUE,
  ...
)

Arguments

Value

An object of class c("besem_fit_ordered","besem_fit","esem_fit"). For method = "rotation", lavaan_fit is auxiliary (one-factor CFA); use std_rotated_loadings, wlsmv_stats, std_loadings, and fitMeasures(x). For method = "set-esem", lavaan_fit is the fitted bifactor CFA.

Method "set-esem" vs "rotation"

rotation estimates cross-loadings (targeted toward zero), then rotates— same estimand as Mplus B-ESEM WLSMV. set-esem **fixes** non-primary specific loadings at zero in lavaan; G and primary loadings are re-optimized under that harder constraint, so loadings and fit differ from Mplus.

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_ordered, besem, run_comparison

Examples

data("HolzingerSwineford1939", package = "lavaan")

# Derive ordered (5-category Likert) versions of the 9 continuous items
items <- paste0("x", 1:9)
ord <- as.data.frame(lapply(HolzingerSwineford1939[, items], function(v) {
  as.integer(cut(v, breaks = quantile(v, probs = seq(0, 1, 0.2)),
                 include.lowest = TRUE))
}))
names(ord) <- items


fit_b_ord <- besem_ordered(
  data = ord,
  specific_factors = list(
    Visual  = c("x1", "x2", "x3"),
    Textual = c("x4", "x5", "x6"),
    Speed   = c("x7", "x8", "x9")
  ),
  n_starts = 5L
)
summary(fit_b_ord, fit.measures = TRUE, standardized = TRUE)

Open or Copy the bifactory Analysis Template

Description

Opens the shipped analysis template in your editor, or copies it to a destination of your choice. The template walks through the full bifactory pipeline (CFA / ESEM / B-ESEM comparison, reliability, ESEM-within-CFA, factor scores, and optional multi-group invariance) using psych::bfi as the demo dataset.

Usage

bifactory_template(to = NULL, overwrite = FALSE)

Arguments

Value

Invisibly returns the path to the template (or the destination when to is supplied).

Examples

# Find the template path
system.file("templates", "template.R", package = "bifactory")

# Copy the template to a file (here a temporary one)
dest <- file.path(tempdir(), "my_analysis.R")
bifactory_template(to = dest, overwrite = TRUE)

## Not run: 
# Open the template directly in your editor (interactive session only)
bifactory_template()

## End(Not run)

Per-Group Chi-Square Decomposition for Invariance Fits

Description

Extracts the per-group contribution to the WLSMV (or ML) chi-square at each invariance level of an esem_invariance object. Optionally compares against Mplus's Chi-Square Contribution From Each Group table, which is useful for pinpointing whether a basin difference between R and Mplus is driven by a single group or distributed across all of them.

Usage

chisq_decomp(x, mplus_dir = NULL, levels = NULL)

Arguments

Details

lavaan's stat.group is the raw (unscaled) per-group chi-square; Mplus prints the scaled per-group contribution that sums to the scaled total. To compare on the same scale, R per-group values are rescaled by the global ratio chisq.scaled / chisq from fitMeasures.

Value

A data frame of class "chisq_decomp" with columns level, group, chisq_R (rescaled to scaled scale), and when mplus_dir is supplied, mplus_group, chisq_M, delta.

Examples

data("HolzingerSwineford1939", package = "lavaan")

spec <- specify_model(
  Visual  = c("x1", "x2", "x3"),
  Textual = c("x4", "x5", "x6"),
  Speed   = c("x7", "x8", "x9"),
  data  = HolzingerSwineford1939,
  group = "school"
)


inv <- esem_invariance(spec)
chisq_decomp(inv)


## Not run: 
# Supply a folder of Mplus .out files to add side-by-side deltas.
chisq_decomp(inv, mplus_dir = "validation/_bfi_g4_mplus_inv")

## End(Not run)

Extract Parameter Estimates from an esem_fit

Description

Forwards to lavaan on lavaan_fit(x). On custom WLSMV rotation paths, use parameters or std_loadings instead; see The lavaan_fit slot on custom WLSMV fits.

Usage

## S3 method for class 'esem_fit'
coef(object, standardized = FALSE, ...)

Arguments

Value

A data.frame of parameter estimates from lavaan::parameterEstimates(), or the standardized solution from lavaan::standardizedsolution() when standardized = TRUE.

Compare Pipeline Results with an EWC Model

Description

Appends fit indices from a fit_ewc result to the pipeline comparison table, producing a unified data frame with CFA, ESEM, B-ESEM, and EWC columns side by side.

Usage

compare_ewc(results, ewc)

Arguments

Value

An object of class c("ewc_comparison","data.frame"). Print with print() for a formatted table.

See Also

fit_ewc

Compare R and Mplus Standardised Loadings

Description

Builds a long-format data frame with STDYX loadings, standard errors, z-scores, p-values, and R - Mplus differences for every loading (primary and cross) across CFA, ESEM, and B-ESEM. Mplus columns are NA when no Mplus results are present in the pipeline object.

Usage

compare_loadings(results)

Arguments

Value

A data frame with columns:

model

CFA, ESEM, or BESEM.

factor

Factor name (lowercase).

item

Item name (lowercase).

loading_type

"primary" or "cross".

R_std, R_se, R_z, R_p

R estimates.

Mplus_std, Mplus_se, Mplus_z, Mplus_p

Mplus estimates (NA if unavailable).

diff_std

R - Mplus standardised loading difference (NA if unavailable).

Examples

data("HolzingerSwineford1939", package = "lavaan")

spec <- specify_model(
  Visual  = c("x1", "x2", "x3"),
  Textual = c("x4", "x5", "x6"),
  Speed   = c("x7", "x8", "x9"),
  data    = HolzingerSwineford1939,
  label   = "Holzinger-Swineford"
)


results <- run_comparison(spec, n_starts = 5L)
lc <- compare_loadings(results)
head(lc)

# Primary loadings only (Mplus_* columns are NA without a Mplus run)
lc[lc$loading_type == "primary", ]

Compute Reliability Indices for CFA, ESEM, and B-ESEM

Description

Computes McDonald's omega reliability indices for all three models in a pipeline result, following the bifactor reporting framework recommended by Morin, Arens & Marsh (2016) and Rodriguez, Reise & Haviland (2016).

Usage

compute_indices(results)

Arguments

Details

For the **total composite** (all items):

• omega_total – total reliability (rotation-invariant for orthogonal models).

• omega_H – hierarchical omega: G factor's contribution to total-score reliability (B-ESEM only).

• ECV – explained common variance: G's share of all common variance (B-ESEM only).

• H(G) – construct replicability of G (B-ESEM only).

For each **subscale** (items of specific factor F_s):

• omega_S – specific factor's unique contribution to subscale reliability: (\sum_{i \in s} \lambda_{s,i})^2 / [(\sum_{i \in s} \lambda_{s,i})^2 + \sum_{i \in s} \psi_i]. Uses target loadings and item residuals only; G does not enter the denominator because \psi_i already has G partialled out.

• omega_sub – total reliability of the subscale sum score (G + specific factor combined).

• omega_H_sub – G's contribution to subscale reliability.

• ECV_s – G's share of common variance within the subscale.

• H(Fs) – construct replicability of the specific factor.

Value

An object of class "reliability_indices" – a named list with elements cfa, esem, besem, each containing the computed indices for that model, plus alpha (Cronbach's alpha per subscale and for G). Pass to print() for a formatted table.

References

McDonald, R. P. (1999). Test theory: A unified treatment. Erlbaum.

Rodriguez, A., Reise, S. P., & Haviland, M. G. (2016). Evaluating bifactor models: Calculating and interpreting statistical indices. Psychological Methods, 21(2), 137-150.

Morin, A. J. S., Arens, A. K., & Marsh, H. W. (2016). A bifactor exploratory structural equation modeling framework for the identification of distinct sources of construct-relevant psychometric multidimensionality. Structural Equation Modeling, 23(1), 116-139.

Examples

data("HolzingerSwineford1939", package = "lavaan")

spec <- specify_model(
  Visual  = c("x1", "x2", "x3"),
  Textual = c("x4", "x5", "x6"),
  Speed   = c("x7", "x8", "x9"),
  data    = HolzingerSwineford1939,
  label   = "Holzinger-Swineford"
)


results <- run_comparison(spec, n_starts = 5L)
indices <- compute_indices(results)
print(indices)

Deprecated: use compute_indices instead

Description

Deprecated: use compute_indices instead

Usage

compute_omega(results)

Arguments

Value

A reliability_indices object.

Shared documentation: estimator paths and defaults

Description

Shared documentation: estimator paths and defaults

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

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.

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

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 ordered-categorical 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")

Compare ESEM Against a Standard CFA

Description

Fits a user-specified CFA model on the same data and compares it against an esem_fit object using fit indices and, where applicable, a chi-square difference test (when models are nested).

Usage

esem_compare(esem_model, cfa_model, data = NULL, estimator = NULL, ...)

Arguments

Value

A list of class "esem_comparison" with:

fit_table

A data.frame of fit indices for both models.

esem_fit

The original esem_fit object.

cfa_fit

The fitted lavaan CFA object.

lavtest

Output of lavaan::lavTestLRT() if models are nested, else NULL.

See Also

esem

Examples

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


# Fit ESEM
esem_result <- esem(d, nfactors = 3)

# Define comparison CFA (no cross-loadings)
cfa_model <- "
  Visual  =~ x1 + x2 + x3
  Textual =~ x4 + x5 + x6
  Speed   =~ x7 + x8 + x9
"
comparison <- esem_compare(esem_result, cfa_model, data = d)
print(comparison)

Measurement Invariance Testing for ESEM Models

Description

Tests configural, weak (metric), strong (scalar), and strict invariance for an ESEM model across groups. Returns a formatted table of fit indices and chi-square difference tests, analogous to Mplus's multi-group output.

Usage

esem_invariance(
  spec,
  model = c("esem", "besem"),
  missing = NULL,
  verbose = TRUE,
  ...
)

Arguments

Value

An object of class "esem_invariance" containing:

table

Data frame with fit indices and D-statistics. print() renders it formatted.

models

Named list of fit objects: configural, weak, strong, strict. NULL entries indicate a model that failed to fit.

lrt

Named list of lavTestLRT() outputs: weak, strong, strict.

spec

The original model specification.

model

Character: "esem" or "besem".

fallback_from, fallback_note

If B-ESEM configural fails to converge, the result may be from an **ESEM** fallback (no general factor); not comparable df-for-df to Mplus B-ESEM.

Scope and caveats

• **Ordered B-ESEM:** Validated against Mplus on BFI for G \in \{2, 3, 4\} (max |\Delta\mathrm{CFI}| \le 0.001, |\Delta\mathrm{SRMR}| \le 0.002, df match). At G \ge 5 with \ge 18 items a warning is issued; verify externally before reporting.

• **WLSMV SRMR:** The invariance table applies an Mplus-style SRMR denominator correction for ordered models where implemented.

• **B-ESEM configural failure:** After retries, may fall back to ESEM (fallback_from = "besem"). The printed table is structurally different from B-ESEM; do not compare to Mplus B-ESEM output.

• **Empty categories:** Ordered multi-group fits may recode categories absent in one group (see verbose output).

See system.file("VALIDATION.md", package = "bifactory") for a summary.

Constraint mapping

| Level | Continuous (MLR) | Ordered (WLSMV/Theta) | |—|—|—| | Configural | free | free | | Weak | loadings | loadings | | Strong | loadings + intercepts | loadings + thresholds | | Strict | + residuals | + residuals |

Scaled chi-square difference tests

Simple subtraction of scaled chi-square values is not valid for MLR or WLSMV. lavaan::lavTestLRT() is used, which applies the Satorra-Bentler (2001) correction for MLR and a mean-variance-adjusted difference test for WLSMV. WLSMV difference test results may differ numerically from Mplus's DIFFTEST procedure.

Identification across groups

At the configural level std.lv = TRUE fixes factor variances to 1 in all groups. Under weak invariance lavaan automatically frees factor variances in Group 2+ and keeps them at 1 in Group 1. Under strong invariance factor means are freed in Group 2+ and fixed to 0 in Group 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.

See Also

specify_model, esem, esem_ordered

Examples

data("HolzingerSwineford1939", package = "lavaan")

spec <- specify_model(
  Visual  = c("x1", "x2", "x3"),
  Textual = c("x4", "x5", "x6"),
  Speed   = c("x7", "x8", "x9"),
  data  = HolzingerSwineford1939,
  group = "school",
  label = "Holzinger-Swineford"
)


# Continuous ESEM measurement invariance across schools
inv <- esem_invariance(spec)
print(inv)

# Access individual model fits
summary(inv$models$strong, fit.measures = TRUE, standardized = TRUE)
lavaan::lavTestScore(inv$models$strong$lavaan_fit)

# Ordered data (WLSMV) and B-ESEM invariance follow the same pattern:
# add `ordered = TRUE` to specify_model(), then optionally `model = "besem"`.
items <- paste0("x", 1:9)
ord <- as.data.frame(lapply(HolzingerSwineford1939[, items], function(v) {
  as.integer(cut(v, breaks = quantile(v, probs = seq(0, 1, 0.2)),
                 include.lowest = TRUE))
}))
names(ord) <- items
ord$school <- HolzingerSwineford1939$school

spec_ord <- specify_model(
  Visual = c("x1", "x2", "x3"), Textual = c("x4", "x5", "x6"),
  Speed = c("x7", "x8", "x9"),
  data = ord, group = "school", ordered = TRUE
)
inv_besem <- esem_invariance(spec_ord, model = "besem")
print(inv_besem)

ESEM for Ordered-Categorical (Likert) Data

Description

Fits ESEM on ordered-categorical indicators (WLSMV / Theta parameterization). Called automatically from esem when ordered is set.

Usage

esem_ordered(
  data,
  nfactors,
  indicators = NULL,
  rotation = "target",
  target = NULL,
  factor_names = NULL,
  method = c("lavaan", "rotation"),
  n_starts = 100L,
  r_obs_override = NULL,
  std.lv = TRUE,
  missing = "pairwise",
  group = NULL,
  group_equal = NULL,
  heywood_fix = TRUE,
  n_obs = NULL,
  ...
)

Arguments

Details

**Default** (method = "lavaan"): lavaan's efa() block with WLSMV (unrestricted model, post-hoc rotation, delta-method SEs)—aligned with Mplus ESTIMATOR = WLSMV; ROTATION = TARGET (oblique).

**Alternate** (method = "rotation"): custom DWLS pipeline (polychoric correlations via psych, GPArotation, sandwich + rotation Jacobian SEs).

## Mplus equivalence

Mplus:
  ANALYSIS:
    ESTIMATOR = WLSMV;
    ROTATION  = TARGET;
    PARAMETERIZATION = THETA;
  MODEL:
    EX MD CI BY item1-item18 (*1);

R (this function):
  esem_ordered(data, nfactors = 3, rotation = "target", target = tgt,
               factor_names = c("EX","MD","CI"))

## Algorithm (method = "lavaan") lavaan's efa() block fits an unrestricted k-factor model under WLSMV, applies post-hoc rotation, and propagates standard errors through the rotation transformation (Asparouhov & Muthen, 2009).

Value

An object of class c("esem_fit_ordered", "esem_fit"). For method = "lavaan", lavaan_fit is the fitted ESEM model. For method = "rotation", lavaan_fit is an auxiliary one-factor CFA; use std_rotated_loadings, wlsmv_stats, and std_loadings. Slot estimator is "WLSMV" or "DWLS" depending on path.

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

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.

References

Asparouhov, T., & Muthen, B. (2009). Exploratory structural equation modeling. Structural Equation Modeling, 16(3), 397–438.

See Also

esem for continuous data, besem_ordered for bifactor ordered ESEM.

Examples

data("HolzingerSwineford1939", package = "lavaan")

# Derive ordered (5-category Likert) versions of the 9 continuous items
items <- paste0("x", 1:9)
ord <- as.data.frame(lapply(HolzingerSwineford1939[, items], function(v) {
  as.integer(cut(v, breaks = quantile(v, probs = seq(0, 1, 0.2)),
                 include.lowest = TRUE))
}))
names(ord) <- items

tgt <- make_target(
  list(Visual = c("x1", "x2", "x3"),
       Textual = c("x4", "x5", "x6"),
       Speed = c("x7", "x8", "x9")),
  item_names = items
)


fit_ord <- esem_ordered(
  data         = ord,
  nfactors     = 3,
  indicators   = items,
  rotation     = "target",
  target       = tgt,
  factor_names = c("Visual", "Textual", "Speed"),
  n_starts     = 5L
)

summary(fit_ord, fit.measures = TRUE, standardized = TRUE)
std_loadings(fit_ord)
factor_correlations(fit_ord)

Generate ESEM-within-CFA Lavaan Syntax

Description

Converts a fitted ESEM solution into explicit lavaan CFA syntax following the ESEM-within-CFA approach (Marsh et al. 2014). The syntax uses unstandardised ESEM loadings as starting values (start(v)*item) or fixed values (v*item) depending on the referent scheme.

Usage

ewc_syntax(esem_fit, spec, referents = NULL, var_fixed = TRUE)

Arguments

Value

A character string of lavaan model syntax. Use cat() to inspect, or pass directly to fit_ewc via custom_syntax.

References

Marsh, H. W., Morin, A. J. S., Parker, P. D., & Kaur, G. (2014). Exploratory structural equation modeling. Annual Review of Clinical Psychology, 10, 85-110.

See Also

find_ewc_referents, fit_ewc

Extract standardized loadings from MplusAutomation output

Description

Convenience helper for use as the target of align_loadings(): pulls the STDYX-standardized loadings from a MplusAutomation::readModels() result and returns one matrix per group, in items-by-factors form with names matching common conventions.

Usage

extract_mplus_loadings(
  mp_out,
  standardized = "stdyx.standardized",
  items = NULL,
  factors = NULL
)

Arguments

Value

A named list of matrices, one per group, each with rownames (items) and colnames (factors).

Extract Factor Correlations from an esem_fit

Description

Extracts the estimated correlation matrix among latent factors.

Usage

factor_correlations(x, digits = 3)

Arguments

Value

A symmetric matrix of factor correlations.

Extract Latent Factor Scores

Description

Returns a data frame of estimated latent factor scores, one column per factor and one row per observation. Works with any esem_fit, plain lavaan object, or esem_invariance result.

Usage

factor_scores(
  x,
  method = c("regression", "bartlett"),
  level = "auto",
  align = NULL,
  dCFI_cutoff = -0.01,
  ...
)

Arguments

Value

A data.frame with one column per latent factor (named by the model's factor labels) and one row per observation. A group column is prepended automatically for multi-group models. Rows for missing observations contain NA in all factor columns.

Examples

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


# Single fit
fit    <- esem(d, nfactors = 3)
scores <- factor_scores(fit)
head(scores)

# Merge back to original data
d_with_scores <- cbind(d, factor_scores(fit))

# Bartlett method
scores_b <- factor_scores(fit, method = "bartlett")

Auto-Select EWC Referent Items

Description

For each factor, selects the primary indicator with the largest absolute unstandardised loading from the ESEM solution. This is the item that will anchor that factor's identification in the EWC model.

Usage

find_ewc_referents(esem_fit, spec)

Arguments

Details

Referents are chosen solely by maximum |loading| on the primary factor. You can override automatic selection by supplying a custom named vector to the referents argument of ewc_syntax or fit_ewc.

Value

A named character vector: names are factor names, values are the selected referent item names (in the original case from spec).

See Also

ewc_syntax, fit_ewc

Extract Fit Measures from an esem_fit

Description

Forwards to lavaan::fitMeasures() when the fit used lavaan's efa() path. When wlsmv_stats is present (custom DWLS/WLSMV from besem_ordered(method = "rotation") or esem_ordered(method = "rotation")), returns those indices instead.

Usage

## S3 method for class 'esem_fit'
fitMeasures(
  object,
  fit.measures = c("cfi", "tli", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper", "srmr",
    "aic", "bic"),
  ...
)

Arguments

Value

A named numeric vector of fit indices (by default CFI, TLI, RMSEA with its confidence interval, SRMR, and, on lavaan paths, AIC and BIC).

Fit an ESEM-within-CFA Model

Description

Generates EWC lavaan syntax from a fitted ESEM solution and estimates it as a standard lavaan::cfa() model – no rotation required.

Usage

fit_ewc(
  esem_fit,
  spec,
  referents = NULL,
  var_fixed = TRUE,
  missing = NULL,
  custom_syntax = NULL,
  ...
)

Arguments

Value

An object of class "ewc_fit":

lavaan_fit

lavaan fit object; all lavaan generics work on it.

syntax

lavaan model string used.

referents

Named vector of referent items.

var_fixed

Identification mode used.

estimator

Estimator ("MLR" for continuous, "DWLS" for ordered).

spec

The spec object.

See Also

ewc_syntax, find_ewc_referents, compare_ewc

Mplus-Matched Fit Indices for a Single Model

Description

Returns CFI, TLI, RMSEA (with 90% CI) and SRMR formatted as a named character vector. Uses the .scaled variants when the fit was estimated with MLR/WLSMV, matching Mplus's default output. For B-ESEM / ESEM WLSMV fits that carry pre-computed Mplus-matched values in $wlsmv_stats, those values are used directly and the RMSEA CI is computed from the scaled chi-square via non-central chi-squared inversion (MacCallum, Browne & Sugawara, 1996).

Usage

fit_indices(fit)

Arguments

Value

A named character vector with elements "CFI", "TLI", "RMSEA [90% CIs]", and "SRMR". Stack rows with rbind() (wrapped in noquote()) to build a comparison table.

Examples

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


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

noquote(rbind(
  ESEM  = fit_indices(fit_e),
  BESEM = fit_indices(fit_b)
))

Generate Mplus B-ESEM Syntax

Description

Generates the correct Mplus syntax for Bifactor ESEM with orthogonal target rotation, using separate BY statements per factor with (*1).

Usage

generate_mplus_besem_syntax(
  specific_factors,
  all_indicators,
  g_name = "G",
  cfa_factors = NULL,
  regressions = NULL,
  covariances = NULL,
  data_file = "mydata.dat",
  missing_code = 999,
  output_path = NULL
)

Arguments

Value

Mplus syntax string (invisibly). Writes file if output_path supplied.

Generate Mplus Syntax for ESEM with Target Rotation

Description

Generates a complete Mplus .inp file for ESEM with target rotation, using the correct Mplus syntax: separate BY statements per factor, each with its own (*1) label. Primary items always appear first, followed by cross-loading items targeted to zero.

Usage

generate_mplus_syntax(
  factors,
  cfa_factors = NULL,
  regressions = NULL,
  covariances = NULL,
  data_file = "mydata.dat",
  missing_code = 999,
  estimator = "MLR",
  output_path = NULL
)

Arguments

Details

Correct output format:

EX BY
  batEX1 batEX2 ... batEX8       <- primary items first
  batMD1~0 ... batCI5~0 (*1);    <- cross-loadings after

MD BY
  batMD1 batMD2 ... batMD5       <- primary items first
  batEX1~0 ... batCI5~0 (*1);    <- cross-loadings after

Value

The Mplus syntax as a character string (invisibly). Writes file if output_path is supplied.

Examples

# Pure syntax generation: builds a Mplus .inp string, no Mplus install needed.
syntax <- generate_mplus_syntax(
  factors = list(
    Visual  = c("x1", "x2", "x3"),
    Textual = c("x4", "x5", "x6"),
    Speed   = c("x7", "x8", "x9")
  ),
  output_path = file.path(tempdir(), "esem_model.inp")
)
cat(syntax)

Extract the Generated lavaan Syntax

Description

Extract the Generated lavaan Syntax

Usage

get_syntax(x, cat = TRUE)

Arguments

Value

The lavaan model string (invisibly).

Extract the Underlying lavaan Fit Object

Description

Returns the raw lavaan object stored in x$lavaan_fit. For continuous ESEM/B-ESEM and method = "lavaan" ordered fits, this is the fitted model. For besem_ordered(method = "rotation") (and esem_ordered(method = "rotation")), it is often an **auxiliary** one-factor CFA used for DWLS weights only—see The lavaan_fit slot on custom WLSMV fits. Use lavaan::lavTestScore(), lavaan::lavInspect(), etc. only when that slot represents the model you intend to analyse.

Usage

lavaan_fit(x)

Arguments

Value

A lavaan object.

Examples

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


fit <- esem(d, nfactors = 3)
lav <- lavaan_fit(fit)

# Inspect the underlying lavaan model
lavaan::lavInspect(lav, "cor.lv")   # factor correlations
lavaan::fitMeasures(lav, c("cfi", "rmsea"))

Create a Bifactor Target Matrix

Description

Builds the target loading matrix for B-ESEM. The general factor G has target 1 for all items (free to load everywhere). Each specific factor has target 1 for its primary items and 0 for all others.

Usage

make_bifactor_target(specific_factors, indicators, g_name = "G")

Arguments

Value

A numeric matrix (items x factors) of class "esem_target". Column order: G first, then specific factors in list order.

Examples

tgt <- make_bifactor_target(
  specific_factors = list(EX = c("y1","y2","y3"),
                          MD = c("y4","y5","y6")),
  indicators       = paste0("y", 1:6)
)
print(tgt)

Create a Target Loading Matrix for Target Rotation

Description

Constructs a target matrix suitable for use with rotation = "target" in esem. Items listed in keys are assigned a target value of 1 for their primary factor; all other cells are set to 0 (penalised to be near zero) or NA (free, no penalty).

Usage

make_target(keys, nitems = NULL, item_names = NULL, cross_loading_value = 0)

Arguments

Details

## Target vs. Soft Target Rotation

- **Hard target** (cross_loading_value = 0): Cross-loadings are penalised toward zero. Use this when you have strong theory. - **Soft target** (cross_loading_value = NA): Only primary loadings are targeted; cross-loadings are completely free. Use this when you are less certain about the zero pattern.

## Items Loading on Multiple Factors

An item can appear in multiple keys entries if it is expected to have meaningful loadings on more than one factor. In that case its target value will be 1 for both factors and cross_loading_value elsewhere.

Value

A numeric matrix (items x factors) with row names set to item names and column names set to factor names from keys.

See Also

esem

Examples

# Simple 15-item, 3-factor target (integer keys)
tgt <- make_target(
  keys   = list(Extrav = 1:5, Agree = 6:10, Open = 11:15),
  nitems = 15
)

# Named items
tgt2 <- make_target(
  keys       = list(F1 = c("y1", "y2", "y3"), F2 = c("y4", "y5", "y6")),
  item_names = paste0("y", 1:6)
)

# Soft target (cross-loadings free)
tgt_soft <- make_target(
  keys                 = list(F1 = 1:5, F2 = 6:10),
  nitems               = 10,
  cross_loading_value  = NA
)

Extract Modification Indices from an esem_fit

Description

Forwards to lavaan::modindices() on lavaan_fit(x). For custom WLSMV paths where lavaan_fit is an auxiliary one-factor CFA, indices refer to that auxiliary model, not the rotated ESEM/B-ESEM solution; see The lavaan_fit slot on custom WLSMV fits.

Usage

## S3 method for class 'esem_fit'
modindices(object, sort. = TRUE, maximum.number = 20, ...)

Arguments

Value

A data.frame of modification indices, as returned by lavaan::modindices().

Display Model Parameters

Description

Prints a formatted parameter table – standardized loadings, standard errors, z-values, and p-values – for one fit object or all three models in a pipeline. Analogous to the STDYX section of Mplus output.

Usage

parameters(
  x,
  model = "all",
  type = "loadings",
  suppress = 0,
  digits = 3,
  highlight_primary = TRUE
)

Arguments

Value

A data frame of class bifactory_parameters holding the standardized parameter table (plus a model column when x is a pipeline). Printing the object renders the formatted, colour-coded table; assign the result to use the values without console output.

Parse Mplus Polychoric Correlation Matrix from Output File

Description

Reads a Mplus .out file and assembles the polychoric (tetrachoric) correlation matrix that Mplus prints under SAMPLE STATISTICS when SAMPSTAT is requested in the OUTPUT section.

Usage

parse_mplus_polychoric(out_file)

Arguments

Details

The matrix can be passed directly to besem_ordered via the r_obs_override argument so that R's DWLS optimisation uses the same polychoric estimates as Mplus, eliminating rotational discrepancies caused by differences in the polychoric estimators.

Value

A named, symmetric numeric matrix of polychoric correlations with 1s on the diagonal. Row and column names are the item names as reported by Mplus (typically upper-case).

Examples

## Not run: 
# Reads a polychoric matrix from a Mplus .out file (SAMPSTAT output),
# so it needs a .out produced by a licensed Mplus run.
R_mplus <- parse_mplus_polychoric(
  "path/to/besem_measurement.out"
)
dim(R_mplus)   # 18 x 18

# Use as input to besem_ordered so R rotates from the same matrix as Mplus
fit_b <- besem_ordered(
  data             = mydata,
  specific_factors = factor_items,
  r_obs_override   = R_mplus
)

## End(Not run)

Partial Measurement Invariance Testing

Description

Given a failed invariance level, uses a greedy score-test loop (lavTestScore + group.partial) to identify and free the minimum set of non-invariant equality constraints needed to restore acceptable fit (DCFI >= -0.010). Downstream invariance levels are then re-tested under the same partial constraints.

Usage

partial_invariance(
  inv,
  level,
  max_free = 10L,
  delta_cfi_cutoff = -0.01,
  verbose = TRUE
)

Arguments

Value

An object of class "esem_partial_invariance" with:

$freed_params

Data frame: round, label, group_name, score (LM), delta_cfi (baseline-relative), converged.

$partial_fit

esem_fit at level with freed params.

$downstream

Named list of esem_fit for levels above target.

$downstream_lrt

lavTestLRT() comparisons among downstream.

$table

Combined fit table (original + partial + downstream).

$converged

Logical: did DCFI pass cutoff?

$group_partial

Character vector: final group.partial labels.

Limitations

B-ESEM models are not supported. Use lavaan::lavTestScore(inv$models[[level]]$lavaan_fit) directly.

See Also

esem_invariance

Plot Sorted Alignment Ratios

Description

Produces a sorted alignment ratio plot (as in Figure 2 of Mehrvarz & Rouder 2026) for one or all factor pairs. Flat horizontal spread = invariant (ICM-CFA tenable). Dispersion = misspecification.

Usage

## S3 method for class 'alignment_check'
plot(x, pair = NULL, ...)

Arguments

Value

Called for its side effect of drawing the alignment-ratio plot; invisibly returns NULL.

Print method for aligned_loadings

Description

Print method for aligned_loadings

Usage

## S3 method for class 'aligned_loadings'
print(x, digits = 3, ...)

Arguments

Value

Invisibly returns x; called for the side effect of printing the aligned loading matrices.

Print Method for alignment_check

Description

Print Method for alignment_check

Usage

## S3 method for class 'alignment_check'
print(x, ...)

Arguments

Value

Invisibly returns x; called for the side effect of printing the alignment-ratio summary.

Print Method for besem_fit

Description

Print Method for besem_fit

Usage

## S3 method for class 'besem_fit'
print(x, ...)

Arguments

Value

Invisibly returns x; called for the side effect of printing a formatted B-ESEM model overview to the console.

Print a Parameter Table

Description

Renders the formatted, colour-coded parameter table produced by parameters. Called automatically when the result of parameters() is not assigned.

Usage

## S3 method for class 'bifactory_parameters'
print(x, ...)

Arguments

Value

x, invisibly.

Print Method for esem_comparison

Description

Print Method for esem_comparison

Usage

## S3 method for class 'esem_comparison'
print(x, ...)

Arguments

Value

Invisibly returns x; called for the side effect of printing the comparison table.

Print Method for esem_comparison_pipeline

Description

Print Method for esem_comparison_pipeline

Usage

## S3 method for class 'esem_comparison_pipeline'
print(x, hints = TRUE, ...)

Arguments

Value

Invisibly returns x; called for the side effect of printing the comparison-pipeline overview.

Print Method for esem_fit

Description

Print Method for esem_fit

Usage

## S3 method for class 'esem_fit'
print(x, ...)

Arguments

Value

Invisibly returns x; called for the side effect of printing a formatted model overview to the console.

Print Method for esem_invariance

Description

Print Method for esem_invariance

Usage

## S3 method for class 'esem_invariance'
print(x, ...)

Arguments

Value

Invisibly returns x; called for the side effect of printing the invariance comparison table.

Print Method for esem_partial_invariance

Description

Print Method for esem_partial_invariance

Usage

## S3 method for class 'esem_partial_invariance'
print(x, ...)

Arguments

Value

Invisibly returns x; called for the side effect of printing the partial-invariance summary.

Print an esem_spec Object

Description

Print an esem_spec Object

Usage

## S3 method for class 'esem_spec'
print(x, ...)

Arguments

Value

Invisibly returns x; called for the side effect of printing the model specification.

Print a Target Rotation Matrix

Description

Displays a make_target() result in a compact, readable format, marking primary loadings (1), penalised cells (0), and free cells (NA).

Usage

## S3 method for class 'esem_target'
print(x, ...)

Arguments

Value

Invisibly returns x; called for the side effect of printing the target loading matrix.

Print Method for ewc_comparison

Description

Displays the fit-index comparison table returned by compare_ewc.

Usage

## S3 method for class 'ewc_comparison'
print(x, ...)

Arguments

Value

Invisibly returns x.

Print Method for ewc_fit

Description

Compact fit summary (CFI/TLI/RMSEA/SRMR and \chi^2) for an ESEM-within-CFA fit.

Usage

## S3 method for class 'ewc_fit'
print(x, ...)

Arguments

Value

Invisibly returns x.

Print a reliability_indices Object

Description

Displays McDonald's omega and Cronbach's alpha indices following the Morin / Rodriguez et al. reporting framework.

Usage

## S3 method for class 'reliability_indices'
print(x, ...)

## S3 method for class 'omega_result'
print(x, ...)

Arguments

Value

Invisibly returns x; called for the side effect of printing the reliability table.

Refine B-ESEM Rotation Using Mplus Solution as Warm Start

Description

The bifactor target rotation criterion surface has many local minima. R and Mplus sometimes converge to different ones, leading to different partitioning of variance between G and specific factors (while fit indices remain identical because they are rotation-invariant).

Usage

refine_rotation(results)

Arguments

Details

This function uses orthogonal Procrustes rotation to compute the rotation matrix \mathbf{T} that maps R's unrotated loading matrix toward the Mplus STDYX solution, then evaluates whether the target criterion at that \mathbf{T} is lower than R's best random-start criterion. If so, the better rotation is adopted. This is methodologically sound because:

1. The same published criterion function is still being minimised.

2. Procrustes only provides a well-informed starting point; the final solution is the converged GPArotation::targetT optimum from that start.

3. The approach is analogous to Mansolf and Reise's (2016) recommendation to use Schmid-Leiman solutions as warm starts.

Value

The same results object with fit_besem updated if a better rotation was found, otherwise unchanged. A message reports whether the criterion improved. The comparison_table is also updated.

References

Mansolf, M., and Reise, S. P. (2016). Exploratory bifactor analysis: The Schmid-Leiman orthogonalization and Jennrich-Bentler analytic rotations. Multivariate Behavioral Research, 51(5), 698–717.

Examples

## Not run: 
# Requires run_comparison() results that include a Mplus rotation reference
# (results$mplus_results non-NULL), so a licensed Mplus install is needed.
results <- run_comparison(spec, mplus_folder = tempfile("mplus_"))
results <- refine_rotation(results)
omega   <- compute_omega(results)
print(omega)

## End(Not run)

Run the Full CFA / ESEM / B-ESEM Comparison Pipeline

Description

Takes a model specification from specify_model and automatically fits CFA, ESEM, and B-ESEM in R, optionally runs all three in Mplus, and returns a comparison table plus all fitted objects.

Usage

run_comparison(
  spec,
  mplus_folder = NULL,
  mplus_command = "Mplus",
  run_alignment = TRUE,
  group_equal = NULL,
  n_starts = 30L
)

Arguments

Value

An object of class "esem_comparison_pipeline" containing:

spec

The original model specification.

fit_cfa

lavaan CFA fit object.

fit_esem

esem_fit object.

fit_besem

besem_fit object.

alignment

alignment_check result (if run).

comparison_table

Data frame of fit indices.

mplus_results

List of Mplus readModels results (if run).

Examples

data("HolzingerSwineford1939", package = "lavaan")

spec <- specify_model(
  Visual  = c("x1", "x2", "x3"),
  Textual = c("x4", "x5", "x6"),
  Speed   = c("x7", "x8", "x9"),
  data    = HolzingerSwineford1939,
  label   = "Holzinger-Swineford"
)


# Fit CFA, ESEM, and B-ESEM in R and build the comparison table
results <- run_comparison(spec, n_starts = 5L)
print(results)

# Access individual fits
summary(results$fit_esem, fit.measures = TRUE, standardized = TRUE)
std_loadings(results$fit_besem, suppress = 0.10)
factor_correlations(results$fit_esem)


## Not run: 
# Also run all three models in Mplus and compare side by side.
# Requires a licensed Mplus installation reachable via `mplus_command`.
results <- run_comparison(spec, mplus_folder = tempfile("mplus_"))
print(results)

## End(Not run)

Generate, Run, and Compare B-ESEM Invariance Models Against Mplus

Description

Creates complete Mplus .inp files for configural, weak, strong, and strict invariance, runs them via MplusAutomation, then prints a side-by-side comparison of fit statistics against the R results from esem_invariance.

Usage

run_mplus_besem_invariance(
  inv,
  output_folder,
  mplus_command = "mplus",
  group_labels = NULL,
  missing_code = -999,
  difftest = TRUE
)

Arguments

Value

A data frame (invisibly) with _R, _Mplus, and delta_ columns for each fit statistic. Printed as a table.

See Also

esem_invariance

Examples

## Not run: 
# Requires a licensed Mplus installation reachable via `mplus_command`.
inv <- esem_invariance(spec, model = "besem")
cmp <- run_mplus_besem_invariance(
  inv,
  output_folder = tempfile("mplus_inv_"),
  mplus_command = "C:/Program Files/Mplus/mplus.exe",
  group_labels  = c("1" = "MALE", "2" = "FEMALE")
)

## End(Not run)

Save Pipeline Results to CSV or xlsx

Description

Writes all standard output files to output_folder:

• <label>_fit_indices.csv – CFI/TLI/RMSEA/SRMR for all models.

• <label>_CFA_loadings.csv – standardised loading matrix (wide).

• <label>_ESEM_loadings.csv – standardised loading matrix (wide).

• <label>_BESEM_loadings.csv – standardised loading matrix (wide).

• <label>_loadings_comparison.csv – R vs Mplus loadings, long format with SEs, z-scores, p-values, and diff_std.

• <label>_omega.csv – Reliability indices (if indices is supplied).

• <label>_results.xlsx – Single xlsx workbook with all results (when xlsx = TRUE).

Usage

save_results(results, omega = NULL, output_folder, label = NULL, xlsx = FALSE)

Arguments

Value

Invisibly returns a character vector of file paths written.

Examples

data("HolzingerSwineford1939", package = "lavaan")

spec <- specify_model(
  Visual  = c("x1", "x2", "x3"),
  Textual = c("x4", "x5", "x6"),
  Speed   = c("x7", "x8", "x9"),
  data    = HolzingerSwineford1939,
  label   = "Holzinger-Swineford"
)


results <- run_comparison(spec, n_starts = 5L)
indices <- compute_indices(results)
save_results(results, indices,
             output_folder = file.path(tempdir(), "esem_results"))

Modify an Existing Target Matrix

Description

Convenience function to set specific cells of a target matrix after initial construction with make_target.

Usage

set_target(target, items, factors, value)

Arguments

Value

The modified target matrix.

Specify an ESEM Model Structure

Description

Creates a model specification object that flows automatically through all subsequent steps – CFA, ESEM, B-ESEM, alignment check, Mplus syntax, and comparison tables. Define your factor structure once at the top; everything else is derived automatically.

Usage

specify_model(
  ...,
  data,
  label = "ESEM Model",
  ordered = FALSE,
  group = NULL,
  estimator = NULL,
  missing = NULL
)

Arguments

Value

An object of class "esem_spec" containing:

factors

Named list of factor -> indicator assignments.

factor_names

Character vector of factor names.

all_items

Character vector of all indicators in order.

nfactors

Number of specific factors.

data

The supplied data frame.

label

Model label.

cfa_syntax

Ready-to-use lavaan CFA model string.

target

Target matrix for ESEM target rotation.

bifactor_target

Target matrix for B-ESEM.

Examples

data("HolzingerSwineford1939", package = "lavaan")

# Continuous indicators (MLR estimator by default)
spec <- specify_model(
  Visual  = c("x1", "x2", "x3"),
  Textual = c("x4", "x5", "x6"),
  Speed   = c("x7", "x8", "x9"),
  data    = HolzingerSwineford1939,
  label   = "Holzinger-Swineford 3-factor"
)
print(spec)

# Multi-group spec (for invariance testing)
spec_mg <- specify_model(
  Visual  = c("x1", "x2", "x3"),
  Textual = c("x4", "x5", "x6"),
  Speed   = c("x7", "x8", "x9"),
  data    = HolzingerSwineford1939,
  group   = "sex",
  label   = "HS 3-factor (multi-group)"
)


# Run the CFA/ESEM/B-ESEM comparison from a spec
results <- run_comparison(spec)
print(results)

Extract Standardized Loadings from an esem_fit

Description

Returns the standardized factor loading matrix. Uses std_rotated_loadings when present (custom WLSMV rotation path or Heywood-corrected loadings); otherwise extracts from lavaan_fit(x). Analogous to Mplus STDYX loadings.

Usage

std_loadings(x, digits = 3, suppress = 0)

Arguments

Value

A matrix of standardized loadings (items x factors).

Summary Method for esem_fit

Description

Prints a full summary for an esem_fit. For models fit through the DWLS-from-scratch path (B-ESEM ordered, or esem_ordered(method = "rotation")), the underlying lavaan_fit slot holds an auxiliary 1-factor CFA used only for weight-matrix extraction; in that case we print a self-contained bifactory summary built from the stored rotated loadings, standard errors, and WLSMV statistics. Otherwise the call is forwarded to lavaan::summary().

Usage

## S3 method for class 'esem_fit'
summary(object, fit.measures = TRUE, standardized = TRUE, rsquare = FALSE, ...)

Arguments

Value

Invisibly returns object; called for the side effect of printing the model summary.

Summary for ESEM-within-CFA Fits

Description

Forwards to lavaan::lavaan-class summary on the underlying lavaan fit, so you do not need to library(lavaan) separately.

Usage

## S3 method for class 'ewc_fit'
summary(
  object,
  fit.measures = TRUE,
  standardized = TRUE,
  show_loadings = TRUE,
  ...
)

Arguments

Value

Invisibly returns the lavaan summary object.