R/estimate_model.R

In EFAtools: Fast and Flexible Implementations of Exploratory Factor Analysis Tools

# Single estimation core: dispatch to a method fitter, then post-process once.
#
# `.PAF()`, `.ML()`, and `.ULS()` are thin fitters that return only their raw
# results (`L`, `h2`, objective `Fm`, `iter`, `convergence`, the original
# correlation matrix `orig_R`, and the matrix `R_final` whose eigenvalues are the
# final eigenvalues), plus any method-specific extras. `.finalize_fit()` performs
# the shared post-processing common to all methods, and `.estimate_model()`
# assembles the method-specific output object.

# Lower bound the ML and ULS optimisers impose on the uniquenesses. A solution
# whose uniqueness is pinned at this floor is an improper (boundary) solution:
# the unconstrained optimum would place it at or below zero. Shared by the ML/ULS
# fitters and the Heywood detector in .finalize_fit() so the value has one source.
.uniqueness_floor <- 0.005

# Squared multiple correlations (1 - 1/diag(R^-1)) used as starting communalities by the
# estimators. Shared by the PAF, ML, and ULS fitters so the start has one source.
#
# A squared multiple correlation is a squared correlation and so lies in [0, 1], but the
# expression only respects that bound for a positive definite R: on an indefinite or
# ill-conditioned matrix (an unsmoothed bootstrap resample, say) it escapes the range and
# would otherwise be written straight onto a correlation diagonal or handed to the
# optimiser as a start outside its box. Clamp it, as psych::smc() does.
#
# The inverse diagonal comes from a Cholesky factorisation, which uses the symmetry and
# positive definiteness that solve()'s LU ignores. Two cases fall back to solve(): a
# matrix the Cholesky rejects (not positive definite), and one whose implied uniqueness
# 1/diag(R^-1) has underflowed to p^2 * eps, where the Cholesky still succeeds but the
# matrix is singular to working precision. Both then reach solve()'s own singularity
# error, which callers rely on -- .parallel_sim_eig() rejects a simulated draw on it.
# The p^2 factor is deliberately more generous than the p * eps rank tolerance used on
# eigenvalue ratios elsewhere (.prepare_cor_input()): the two quantities are not on the
# same scale, and the asymmetry favours caution, because a needless fallback only repeats
# the calculation solve() would have done anyway while a missed one turns a rejected draw
# into an accepted degenerate one.
#
# The result is unnamed on both routes (chol2inv() drops dimnames, solve() keeps them), so
# the return shape does not depend on which one ran. Callers that surface the values name
# them from the correlation matrix themselves.
#
# `R_inv` lets a caller that has already formed R^-1 hand it over rather than pay for a
# second factorisation (efa_screen() builds one inverse for several measures at once). The
# underflow check and the solve() fallback below still run, so the singularity contract is
# the same on either route.
.smc_start <- function(R, R_inv = NULL) {
  d <- if (!is.null(R_inv)) {
    diag(R_inv)
  } else {
    tryCatch(diag(chol2inv(chol(R))), error = function(e) NULL)
  }
  if (is.null(d) || min(1 / d) <= ncol(R)^2 * .Machine$double.eps) {
    d <- diag(solve(R))
  }
  unname(pmin(pmax(1 - 1 / d, 0), 1))
}

# Sign convention shared by the unrotated (.finalize_fit) and rotated (.reflect_and_order)
# solutions: reflect each factor so its loadings sum to a non-negative value (as in the
# psych package and SPSS). Returns the +/-1 sign vector; apply it with
# `L %*% diag(signs, nrow = length(signs))` -- the `nrow` guards the single-factor case,
# where diag() would otherwise read a length-1 vector as a matrix dimension.
.reflect_signs <- function(L) {
  signs <- sign(colSums(L))
  signs[signs == 0] <- 1
  signs
}

# Shared post-processing for an unrotated solution. Reflects the loadings to a
# consistent sign, names them (with a V-fallback when the input is unnamed),
# computes the explained variances, fit indices, communalities, model-implied
# correlation matrix, residuals, and the original/final eigenvalues.
.finalize_fit <- function(fit, N, method, lean = FALSE, fiml = NULL) {

  L <- fit$L
  orig_R <- fit$orig_R
  h2 <- fit$h2
  n_factors <- ncol(L)

  # reverse the sign of loadings as done in the psych package and SPSS
  signs <- .reflect_signs(L)
  L <- L %*% diag(signs, nrow = length(signs))

  # Bootstrap replicate path: .boot_se_ci() aggregates only each replicate's
  # unrotated loadings, fit indices, and residuals, so only those are computed.
  # The two eigendecompositions, variable naming, explained variances, and
  # Heywood detection are skipped (none are aggregated per replicate), and the
  # analytic RMSEA confidence bounds are not solved (.gof(ci = FALSE)). The sign
  # reflection above is kept so the loadings are identical to the full path,
  # leaving the downstream target alignment unchanged.
  if (lean) {
    model_implied_R <- L %*% t(L) + diag(1 - h2)
    return(list(
      unrot_loadings = L,
      fit_indices = .gof(L, orig_R, N, method, fit$Fm, ci = FALSE, fiml = fiml),
      residuals = orig_R - model_implied_R,
      convergence = fit$convergence
    ))
  }

  if (!is.null(colnames(orig_R))) {
    # name the loading matrix so the variables can be identified
    rownames(L) <- colnames(orig_R)
  } else {
    varnames <- paste0("V", seq_len(ncol(orig_R)))
    colnames(orig_R) <- varnames
    rownames(orig_R) <- varnames
    rownames(L) <- varnames
  }

  colnames(L) <- paste0("F", seq_len(n_factors))

  vars_accounted <- .compute_vars(L_unrot = L, L_rot = L)
  colnames(vars_accounted) <- colnames(L)

  fit_ind <- .gof(L, orig_R, N, method, fit$Fm, fiml = fiml)

  # calculate model implied R
  model_implied_R <- L %*% t(L) + diag(1 - h2)

  # create the output object
  class(L) <- c("efa_loadings", "LOADINGS")

  # Name communalities
  names(h2) <- colnames(orig_R)

  # Detect Heywood (improper) cases. Named integer vector of the affected
  # variables (empty if none); surfaced to the user by EFA() and shown in
  # summary(). Under PAF the communality can reach or exceed 1 directly. Under
  # ML/ULS the optimiser constrains the uniquenesses to [floor, 1], so an improper
  # solution instead shows up as a uniqueness pinned at the lower floor (the
  # boundary case); flag those too so detection is consistent across estimators.
  # The ML/ULS fitters return psi as a p x 1 matrix; `|` propagates a dim
  # attribute in preference to names, which would strip the variable names off
  # the result, so drop the dim with as.vector() before combining.
  heywood_comm <- h2 >= 1
  heywood_boundary <- if (!is.null(fit$psi)) {
    as.vector(fit$psi) <= .uniqueness_floor + sqrt(.Machine$double.eps)
  } else {
    rep(FALSE, length(h2))
  }
  heywood <- which(heywood_comm | heywood_boundary)

  list(
    orig_R = orig_R,
    h2 = h2,
    # only the eigenvalues are reported, so skip the eigenvectors LAPACK would
    # otherwise be asked for
    orig_eigen = eigen(orig_R, symmetric = TRUE, only.values = TRUE)$values,
    final_eigen = eigen(fit$R_final, symmetric = TRUE, only.values = TRUE)$values,
    iter = fit$iter,
    convergence = fit$convergence,
    heywood = heywood,
    unrot_loadings = L,
    vars_accounted = vars_accounted,
    fit_indices = fit_ind,
    model_implied_R = model_implied_R,
    residuals = orig_R - model_implied_R
  )
}

# Dispatch to the requested fitter, run the shared post-processor, and assemble
# the method-specific output object (field set and order differ per method).
.estimate_model <- function(R, method, n_factors, N = NA,
                            type = "none", max_iter = NA, init_comm = NA,
                            criterion = NA, criterion_type = NA, abs_eigen = NA,
                            start_method = NA, weights = NULL, lean = FALSE,
                            fiml = NULL) {

  fit <- switch(
    method,
    PAF = .PAF(R, n_factors = n_factors, type = type, max_iter = max_iter,
               init_comm = init_comm, criterion = criterion,
               criterion_type = criterion_type, abs_eigen = abs_eigen),
    ML = .ML(R, n_factors = n_factors, start_method = start_method),
    ULS = .ULS(R, n_factors = n_factors),
    DWLS = .DWLS(R, n_factors = n_factors, weights = weights)
  )

  common <- .finalize_fit(fit, N = N, method = method, lean = lean, fiml = fiml)

  # The bootstrap replicate path needs only the post-processed common fields;
  # skip the method-specific output assembly, which the aggregation never reads.
  if (lean) {
    return(common)
  }

  if (method == "PAF") {

    h2_init <- fit$h2_init
    names(h2_init) <- colnames(common$orig_R)

    output <- list(
      orig_R = common$orig_R,
      h2_init = h2_init,
      h2 = common$h2,
      orig_eigen = common$orig_eigen,
      init_eigen = fit$init_eigen,
      final_eigen = common$final_eigen,
      iter = common$iter,
      convergence = common$convergence,
      heywood = common$heywood,
      unrot_loadings = common$unrot_loadings,
      vars_accounted = common$vars_accounted,
      fit_indices = common$fit_indices,
      model_implied_R = common$model_implied_R,
      residuals = common$residuals,
      settings = fit$settings
    )

  } else if (method == "ML") {

    output <- c(common, list(settings = fit$settings))

  } else {

    output <- common

  }

  output
}