OptimalFactor

A comprehensive package for optimal factor model refinement in both EFA and CFA using machine-learning inspired boosting algorithms.
Visit the author’s website

Features

• EFA-Boosting Algorithm: Advanced iterative optimization for Exploratory Factor Analysis with smart-pruning of the candidate space (large speed-ups on long instruments) and a canonical stop_reason field that explains exactly why the loop ended
• Adaptive Fit Indices: Dynamic weights based on df x N following Kenny & McCoach (2003) and Shi, Lee & Maydeu-Olivares (2019)
• Automatic Problem Detection: Heywood cases, cross-loadings, and low loadings
• Admissibility over fit (CFA): cfa_boosting() removes items below the loading floor, and items loading on a foreign factor, even when the fit targets are already met — because global fit does not reveal either of them
• Monte Carlo validation: simulate_recovery() and simulate_cfa_recovery() measure how often each pipeline recovers a structure that is known in advance
• Global Search: Multi-item removal optimization
• AI Integration: Optional GPT-powered conceptual analysis of removed items
• Specification Search (CFA): Heuristic hill-climbing over seed configurations following MacCallum (1986), with move / drop / cov operations and optional bifactor variant

Looking for the interactive app? The OptimalFactor Wizard (guided 5-phase Shiny interface with reliability, convergent/discriminant validity, AI autopilot and APA-7 report downloads) is distributed separately as a web application on Posit Connect Cloud. This package contains the analysis engine only.

Installation

You can install the latest version of OptimalFactor from GitHub:

if (!require("devtools")) {
  install.packages("devtools")
}
devtools::install_github("jventural/OptimalFactor")

Quick Start

EFA-Boosting

library(OptimalFactor)

result <- efa_boosting(
  data = your_data,
  name_items = "item",
  n_factors = 3,
  verbose = TRUE
)

Main Functions

Function	Description
efa_boosting()	EFA optimization with adaptive composite fit
cfa_boosting()	CFA optimization with modification indices
specification_search_theory()	Theory-guided CFA specification search (MacCallum, 1986 + theory-congruence loss)
cross_validate_cfa()	Split-half cross-validation of a factor model
bifactor_indices()	Bifactor statistical indices
redundancy_short_form()	Redundancy-guided short form of a unidimensional scale
item_stability()	Resampling stability of the EFA-boosting item selection
simulate_recovery()	Monte Carlo recovery of a known factor structure (EFA)
simulate_cfa_recovery()	Monte Carlo item selection against a known population (CFA)
report_efa_results()	Structured + console report of an EFA-boosting run
report_cfa_results()	Structured + console report of a CFA-boosting run
print_conceptual_analysis()	Display AI-generated item analyses

Is the item selection stable, or is it chance?

Data-driven item purification invites the classic objection of capitalization on chance (MacCallum, Roznowski & Necowitz, 1992). item_stability() answers it empirically: it re-runs the whole pipeline over resamples and reports how often each item survives and how stable its factor assignment is.

st <- item_stability(your_data, "item", n_factors = 3, R = 100)
st              # retention rates, factor agreement, unstable decisions

An item dropped in 97 of 100 resamples is a defensible removal; one dropped in 55 is a coin flip. The printed summary flags every decision taken in 25–75 % of resamples as unstable.

When does the pipeline recover the true structure?

simulate_recovery() simulates ordinal data from a population model with good, cross-loading and weak items, runs the pipeline on each replication and crosses conditions over sample size, loading size and items per factor.

sim <- simulate_recovery(n = c(200, 500, 1000), loading = c(0.50, 0.70),
                         n_reps = 200, n_cores = 8)
sim$summary     # recovery rate, sensitivity, specificity per condition
plot(sim)       # the recovery curve against sample size

Sensitivity (contaminated items removed) and specificity (good items kept) must be read together: an algorithm that removes everything scores a perfect sensitivity and a dismal specificity.

Replications are fitted in parallel through n_cores, while the datasets are always drawn sequentially from seed, so the results do not depend on how many cores ran them. Each replication is capped by timeout (120 s by default): resampling puts the pipeline on datasets nobody would analyse by hand, and a single one of them can otherwise grind for minutes and stall the whole run.

What if the number of factors is wrong?

Every condition above hands the algorithm the true number of factors, and in practice nobody knows it. n_factors_fitted crosses that too, and the answer is asymmetric:

simulate_recovery(n = c(300, 600), n_factors = 3,
                  n_factors_fitted = c(2, 3, 4), n_reps = 200)

Extracting too few factors is destructive: specificity falls to around .65–.73, so good items are discarded to force the data into a structure that cannot hold them, and the loop stops satisfied with the result. Extracting too many makes the pipeline inert instead: it removes nothing, because almost any removal would breach min_items_per_factor, and the scale comes out untouched. When the dimensionality is uncertain, erring towards more factors is the harmless mistake.

The confirmatory side

simulate_cfa_recovery() asks a different question, and the difference is not cosmetic. In the exploratory case the structure is discovered, so recovery means every retained item landed on its true factor. Here the structure is imposed by the researcher, so the assignment cannot be wrong: what is measured is the item selection. The specified model deliberately contains the contaminated items.

sim <- simulate_cfa_recovery(n = c(300, 600), loading = 0.60, n_reps = 200)
sim             # exact_rate, sensitivity, specificity, covariances added

mean_covs_added counts the residual covariances the pipeline introduced. The population has none, so anything above zero is capitalization on chance turned into a number.

This simulation was written to evaluate cfa_boosting() and found two defects in it, both now fixed. The loop tested the fit targets first and stopped as soon as they were met, leaving the low-loading check unreachable: an item loading .20 in the population sat comfortably inside a model with RMSEA = .03. And global fit does not betray a cross-loading either, because the omitted path is absorbed by the interfactor correlation — six items cross-loading at .60 left RMSEA at .053 and pushed CFI up to .997. Enforcing both criteria as admissibility rather than fit took sensitivity from .167 to .958 and exact recovery from 0 to .917, with specificity rising to 1.000 and fewer residual covariances needed.

# Both are on by default; set either to FALSE for the pre-1.3.0 behaviour
cfa_boosting(your_data, model,
             thresholds = list(enforce_loading = TRUE,
                               enforce_simple_structure = TRUE,
                               cross_loading = 0.30))

Cross-loading detection uses the modification index of the absent loading, but never significance alone: the standardized expected parameter change must also reach cross_loading, since a significant index with a trivial EPC is precisely the capitalization on chance the whole exercise exists to avoid.

Reliability floor

Removing items to improve fit tends to cost internal consistency. Setting min_omega makes that trade-off explicit: an item is removed only if every factor keeps McDonald’s omega at or above the floor.

efa_boosting(your_data, "item", n_factors = 3,
             thresholds = list(min_omega = 0.80))

The floor is a constraint, not a term in the loss — reliability is not interpretable in the ill-fitting intermediate models the loss visits, so omega vetoes individual removals instead of trading against the fit indices. Omega per factor (omega_final) is now reported in every run, floor or no floor.

Theory-Guided Specification Search (CFA)

specification_search_theory() runs a heuristic hill-climbing search over CFA model configurations in the spirit of MacCallum (1986), extended with a theory-congruence term in the loss: candidate models are penalized for moving items away from their theoretical factor or for dropping theoretical items, so the search cannot drift toward models that fit well but break the intended structure. A single theory_weight parameter grades how much theory counts relative to fit (0 reproduces the classic fit-only search).

For each seed (a candidate number of factors with an initial item-to-factor assignment — the theoretical structure itself is always included as a seed) the algorithm iteratively tries three local operations and keeps the change with the lowest composite loss:

• move an item from its current factor to another;
• drop an item with low loading or causing strain;
• cov add a residual covariance suggested by modification indices.

A bifactor variant (one orthogonal general factor plus the configured group factors) can be tried automatically for every seed with k >= 2. The procedure returns every model evaluated, the subset that meets the user-supplied fit targets (CFI / RMSEA), and the best model under the composite loss.

library(lavaan)
library(OptimalFactor)

# Simulated 3-factor data
sim <- '
  F1 =~ 0.75*x1 + 0.70*x2 + 0.65*x3
  F2 =~ 0.75*x4 + 0.70*x5 + 0.65*x6
  F3 =~ 0.75*x7 + 0.70*x8 + 0.65*x9
  F1 ~~ 0.30*F2 + 0.30*F3
  F2 ~~ 0.30*F3
'
df    <- simulateData(sim, sample.nobs = 500, seed = 2026)
items <- paste0("x", 1:9)

res <- specification_search_theory(
  data          = df,
  items         = items,
  theory        = list(F1 = c("x1","x2","x3"),
                        F2 = c("x4","x5","x6"),
                        F3 = c("x7","x8","x9")),
  theory_weight = 0.5,   # 0 = fit-only (classic search); higher = more conservative
  estimator     = "ML",
  ordered       = FALSE,
  try_bifactor  = TRUE,
  verbose       = TRUE
)

print(res, top = 8)         # ranking by composite loss (+ theory congruence)
res$successful              # configurations meeting CFI / RMSEA targets
summary(res$best$fit)       # fitted lavaan object for the best model

Deprecated: the original fit-only specification_search() is kept for backward compatibility but emits a deprecation warning; use specification_search_theory() (with theory_weight = 0 if you really want the fit-only behaviour).

Important caveat — MacCallum (1986). Specification search capitalizes on chance. Use it as an exploratory device only and:

1. Report the procedure transparently as exploratory.
2. Cross-validate the resulting model with an independent sample or via bootstrap.
3. Justify every accepted modification on substantive theoretical grounds.

The function prints this warning at the start of every run (suppress with verbose = FALSE).

Examples

Basic tutorial of the R OptimalFactor library

References

• MacCallum, R. C. (1986). Specification searches in covariance structure modeling. Psychological Bulletin, 100(1), 107–120.
• Saris, W. E., Satorra, A., & van der Veld, W. M. (2009). Testing structural equation models or detection of misspecifications? Structural Equation Modeling, 16(4), 561–582.
• Kenny, D. A., & McCoach, D. B. (2003). Effect of the number of variables on measures of fit in structural equation modeling. Structural Equation Modeling, 10(3), 333-351. https://doi.org/10.1207/S15328007SEM1003_1
• Shi, D., Lee, T., & Maydeu-Olivares, A. (2019). Understanding the model size effect on SEM fit indices. Educational and Psychological Measurement, 79(2), 310-334. https://doi.org/10.1177/0013164418783530

License

GPL-3

Citation

Ventura-Leon, J. (2026). OptimalFactor: Item Selection in Factor Analysis by Boosting and Specification Search [R package]. GitHub. https://github.com/jventural/OptimalFactor

Author

Jose Ventura-Leon jventuraleon@gmail.com