R/efa_retention.R

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

# Unified result class for the factor-retention criteria. Every criterion
# returns an `efa_retention` object built by `.new_efa_retention()`, printed by
# one `print.efa_retention()`, and plotted by one `plot.efa_retention()` that
# dispatches to the two ggplot helpers below. The registry is the single source
# of truth mapping each criterion id to its display label.

# Registry of the factor-retention criteria (id -> metadata). Each `fun` runs
# its criterion from the efa_retain() control list `ctl`; `x` is the raw data
# when `needs_raw`, the prepared correlation matrix otherwise. Every criterion
# that fits a model through efa_fit() (HULL, KGC, PARALLEL, SCREE, NEST, SMT)
# receives `ctl$estimate_control`, so the estimation settings reach all of them
# and not just some; those that also forward arguments to efa_fit() receive them
# via `ctl$dots`. CD, EKC and MAP run no efa_fit() model (CD's `max_iter` caps
# its own comparison-data generation, not an EFA), so neither applies to them.
# `poly_ok = FALSE` marks the criteria that do not support polychoric/tetrachoric
# correlations -- either because they compare the data against continuous
# reference data (CD, PARALLEL, NEST, HULL) or because their normal-theory
# chi-square test is not valid for such correlations (SMT); these are skipped
# with an informative note by efa_retain().
.retention_registry <- list(
  CD = list(
    label = "Comparison data", needs_raw = TRUE, poly_ok = FALSE,
    fun = function(x, ctl) {
      efa_cd(x, n_factors_max = ctl$n_factors_max, N_pop = ctl$N_pop,
         N_samples = ctl$N_samples, alpha = ctl$alpha,
         cor_method = ctl$cor_method, max_iter = ctl$max_iter_CD)
    }),
  EKC = list(
    label = "Empirical Kaiser Criterion", needs_raw = FALSE,
    fun = function(x, ctl) {
      efa_ekc(x, N = ctl$N, use = ctl$use, cor_method = ctl$cor_method)
    }),
  HULL = list(
    label = "Hull method", needs_raw = FALSE, poly_ok = FALSE,
    fun = function(x, ctl) {
      do.call(efa_hull, c(list(x, N = ctl$N, n_fac_theor = ctl$n_fac_theor,
                           estimator = ctl$estimator, gof = ctl$gof,
                           eigen_type = ctl$eigen_type_HULL, use = ctl$use,
                           cor_method = ctl$cor_method,
                           n_datasets = ctl$n_datasets, percent = ctl$percent,
                           decision_rule = ctl$decision_rule,
                           n_factors = ctl$n_factors,
                           estimate_control = ctl$estimate_control),
                      ctl$dots))
    }),
  KGC = list(
    label = "Kaiser-Guttman criterion", needs_raw = FALSE,
    fun = function(x, ctl) {
      do.call(efa_kgc, c(list(x, eigen_type = ctl$eigen_type_other, use = ctl$use,
                          cor_method = ctl$cor_method,
                          n_factors = ctl$n_factors, estimator = ctl$estimator,
                          estimate_control = ctl$estimate_control),
                     ctl$dots))
    }),
  MAP = list(
    label = "Minimum average partial", needs_raw = FALSE,
    fun = function(x, ctl) {
      efa_map(x, use = ctl$use, cor_method = ctl$cor_method)
    }),
  NEST = list(
    label = "Next Eigenvalue Sufficiency Test", needs_raw = FALSE, poly_ok = FALSE,
    fun = function(x, ctl) {
      do.call(efa_nest, c(list(x, N = ctl$N, use = ctl$use,
                           cor_method = ctl$cor_method, alpha = ctl$alpha_nest,
                           n_datasets = ctl$n_datasets_nest, estimator = ctl$estimator,
                           estimate_control = ctl$estimate_control),
                      ctl$dots))
    }),
  PARALLEL = list(
    label = "Parallel analysis", needs_raw = FALSE, poly_ok = FALSE,
    fun = function(x, ctl) {
      do.call(efa_parallel, c(list(x, N = ctl$N, n_datasets = ctl$n_datasets,
                               percent = ctl$percent,
                               eigen_type = ctl$eigen_type_other, use = ctl$use,
                               cor_method = ctl$cor_method,
                               decision_rule = ctl$decision_rule,
                               n_factors = ctl$n_factors, estimator = ctl$estimator,
                               estimate_control = ctl$estimate_control),
                          ctl$dots))
    }),
  SCREE = list(
    label = "Scree plot", needs_raw = FALSE, visual = TRUE,
    fun = function(x, ctl) {
      do.call(efa_scree, c(list(x, eigen_type = ctl$eigen_type_other, use = ctl$use,
                            cor_method = ctl$cor_method,
                            n_factors = ctl$n_factors, estimator = ctl$estimator,
                            estimate_control = ctl$estimate_control),
                       ctl$dots))
    }),
  SMT = list(
    label = "Sequential model tests", needs_raw = FALSE, poly_ok = FALSE,
    fun = function(x, ctl) {
      efa_smt(x, N = ctl$N, use = ctl$use, cor_method = ctl$cor_method,
          estimate_control = ctl$estimate_control)
    })
)

# Control list consumed by the factor-retention registry funs. The defaults
# mirror the efa_retain() argument defaults so any requested criterion resolves;
# callers pass only what they override. `gof` defaults to the Hull
# goodness-of-fit indices valid for `estimator` (PAF supports only the CAF).
.n_factors_ctl <- function(N = NA, use = "pairwise.complete.obs",
                           cor_method = "pearson", n_factors_max = NA,
                           N_pop = 10000, N_samples = 500, alpha = .30,
                           max_iter_CD = 50, n_fac_theor = NA, estimator = "ML",
                           gof = if (estimator == "PAF") "CAF" else c("CAF", "CFI", "RMSEA"),
                           eigen_type_HULL = "SMC", eigen_type_other = "SMC",
                           n_factors = 1, n_datasets = 1000, percent = 95,
                           decision_rule = "means",
                           n_datasets_nest = 1000, alpha_nest = .05,
                           estimate_control = NULL, dots = list()) {
  list(N = N, use = use, cor_method = cor_method,
       n_factors_max = n_factors_max, N_pop = N_pop,
       N_samples = N_samples, alpha = alpha, max_iter_CD = max_iter_CD,
       n_fac_theor = n_fac_theor, estimator = estimator, gof = gof,
       eigen_type_HULL = eigen_type_HULL,
       eigen_type_other = eigen_type_other, n_factors = n_factors,
       n_datasets = n_datasets, percent = percent,
       decision_rule = decision_rule,
       n_datasets_nest = n_datasets_nest, alpha_nest = alpha_nest,
       estimate_control = estimate_control, dots = dots)
}

# One sample-size rule for the criteria that need more observations than
# variables: parallel analysis and NEST draw reference data with N cases, and
# EKC and SMT rest on normal-theory quantities ((1 + sqrt(J / N))^2 and the
# Bartlett-corrected chi-square) that are not defined below that boundary. Kept
# in one place so the contract, the message, and the condition class stay
# identical across the four criteria. Both arguments must already be resolved:
# call it after the N-required guard (.prepare_cor_input(N_policy = "required"),
# or the caller's own check), because an NA would make the comparison itself fail
# instead of raising the classed condition.
.assert_n_gt_vars <- function(N, n_vars, error_call = rlang::caller_env()) {
  if (N <= n_vars) {
    cli::cli_abort(
      c("{.arg N} must be larger than the number of variables.",
        "x" = "You supplied {.arg N} = {N} for {n_vars} variable{?s}."),
      class = "efa_n_too_small",
      call = error_call
    )
  }
  invisible(N)
}

# Name a criterion's suggested factor counts the way efa_retain() aggregates
# them: a single-variant suggestion keeps the bare id, a multi-variant one
# becomes "<id>_<variant>".
.retention_key <- function(id, nf) {
  names(nf) <- ifelse(names(nf) == id, id, paste(id, names(nf), sep = "_"))
  nf
}

# Construct an efa_retention object from a list of per-sub-variant records. The
# top-level `n_factors` named vector is derived from the records. `subtitle` is an
# optional one-line context string (e.g. the estimation method); `note` is an
# optional vector of cli info lines.
.new_efa_retention <- function(id, results, settings,
                               subtitle = NULL, note = NULL) {

  if (!id %in% names(.retention_registry)) {
    cli::cli_abort("Unknown factor-retention criterion id {.val {id}}.",
                   class = "efa_unknown_criterion")
  }

  label <- .retention_registry[[id]]$label

  n_factors <- vapply(results, function(r) as.numeric(r$n_factors), numeric(1))
  names(n_factors) <- vapply(results, function(r) r$name, character(1))

  structure(
    list(
      criterion = c(id = id, label = label),
      n_factors = n_factors,
      results = results,
      subtitle = subtitle,
      note = note,
      settings = settings
    ),
    class = "efa_retention"
  )
}

# Subtitle naming the eigenvalue types a criterion computed its eigenvalues on,
# shared by the eigenvalue-based criteria so they all say it the same way. `detail`
# appends a clause (e.g. the number of simulated reference datasets) in place of the
# closing full stop.
.eigen_subtitle <- function(eigen_type, detail = NULL) {
  paste0("Eigenvalues found using ", cli::ansi_collapse(eigen_type),
         if (is.null(detail)) "." else paste0("; ", detail, "."))
}

# Render a whole-number count (of factors, or of simulated datasets) for the report, a
# condition message, or a plot label. The criteria differ in the storage mode they put in
# their record: a count derived with cumprod() or which.min() - 1 comes out double, one
# taken from an index integer. Under a negative options(scipen) a double is rendered in
# scientific notation ("3e+00" instead of "3") by as.character() and paste0() alike, and
# also by grid, which coerces a numeric plot label at draw time. A count is exempt from
# that option because it is whole by construction; a genuinely continuous quantity (an
# eigenvalue, a fit index) still honours it.
#
# scientific = FALSE pins the fixed notation whatever the option is set to. trim keeps a
# vector of counts unpadded, so that the point labels of a hull with ten or more solutions
# are not indented to the width of the widest of them.
#
# This covers the counts the criteria compute. A count that cli interpolates from a user
# argument -- `{N}` in .assert_n_gt_vars() below -- does not come through here.
.retention_count <- function(n) {
  format(n, scientific = FALSE, trim = TRUE)
}

# Bullet lines (one per record) shared by format.efa_retention and
# format.efa_retain.
.retention_bullets <- function(results) {
  vapply(results, function(r) {
    value <- if (is.na(r$n_factors)) "not applicable" else .retention_count(r$n_factors)
    paste0(r$label, ": ", value)
  }, character(1))
}

#' Print method for efa_retention objects
#'
#' @param x an object of class efa_retention, returned by a factor-retention
#'   criterion (e.g. [efa_ekc()] or [efa_hull()]).
#' @param ... not used.
#'
#' @returns `print()` returns its argument `x` invisibly; it is
#'   `cat(format(x), sep = "\n")`.
#'
#' @export
#' @method print efa_retention
#'
#' @examples
#' efa_ekc(test_models$baseline$cormat, N = 500)
print.efa_retention <- function(x, ...) {
  cat(format(x, ...), sep = "\n")
  invisible(x)
}

#' Format method for efa_retention objects
#'
#' @param x an object of class efa_retention, returned by a factor-retention
#'   criterion (e.g. [efa_ekc()] or [efa_hull()]).
#' @param ... not used.
#'
#' @returns A character vector with the report lines (styled to the active
#'   console theme; plain when colours are disabled).
#'
#' @export
#' @method format efa_retention
#'
#' @examples
#' writeLines(format(efa_ekc(test_models$baseline$cormat, N = 500)))
format.efa_retention <- function(x, ...) {
  cli::cli_format_method({
    cli::cli_rule(left = x$criterion[["label"]])
    if (!is.null(x$subtitle)) {
      cli::cli_text("{x$subtitle}")
    }
    # criteria with no numeric suggestion (e.g. the visual scree plot) skip the
    # bullets and rely on their subtitle/note
    if (any(!is.na(x$n_factors))) {
      cli::cli_text("")
      cli::cli_ul(.retention_bullets(x$results))
    }
    if (!is.null(x$note)) {
      cli::cli_text("")
      for (msg in x$note) {
        cli::cli_alert_info("{msg}", wrap = TRUE)
      }
    }
  })
}

#' Plot method for efa_retention objects
#'
#' Plots the result of a factor-retention criterion. Eigenvalue-based criteria
#' (e.g. [efa_ekc()]) are shown as an eigenvalue plot, the Hull method ([efa_hull()]) as
#' a convex-hull plot. Criteria with more than one sub-variant are faceted.
#'
#' @param x an object of class efa_retention, returned by a factor-retention
#'   criterion (e.g. [efa_ekc()] or [efa_hull()]).
#' @param ... not used.
#'
#' @returns A [ggplot2::ggplot] object, or invisibly `NULL` if the criterion has
#'   no plottable result.
#'
#' @export
#' @method plot efa_retention
#'
#' @examples
#' plot(efa_ekc(test_models$baseline$cormat, N = 500))
plot.efa_retention <- function(x, ...) {

  plot_types <- unique(vapply(x$results, function(r) r$plot_type, character(1)))
  plot_types <- setdiff(plot_types, "none")

  if (length(plot_types) == 0) {
    cli::cli_inform("No plot is available for {x$criterion[['label']]}.",
                    class = "efa_no_plot")
    return(invisible(NULL))
  }

  if ("hull" %in% plot_types) {
    .gg_hull_plot(x)
  } else {
    .gg_eigen_plot(x)
  }

}

# Eigenvalue plot (factor index on x, eigenvalues on y) with an optional dashed
# reference series, an optional horizontal threshold, and an optional highlighted
# retained point. Covers the EKC/KGC/PARALLEL/SCREE/CD plots. Returns a ggplot.
#' @importFrom rlang .data
.gg_eigen_plot <- function(x) {

  # drop records with no plottable points (e.g. CD when it suggests 0 factors)
  records <- Filter(function(r) length(r$x) > 0, x$results)
  if (length(records) == 0) {
    cli::cli_inform("No plot is available for {x$criterion[['label']]}.",
                    class = "efa_no_plot")
    return(invisible(NULL))
  }

  variant_levels <- vapply(records, function(r) r$label, character(1))

  dat <- do.call(rbind, lapply(records, function(r) {
    data.frame(
      variant = r$label,
      factor = r$x,
      # a record may have no primary series (e.g. PARALLEL without real data
      # plots only its reference series)
      value = if (is.null(r$y)) NA_real_ else r$y,
      reference = if (is.null(r$reference)) NA_real_ else r$reference,
      stringsAsFactors = FALSE
    )
  }))
  dat$variant <- factor(dat$variant, levels = variant_levels)

  highlights <- do.call(rbind, lapply(records, function(r) {
    if (is.null(r$highlight) || is.na(r$highlight) || r$highlight < 1) return(NULL)
    # `factor` positions the point, `label` names it (see .retention_count())
    data.frame(variant = r$label, factor = r$highlight,
               value = r$y[r$highlight], label = .retention_count(r$highlight),
               stringsAsFactors = FALSE)
  }))

  thresholds <- do.call(rbind, lapply(records, function(r) {
    if (is.null(r$threshold)) return(NULL)
    data.frame(variant = r$label, yintercept = r$threshold,
               stringsAsFactors = FALSE)
  }))

  p <- ggplot2::ggplot(dat, ggplot2::aes(.data$factor, .data$value)) +
    ggplot2::geom_line(na.rm = TRUE) +
    ggplot2::geom_point(na.rm = TRUE)

  if (any(!is.na(dat$reference))) {
    p <- p + ggplot2::geom_line(ggplot2::aes(y = .data$reference),
                                linetype = 2, colour = "darkgray", na.rm = TRUE)
  }

  if (!is.null(thresholds)) {
    thresholds$variant <- factor(thresholds$variant, levels = variant_levels)
    p <- p + ggplot2::geom_hline(data = thresholds,
                                 ggplot2::aes(yintercept = .data$yintercept),
                                 linetype = 2, colour = "darkgray")
  }

  # records can carry several named reference series (e.g. PARALLEL's simulated
  # means and percentile eigenvalues), each drawn as a dashed coloured line with
  # a shared legend; criteria without `references` are unaffected
  ref_records <- Filter(function(r) !is.null(r$references), records)
  if (length(ref_records) > 0) {
    ref_dat <- do.call(rbind, lapply(ref_records, function(r) {
      data.frame(variant = r$label,
                 factor = rep(r$x, length(r$references)),
                 series = rep(names(r$references), each = length(r$x)),
                 value = unlist(r$references, use.names = FALSE),
                 stringsAsFactors = FALSE)
    }))
    ref_dat$variant <- factor(ref_dat$variant, levels = variant_levels)
    ref_dat$series <- factor(ref_dat$series, levels = unique(ref_dat$series))
    # label the primary series in the legend alongside the reference series
    # (drawn over the identical primary line, so only the legend key is added)
    ref_labels <- vapply(ref_records, function(r) r$label, character(1))
    real_dat <- dat[dat$variant %in% ref_labels & !is.na(dat$value), ,
                    drop = FALSE]
    if (nrow(real_dat) > 0) {
      p <- p +
        ggplot2::geom_line(data = real_dat,
                           ggplot2::aes(linetype = "Real Eigenvalues")) +
        ggplot2::scale_linetype_manual(values = c(`Real Eigenvalues` = 1))
    }
    p <- p +
      ggplot2::geom_line(data = ref_dat,
                         ggplot2::aes(.data$factor, .data$value,
                                      colour = .data$series),
                         linetype = 2, na.rm = TRUE) +
      ggplot2::scale_colour_viridis_d(end = 0.8) +
      ggplot2::labs(colour = NULL, linetype = NULL)
  }

  if (!is.null(highlights)) {
    highlights$variant <- factor(highlights$variant, levels = variant_levels)
    p <- p +
      ggplot2::geom_point(data = highlights, shape = 1, size = 4, colour = "red") +
      ggplot2::geom_text(data = highlights, ggplot2::aes(label = .data$label),
                         colour = "red", vjust = -1)
  }

  if (length(variant_levels) > 1) {
    p <- p + ggplot2::facet_wrap(ggplot2::vars(.data$variant), scales = "free_y")
  }

  # y-axis label defaults to "Eigenvalues" but a criterion can override it (e.g.
  # CD plots mean RMSE of the eigenvalues)
  y_label <- records[[1]]$y_label
  if (is.null(y_label)) y_label <- "Eigenvalues"

  p +
    ggplot2::scale_x_continuous(breaks = seq_len(max(dat$factor))) +
    ggplot2::labs(x = "Factor", y = y_label, title = x$criterion[["label"]]) +
    .gg_theme()

}

# Convex-hull plot (degrees of freedom on x, goodness-of-fit on y). Points on the
# hull are emphasised and connected; the retained solution is highlighted in red;
# each point is labelled with its number of factors. Records whose hull is empty
# (the degenerate fewer-than-three-solutions case) are dropped. Returns a ggplot.
#' @importFrom rlang .data
.gg_hull_plot <- function(x) {

  results <- Filter(function(r) any(r$on_hull), x$results)

  if (length(results) == 0) {
    cli::cli_inform("No hull plot is available for {x$criterion[['label']]}.")
    return(invisible(NULL))
  }

  # The RMSEA record carries 1 - RMSEA on the y-axis (see the HULL details), so its
  # facet is labelled accordingly; the record/bullet label stays "RMSEA".
  hull_facet <- function(label) ifelse(label == "RMSEA", "1 - RMSEA", label)

  variant_levels <- hull_facet(vapply(results, function(r) r$label, character(1)))

  dat <- do.call(rbind, lapply(results, function(r) {
    # `nfac` is only ever drawn as a text label (see .retention_count())
    data.frame(variant = hull_facet(r$label), df = r$x, fit = r$y,
               nfac = .retention_count(r$point_labels), on_hull = r$on_hull,
               stringsAsFactors = FALSE)
  }))
  dat$variant <- factor(dat$variant, levels = variant_levels)

  retained <- do.call(rbind, lapply(results, function(r) {
    idx <- which(r$point_labels == r$highlight)
    if (length(idx) == 0) return(NULL)
    data.frame(variant = hull_facet(r$label), df = r$x[idx], fit = r$y[idx],
               stringsAsFactors = FALSE)
  }))

  p <- ggplot2::ggplot(dat, ggplot2::aes(.data$df, .data$fit)) +
    ggplot2::geom_line(data = dat[dat$on_hull, , drop = FALSE]) +
    ggplot2::geom_point(ggplot2::aes(colour = .data$on_hull)) +
    ggplot2::geom_text(ggplot2::aes(label = .data$nfac), vjust = -1, size = 3) +
    ggplot2::scale_colour_manual(values = c(`TRUE` = "black", `FALSE` = "darkgray"),
                                 guide = "none")

  if (!is.null(retained)) {
    retained$variant <- factor(retained$variant, levels = variant_levels)
    p <- p + ggplot2::geom_point(data = retained, shape = 1, size = 4,
                                 colour = "red")
  }

  if (length(variant_levels) > 1) {
    p <- p + ggplot2::facet_wrap(ggplot2::vars(.data$variant), scales = "free_y")
  }

  p +
    ggplot2::labs(x = "Degrees of freedom", y = "Goodness of fit",
                  title = x$criterion[["label"]]) +
    .gg_theme()

}