R/print.efa_screen.R

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

#' Print and format an efa_screen object
#'
#' `print()` turns the factor-analysis screening diagnostics computed by
#' [efa_screen()] into a sectioned report with banded, colour-coded verdicts:
#' sampling adequacy and sphericity (the Kaiser-Meyer-Olkin measure and Bartlett's
#' test of sphericity), multicollinearity (the determinant and condition number of
#' the correlation matrix), the per-variable diagnostics, and, when raw data were
#' supplied, multivariate normality and multivariate outliers. It closes with a
#' consolidated list of actionable recommendations (for example, which items to
#' consider dropping, whether to prefer an ordinal or a robust estimator, and a
#' caveat that keeps an over-powered Bartlett's test from being over-trusted).
#' `format()` assembles the same report and returns it as a character vector;
#' `print()` is `cat(format(x), sep = "\n")`. The lines follow the active console
#' theme, so they are plain when colours are disabled (for example when captured
#' into a file or stripped with [cli::ansi_strip()]). `print()` does not draw a
#' plot.
#'
#' @param x An object of class `efa_screen` (output from [efa_screen()]).
#' @param digits Integer. The number of decimal places the reported values are
#'   rounded to. Default is 3.
#' @param ... Not used; for consistency with the generic.
#'
#' @returns `print()` returns its argument `x` invisibly. `format()` returns a
#'   character vector with the report lines.
#'
#' @family factor analysis suitability
#'
#' @export
#'
#' @method print efa_screen
#'
#' @examples
#' # From raw data
#' efa_screen(iris[, 1:4])
#'
#' # From a correlation matrix (supply N for Bartlett's test of sphericity)
#' efa_screen(test_models$baseline$cormat, N = 500)
#'
#' # format() returns the same lines as a character vector:
#' writeLines(format(efa_screen(test_models$baseline$cormat, N = 500)))
#'
print.efa_screen <- function(x, digits = 3, ...) {
  cat(format(x, digits = digits, ...), sep = "\n")
  invisible(x)
}

#' @rdname print.efa_screen
#' @export
#' @method format efa_screen
format.efa_screen <- function(x, digits = 3, ...) {

  raw <- !is.null(x$per_item)

  # Overall KMO band (reused Kaiser & Rice, 1974 cut-offs) and the significance flags
  # that several sections and the recommendations share.
  kmo <- x$kmo$KMO
  bart_sig <- .screen_is_sig(x$bartlett$p_value)
  # The multivariate-normality verdict is a count over the tests that are available. The
  # Henze-Zirkler p-value is withheld when its null approximation degenerates at many
  # variables, so the set is not always all three. `mvn_nonnormal` gates the
  # recommendations as before; the counts let the verdict sentence separate one marginal
  # rejection from a unanimous one.
  mvn_p <- if (raw && !inherits(x$normality, "efa_screen_no_mvn")) {
    c(x$normality$mardia$skewness_p, x$normality$mardia$kurtosis_p,
      x$normality$hz$p_value)
  }
  mvn_n_avail <- sum(!is.na(mvn_p))
  mvn_n_rej <- sum(vapply(mvn_p, .screen_is_sig, logical(1)))
  mvn_nonnormal <- mvn_n_rej > 0L

  cli::cli_format_method({

    # -- Sampling adequacy and sphericity --------------------------------------
    .print_efa_rule("Sampling adequacy and sphericity")

    if (!is.null(kmo) && !is.na(kmo)) {
      band <- .kmo_band(kmo)
      kval <- round(kmo, digits)
      band$alert("The overall KMO value for your data is {band$label} (Overall KMO = {kval}).",
                 wrap = TRUE)
      # The bands are named on screen for the same reason the condition-index rule is:
      # a verbal band is a convention, not a finding, and an unattributed cut-off
      # standing next to an attributed one reads as the latter.
      cli::cli_text(paste("These data are {band$suitability} suitable for factor",
                          "analysis (verbal bands: Kaiser & Rice, 1974)."))
    } else {
      cli::cli_alert_warning("The overall KMO value for your data is not available.",
                             wrap = TRUE)
    }

    cli::cli_text("")

    bart <- x$bartlett
    if (is.null(bart)) {
      cli::cli_alert_warning(
        "Bartlett's test of sphericity was not computed; no sample size (N) was supplied.",
        wrap = TRUE)
    } else if (!is.null(bart$p_value) && !is.na(bart$p_value)) {
      # Bartlett significance wording reused verbatim from format.efa_bartlett.
      if (bart$p_value < .05) {
        cli::cli_alert_success(
          "The Bartlett's test of sphericity was significant at an alpha level of .05.",
          wrap = TRUE)
        cli::cli_text("These data are probably suitable for factor analysis.")
      } else {
        cli::cli_alert_danger(
          "The Bartlett's test of sphericity was not significant at an alpha level of .05.",
          wrap = TRUE)
        cli::cli_text("These data are probably not suitable for factor analysis.")
      }
      cli::cli_verbatim(paste0("\U03C7\U00B2(", bart$df, ") = ", round(bart$chisq, 2),
                               ", ", cli::style_italic("p"), .screen_p_str(bart$p_value)))
    } else {
      .print_bartlett_no_result()
    }

    # -- Multicollinearity -----------------------------------------------------
    .print_efa_rule("Multicollinearity")

    # The determinant is the product of the p eigenvalues of R, so it shrinks
    # geometrically as variables are added even when every eigenvalue stays far from
    # zero. A fixed cut-off on it (below 0.00001; Field, 2018) is therefore a statement
    # about the number of variables as much as about the data: it fires on wide but
    # well-conditioned item pools and contradicts the condition index printed directly
    # under it. The determinant is reported as a number only, and the condition index
    # (above 30; Belsley, 1991) carries this section's verdict and the consolidated
    # recommendation, so the one cut-off lives in one place.
    det_R <- x$determinant
    dstr <- format(det_R, digits = digits, scientific = FALSE)
    cli::cli_alert_info(
      "Determinant: {dstr}. It falls as variables are added, so the condition index below carries the verdict.",
      wrap = TRUE)

    cond <- x$condition
    cstr <- formatC(cond, format = "f", digits = digits)
    # The verdict is on the condition *index*, and its bands live in .ci_band(). The
    # index is shown alongside the rule that governs it: quoting the rule without the
    # quantity invites the value to be read against the wrong scale. Above 30 the
    # strength is named as well, because a flat label would give an index of 35 and an
    # index of 3000 the same verdict.
    ci <- sqrt(cond)
    cistr <- formatC(ci, format = "f", digits = digits)
    ci_band <- .ci_band(ci)
    if (ci_band$flag) {
      ci_band$alert(
        "Condition number: {cstr} (condition index {cistr}). The index is above 30, which indicates a near linear dependency; its relative strength is {ci_band$strength} (Belsley, 1991).",
        wrap = TRUE)
    } else {
      ci_band$alert("Condition number: {cstr} (condition index {cistr}). {ci_band$note}",
                    wrap = TRUE)
    }

    # -- Per-variable diagnostics ----------------------------------------------
    .print_efa_rule("Per-variable diagnostics")
    # The standard display names are MSA (the per-variable measure of sampling adequacy,
    # stored as kmo_i) and SMC (the squared multiple correlation). The per-variable
    # values are laid out as base R's data-frame print, emitted verbatim so cli does not
    # reflow the aligned columns (as format.efa_kmo does for the KMO vector). This is the
    # one table not routed through .efa_num(), so it is also the one place a leading zero
    # is kept: the frame mixes bounded coefficients with variances and missing-data
    # percentages on the data's own scale, for which the leading zero is not redundant.
    disp <- .screen_per_item_display(x, digits)
    cli::cli_verbatim(utils::capture.output(print(disp)))

    # -- Multivariate normality (raw data only) --------------------------------
    if (raw) {
      .print_efa_rule("Multivariate normality")
      nm <- x$normality
      if (inherits(nm, "efa_screen_no_mvn")) {
        cli::cli_alert_warning(
          "Multivariate normality tests were skipped (singular complete-case covariance).",
          wrap = TRUE)
      } else {
        md <- nm$mardia
        .screen_mvn_line(.screen_is_sig(md$skewness_p),
                         paste0("Mardia's skewness: \U03C7\U00B2(", md$skewness_df, ") = ",
                                round(md$skewness, 2), ", p", .screen_p_str(md$skewness_p), "."))
        # An undefined kurtosis statistic is neither a pass nor a rejection, so it gets
        # an info alert, as a withheld Henze-Zirkler p-value does below.
        if (is.na(md$kurtosis)) {
          cli::cli_alert_info(
            "Mardia's kurtosis: not available at this number of complete cases.")
        } else {
          .screen_mvn_line(.screen_is_sig(md$kurtosis_p),
                           paste0("Mardia's kurtosis: z = ", round(md$kurtosis, 2),
                                  ", p", .screen_p_str(md$kurtosis_p), "."))
        }
        # A withheld Henze-Zirkler p-value is neither a pass nor a rejection, so it gets
        # an info alert rather than the success/danger pair the other two lines use.
        if (inherits(nm$hz, "efa_screen_no_hz")) {
          cli::cli_alert_info(
            paste0("Henze-Zirkler: HZ = ", round(nm$hz$statistic, 2),
                   ", p not available at this number of variables."))
        } else {
          .screen_mvn_line(.screen_is_sig(nm$hz$p_value),
                           paste0("Henze-Zirkler: HZ = ", round(nm$hz$statistic, 2),
                                  ", p", .screen_p_str(nm$hz$p_value), "."))
        }
        cli::cli_text(if (isTRUE(mvn_nonnormal)) {
          paste0("These data depart from multivariate normality: ", mvn_n_rej, " of the ",
                 mvn_n_avail, " tests ", .screen_plural(mvn_n_rej, "rejects", "reject"),
                 " it.")
        } else {
          paste0("These data are consistent with multivariate normality: none of the ",
                 mvn_n_avail, " tests rejects it.")
        })
        # The tests use the complete cases, which under non-ignorable missingness are not
        # a random subsample: with missing data the verdict above is a statement about a
        # subset, so say which one. Silent when every row is complete.
        if (.screen_incomplete(nm$n_complete, x$settings$n_obs)) {
          cli::cli_text("(Computed from {nm$n_complete} complete cases of the
                         {x$settings$n_obs} rows supplied.)")
        }
      }
    }

    # -- Outliers (raw data only) ----------------------------------------------
    if (raw) {
      .print_efa_rule("Outliers")
      o <- x$outliers
      if (inherits(o, "efa_screen_no_outliers")) {
        cli::cli_alert_warning("Outlier diagnostics were skipped.")
        # The recorded reason, as in the classical-fallback branch below: which of the
        # three ways the complete-case covariance failed decides what the user should do,
        # and "singular" alone points at collinearity even when the cause is missingness.
        if (!is.null(o$reason)) cli::cli_text(o$reason)
      } else {
        if (identical(o$method, "classical")) {
          cli::cli_alert_warning(paste(
            "A robust (MCD) covariance could not be computed; classical Mahalanobis",
            "distances were used."), wrap = TRUE)
          # The recorded reason (too few complete cases, near-collinear variables, or an
          # exact fit) and the consequence of the fallback: the classical covariance is
          # computed from every observation, outliers included, so the diagnostic no longer
          # has the high-breakdown property the MCD gives it.
          if (!is.null(o$fallback_reason)) cli::cli_text(o$fallback_reason)
          cli::cli_text(paste(
            "These distances come from a covariance the outliers themselves inflate, so",
            "the diagnostic is no longer high-breakdown and tends to under-flag."))
        }
        nf <- length(o$flagged)
        cstr_o <- format(o$cutoff, digits = digits, scientific = FALSE)
        # Name the distance for the estimate actually used (robust MCD, or classical
        # Mahalanobis when the robust covariance could not be formed).
        dist_label <- if (identical(o$method, "classical")) "Mahalanobis distance" else "robust distance"
        cli::cli_alert_info(paste0(
          nf, " of ", .screen_case_base(o$n_complete, x$settings$n_obs), " ",
          .screen_plural(nf, "was", "were"), " flagged as ",
          .screen_plural(nf, "a multivariate outlier", "multivariate outliers"),
          " (", dist_label, " > ", cstr_o, ")."), wrap = TRUE)
      }
    }

    # -- Recommendations -------------------------------------------------------
    .print_efa_rule("Recommendations")
    recs <- .screen_recommendations(x, raw, kmo, bart_sig, mvn_nonnormal,
                                    ci_band$flag, digits)
    # cli treats each bullet as a glue/cli template, so double any braces coming from
    # user variable names to render them literally (as .efa_emit_bullets does).
    msg <- gsub("}", "}}", gsub("{", "{{", recs$message, fixed = TRUE), fixed = TRUE)
    cli::cli_bullets(stats::setNames(msg, recs$symbol))

  })
}

# Kaiser & Rice (1974) verbal bands for the overall KMO value, and the single place they
# are written down: the highest band the value clears gives its label, the suitability it
# implies, and the severity (success >= .7, warning >= .6, danger below) as an alert
# function, the matching colour, and the equivalent cli_bullets() symbol, so that callers
# reporting the same value as an alert, a coloured label, or a bullet cannot drift apart.
.kmo_band <- function(kmo) {
  bands <- list(
    list(min = .9,   label = "marvellous",   colour = cli::col_green,  alert = cli::cli_alert_success, symbol = "v", suitability = "probably"),
    list(min = .8,   label = "meritorious",  colour = cli::col_green,  alert = cli::cli_alert_success, symbol = "v", suitability = "probably"),
    list(min = .7,   label = "middling",     colour = cli::col_green,  alert = cli::cli_alert_success, symbol = "v", suitability = "probably"),
    list(min = .6,   label = "mediocre",     colour = cli::col_yellow, alert = cli::cli_alert_warning, symbol = "!", suitability = "probably"),
    list(min = .5,   label = "miserable",    colour = cli::col_red,    alert = cli::cli_alert_danger,  symbol = "x", suitability = "hardly"),
    list(min = -Inf, label = "unacceptable", colour = cli::col_red,    alert = cli::cli_alert_danger,  symbol = "x", suitability = "not")
  )
  Find(function(b) kmo >= b$min, bands)
}

# Belsley's (1991) reading of the condition index, and the single place it is written
# down: the highest band the index clears gives the severity as an alert function, the
# verdict text, and `flag`, which gates both this section's verdict and the consolidated
# recommendation so the two cannot drift apart. The bands are steps on the progression
# 1, 3, 10, 30, 100, 300, 1000 that Belsley reads relative strength from (p. 38); they
# are not absolute cut-offs. An index of 10 or less is rarely of interest (pp. 43 and
# 46). Belsley gives 30 as one example value for flagging a near dependency, and a band
# must exceed it, not only reach it (p. 38). He calls 30 to 100 moderate (p. 42), and in
# his worked example he grades indexes of 35, 153, and 455 as moderately strong, strong,
# and very strong (pp. 44-45).
.ci_band <- function(ci) {
  bands <- list(
    list(above = 300,  alert = cli::cli_alert_danger,  flag = TRUE, strength = "very strong"),
    list(above = 100,  alert = cli::cli_alert_danger,  flag = TRUE, strength = "strong"),
    list(above = 30,   alert = cli::cli_alert_danger,  flag = TRUE, strength = "moderate"),
    list(above = 10,   alert = cli::cli_alert_warning, flag = FALSE,
         note = "The index is above 10, but not above 30, the value that flags a near linear dependency (Belsley, 1991)."),
    list(above = -Inf, alert = cli::cli_alert_success, flag = FALSE,
         note = "An index of 10 or less is rarely of interest (Belsley, 1991).")
  )
  # An index that is not available takes the lowest band, as the report has always done.
  if (is.na(ci)) return(bands[[length(bands)]])
  Find(function(b) ci > b$above, bands)
}

# TRUE when the complete cases the normality and outlier diagnostics use are fewer than
# the rows supplied, i.e. when the two possible denominators differ.
.screen_incomplete <- function(n_complete, n_obs) {
  !is.null(n_obs) && !is.na(n_obs) && !is.null(n_complete) && n_obs > n_complete
}

# The denominator of the complete-case diagnostics, named as such -- together with the
# number of rows supplied -- whenever missing values make the two differ, so a count or a
# rate cannot be read against the wrong base.
.screen_case_base <- function(n_complete, n_obs) {
  if (.screen_incomplete(n_complete, n_obs)) {
    paste0(n_complete, " complete cases of the ", n_obs, " rows supplied")
  } else {
    paste0(n_complete, " observations")
  }
}

# TRUE when a p-value is present and below .05 (guards NA / NULL p-values).
.screen_is_sig <- function(p) !is.null(p) && !is.na(p) && p < .05

# p-value tail formatting shared by the Bartlett and MVN lines and by format.efa_bartlett.
.screen_p_str <- function(p) {
  if (is.null(p) || is.na(p)) " = NA" else if (p < .001) " < .001" else paste0(" = ", round(p, 3))
}

# Singular/plural picker for a scalar count (avoids cli's numeric pluralizer).
.screen_plural <- function(n, singular, plural) if (n == 1L) singular else plural

# One multivariate-normality test line: a danger alert when significant, a success alert
# otherwise. `text` already carries the test statistic and p-value (no cli braces). Left
# unwrapped, unlike the prose alerts around it: it is a fixed-format statistic line, and
# wrapping would let a break fall inside the statistic.
.screen_mvn_line <- function(sig, text) {
  if (sig) cli::cli_alert_danger(text) else cli::cli_alert_success(text)
}

# Oxford-comma list of names ("a", "a and b", "a, b, and c").
.screen_and_list <- function(x) {
  if (length(x) == 1L) return(x)
  if (length(x) == 2L) return(paste(x, collapse = " and "))
  paste0(paste(x[-length(x)], collapse = ", "), ", and ", x[length(x)])
}

# Build the per-variable display data frame: the full per-item table for raw data
# (columns renamed to the standard MSA/SMC abbreviations), or MSA and SMC alone for a
# correlation-matrix input. Numeric columns are rounded to `digits`.
.screen_per_item_display <- function(x, digits) {
  if (!is.null(x$per_item)) {
    disp <- x$per_item
    num <- vapply(disp, is.numeric, logical(1))
    disp[num] <- lapply(disp[num], round, digits)
    names(disp)[names(disp) == "smc"] <- "SMC"
    names(disp)[names(disp) == "kmo_i"] <- "MSA"
    # `missing` is a percentage; the bare header reads as a count next to `variance`,
    # which is on the data's own scale.
    names(disp)[names(disp) == "missing"] <- "missing%"
    # A `flags` entry is NA for a variable treated as continuous, i.e. "no category
    # screening applies". Printed as <NA> that reads as a failed computation, so drop
    # the column when it says nothing at all and render the remaining NAs as a dash.
    if (all(is.na(disp$flags))) disp$flags <- NULL
    else disp$flags[is.na(disp$flags)] <- "-"
    disp
  } else {
    # make.unique() guards against duplicate variable names (a correlation matrix with
    # duplicate dimnames is a valid input), which base data.frame() row names reject; it
    # matches how the raw per-item table is keyed (efa_screen()).
    data.frame(MSA = round(x$kmo$KMO_i, digits),
               SMC = round(x$smc, digits),
               row.names = make.unique(names(x$kmo$KMO_i)))
  }
}

# Assemble the consolidated recommendations as parallel symbol/message vectors. Each
# check fires only when it applies; when no warning fires a single success line is
# returned, and a correlation-matrix input always closes with a raw-data note. Every
# recommendation is anchored in the literature cited in efa_screen()'s @source.
.screen_recommendations <- function(x, raw, kmo, bart_sig, mvn_nonnormal, multicollinear,
                                    digits) {
  sym <- character(0)
  msg <- character(0)
  push <- function(s, m) { sym[[length(sym) + 1L]] <<- s; msg[[length(msg) + 1L]] <<- m }

  # Few-category items (fewer than 5 response categories) are detected up front because
  # they steer the estimator recommendation away from the continuous-non-normal advice.
  few <- character(0)
  if (raw && !is.null(x$categories)) {
    n_cat <- vapply(x$categories, function(ct) {
      if (length(ct) == 1L && is.na(ct[1])) NA_integer_ else length(ct)
    }, integer(1))
    few <- names(n_cat)[!is.na(n_cat) & n_cat < 5L]
  }

  # Low overall sampling adequacy (KMO < .6): the set may not be factorable.
  if (!is.null(kmo) && !is.na(kmo) && kmo < .6) {
    push("!", paste0("Overall sampling adequacy is ", .kmo_band(kmo)$label,
                     " (KMO = ", round(kmo, digits),
                     "); the variables may not form a factorable set."))
  }

  # Variables with a low individual MSA (< .5) share little variance and are drop
  # candidates (Kaiser & Rice, 1974). make.unique() keys the names as the per-item table
  # and the sparse/empty recommendations do, so a variable is labelled consistently even
  # when the input has duplicate variable names.
  msa <- x$kmo$KMO_i
  low <- make.unique(names(msa))[!is.na(msa) & msa < .5]
  if (length(low)) {
    push("!", paste0(length(low), " ", .screen_plural(length(low), "variable has", "variables have"),
                     " a low individual MSA (< .5): ", .screen_and_list(low),
                     "; consider removing ", .screen_plural(length(low), "it", "them"),
                     " (little shared variance)."))
  }

  # Few-category ordinal data: prefer a categorical estimator (Rhemtulla et al., 2012).
  if (length(few)) {
    push("!", paste0(length(few), " ", .screen_plural(length(few), "item has", "items have"),
                     " fewer than 5 response categories; treating ",
                     .screen_plural(length(few), "it", "them"),
                     " as ordinal (polychoric correlations with a categorical estimator such as DWLS)",
                     " is less biased than normal-theory ML."))
  }

  # Continuous / 5+-category non-normal data: normal-theory standard errors are biased.
  # Suppressed when few-category items already routed the user to an ordinal estimator.
  if (isTRUE(mvn_nonnormal) && !length(few)) {
    push("!", paste("These data depart from multivariate normality; normal-theory standard errors",
                    "and fit statistics may be biased - prefer robust (sandwich) or bootstrapped",
                    "standard errors."))
  }

  # An over-powered Bartlett's test: significant but built on the normality it violates.
  if (bart_sig && isTRUE(mvn_nonnormal)) {
    push("!", paste("Bartlett's test is significant, but it assumes multivariate normality and grows",
                    "more sensitive as N increases; because these data are non-normal, treat it as",
                    "uninformative here and rely on the KMO."))
  }

  # Multicollinearity: a condition index above 30, the value Belsley (1991) gives for
  # flagging a near linear dependency, decided once by the caller from the same band
  # table as the section verdict. A high index does not have to come from one redundant
  # pair - a set of items that is together nearly linearly dependent gives the same
  # result - so the remedy names both.
  if (isTRUE(multicollinear)) {
    push("!", paste("A high condition index indicates multicollinearity; look for a redundant item",
                    "pair (r > .8) or a larger set of items that is nearly linearly dependent, and",
                    "consider removing some of them."))
  }

  if (raw) {
    flags <- x$per_item$flags
    items <- rownames(x$per_item)
    sparse <- items[!is.na(flags) & grepl("sparse", flags)]
    if (length(sparse)) {
      push("!", paste0(length(sparse), " ", .screen_plural(length(sparse), "variable has", "variables have"),
                       " a sparse response category (< 5 responses): ", .screen_and_list(sparse),
                       "; a low-frequency category can destabilise polychoric estimates",
                       " - consider collapsing it into an adjacent category."))
    }
    empty <- items[!is.na(flags) & grepl("empty", flags)]
    if (length(empty)) {
      push("!", paste0(length(empty), " ", .screen_plural(length(empty), "variable has", "variables have"),
                       " an unused interior response category (a gap in the scale): ",
                       .screen_and_list(empty), "; check the coding."))
    }
    if (!inherits(x$outliers, "efa_screen_no_outliers")) {
      nf <- length(x$outliers$flagged)
      nominal <- 1 - (x$settings$outlier_cutoff %||% 0.975)
      frac <- nf / x$outliers$n_complete
      # Far more flags than the cutoff nominally admits is not a contamination count to
      # work through: a high-breakdown estimate fitted to the most concentrated half of a
      # clustered or mixture sample legitimately calls the rest distant, so the excess is
      # evidence that the data are not elliptically distributed. The rate alone is not
      # enough to draw that conclusion - at a tight cutoff, or in a small sample, a couple
      # of genuine outliers clear it - so the flagged set must also be too long to work
      # through case by case, which is what makes the ordinary advice below unhelpful.
      if (nf >= 10L && frac > 5 * nominal) {
        # The rate is over the complete cases, so with missing data both bases are named:
        # read against the rows supplied it would look milder than it is.
        of_those <- if (.screen_incomplete(x$outliers$n_complete, x$settings$n_obs)) {
          "% of those cases"
        } else {
          "%"
        }
        push("!", paste0(nf, " of ",
                         .screen_case_base(x$outliers$n_complete, x$settings$n_obs),
                         " (", round(100 * frac), of_those,
                         ") exceed the outlier cutoff, far above the ",
                         # keep a significant digit: a tight cutoff has a nominal rate well
                         # below 0.1%, which rounds to a self-defeating "0%"
                         format(signif(100 * nominal, 2), trim = TRUE, scientific = FALSE),
                         "% expected under multivariate normality; this usually means the",
                         " data are not elliptically distributed (subgroups or a mixture)",
                         " rather than that this many cases are contaminated."))
      } else if (nf > 0L) {
        push("!", paste0(nf, " ", .screen_plural(nf, "observation was", "observations were"),
                         " flagged as ",
                         .screen_plural(nf, "a potential multivariate outlier", "potential multivariate outliers"),
                         "; inspect ", .screen_plural(nf, "it", "them"),
                         " (see `$outliers$flagged`) before down-weighting or excluding."))
      }
    }
  }

  # A clean bill of health when nothing above fired; the raw-data note always closes a
  # correlation-matrix report (its raw-only sections were absent above).
  if (!length(msg)) {
    push("v", "The data appear suitable for factor analysis.")
  }
  if (!raw) {
    push("i", paste("Per-item variance, missing-data, category, normality, and outlier diagnostics",
                    "require raw data; only a correlation matrix was supplied."))
  }

  list(symbol = sym, message = msg)
}