Factor Analysis

The code that produced every figure in this report, exactly as it ran. Fingerprint 08224f67fd28947c. The same code on the same data gives the same report.

`standard_factor_analysis_v2` <- function(pf) {
  `%||%` <- function(a, b) if (!is.null(a)) a else b
  #' Readable figures (LAT-3180, LAT-3181): whole numbers from a thousand up, one decimal from a hundred, two from
  #' one, three significant figures below one. A cell carries what the value needs, not what R prints.
  tidy <- function(x) {
    x <- as.numeric(x)
    ifelse(is.na(x), NA_real_,
      ifelse(abs(x) >= 1000, round(x, 0),
        ifelse(abs(x) >= 100, round(x, 1),
          ifelse(abs(x) >= 1, round(x, 2), signif(x, 3)))))
  }
  #' A p-value under 1e-4 must LEAVE the tool as a small number, not as 0: results serialise at four decimal digits,
  #' so 1.4e-05 would arrive as 0 and a table would print "p-value 0" (LAT-3181).
  p_cell  <- function(p) if (is.na(p)) NA_real_ else max(p, 1e-12)
  p_text  <- function(p) if (is.na(p)) "" else if (p < 1e-4) "<0.0001" else format(signif(p, 3))
  pct     <- function(x) tidy(100 * x)

  inputs <- pf$taskList$inputs
  params <- inputs$module_parameters %||% list()
  # THE QUESTION this tool answers: the customer's objective, verbatim, when given.
  question <- (inputs$userContext %||% list())$objective %||%
    "What shared factors sit behind these measures, and how much of each measure do they explain?"

  #' ## Column mapping
  #' Three or more numeric `measure_N` columns (a series, any number) describing the same rows. Semantic names
  #' inside; the customer's headers live in `col_map`.
  col_map <- inputs$column_mapping %||% list()
  df <- renderObject.taskFunction.init(inputs, col_map)   # df has SEMANTIC names
  human <- function(sem) {
    v <- col_map[[sem]]
    if (is.null(v) || !nzchar(as.character(v))) sem else as.character(v)
  }
  n_in <- nrow(df)
  meas_cols <- grep("^measure_[0-9]+$", names(df), value = TRUE)
  meas_cols <- meas_cols[order(as.integer(sub("^measure_", "", meas_cols)))]
  if (length(meas_cols) < 3) stop("column_mapping must map at least three numeric columns to measure_1, measure_2, measure_3, ... (the measures whose shared factors you want).")

  raw_names <- local({
    ds <- inputs$dataset %||% inputs$df
    if (is.data.frame(ds)) return(names(ds))
    if (is.list(ds) && length(ds) > 0) {
      rows <- ds[seq_len(min(length(ds), 50))]
      nm <- unique(unlist(lapply(rows, function(r) if (is.list(r)) names(r) else NULL)))
      if (length(nm)) return(nm)
      if (!is.null(names(ds)) && all(nzchar(names(ds)))) return(names(ds))
    }
    character(0)
  })
  mapped_actual <- unique(as.character(unlist(col_map)))
  ignored_cols  <- setdiff(raw_names, unique(c(mapped_actual, make.names(mapped_actual))))

  #' ## Parameters
  #' `factors`: how many factors to fit. Unset (the default) fits the smallest number the likelihood-ratio test does
  #' not reject, which is the usual choice. A number fits exactly that many.
  factors_param <- params$factors %||% NULL
  if (!is.null(factors_param)) {
    factors_param <- suppressWarnings(as.integer(factors_param))
    if (is.na(factors_param) || factors_param < 1) stop("module_parameters$factors must be a whole number of 1 or more, or left unset.")
  }
  #' `rotation`: varimax (the default) turns the factors so each is driven by its own group of measures, which is what
  #' makes them readable; none leaves them as fitted. Rotation needs at least two factors, and the method says when it
  #' did not apply.
  rotation_param <- tolower(as.character(params$rotation %||% "varimax"))
  if (!(rotation_param %in% c("varimax", "none"))) stop("module_parameters$rotation must be varimax or none")

  #' ## Data preparation
  #' Each mapped measure is read as a number (95% rule). A measure that is not numeric, is constant, or is a running
  #' index is excluded and named. Rows missing any kept measure are excluded and counted (a factor is estimated from
  #' the correlations between every measure, so a filled value would move it). At least three measures and 30 complete
  #' rows are required: maximum likelihood needs more rows than a description does.
  excluded_cols <- character(0); why <- character(0); X <- list()
  for (mc in meas_cols) {
    v <- df[[mc]]; nm <- human(mc)
    if (!is.numeric(v)) {
      ch <- trimws(as.character(v)); nb <- !is.na(ch) & ch != ""
      conv <- suppressWarnings(as.numeric(ch))
      if (sum(nb) == 0 || sum(!is.na(conv[nb])) < 0.95 * sum(nb)) { excluded_cols <- c(excluded_cols, nm); why <- c(why, "not numeric"); next }
      v <- conv
    }
    v <- as.numeric(v); v[!is.finite(v)] <- NA; ok <- v[!is.na(v)]
    if (length(ok) < 3 || isTRUE(stats::sd(ok) == 0)) { excluded_cols <- c(excluded_cols, nm); why <- c(why, "constant"); next }
    if (length(unique(ok)) == length(ok) && all(abs(diff(sort(ok)) - 1) < 1e-9)) { excluded_cols <- c(excluded_cols, nm); why <- c(why, "a running index"); next }
    X[[nm]] <- v
  }
  if (length(X) < 3) stop(sprintf("Fewer than three usable numeric measures remained (excluded: %s); factor analysis needs at least three.",
                                  if (length(excluded_cols)) paste(paste0(excluded_cols, " (", why, ")"), collapse = ", ") else "none"))
  X <- as.data.frame(X, check.names = FALSE, stringsAsFactors = FALSE)
  complete <- stats::complete.cases(X)
  n_incomplete <- sum(!complete)
  X <- X[complete, , drop = FALSE]
  M <- as.matrix(X); n <- nrow(M); k_meas <- ncol(M)
  if (n < 30) stop(sprintf("Only %d rows have every measure; maximum-likelihood factor analysis needs at least 30.", n))

  R <- stats::cor(M)
  if (any(!is.finite(R))) stop("Two or more measures could not be correlated (a measure has no variation within the complete rows).")
  #' A singular correlation matrix means one measure is a combination of the others: factor analysis cannot separate
  #' them, and the tool says which rather than failing inside the fit.
  if (abs(det(R)) < 1e-12) {
    ev <- eigen(R, symmetric = TRUE, only.values = TRUE)$values
    stop(sprintf("These measures are collinear (the smallest eigenvalue of their correlation matrix is %s), so no factor model can separate them. Drop a measure that duplicates another and re-run.", signif(min(ev), 3)))
  }

  #' ## Sampling adequacy, computed in base R
  #' The runtime image carries no psych, so KMO and Bartlett are computed from the correlation matrix directly. Both
  #' forms were checked against psych::KMO and psych::cortest.bartlett on the same matrix and agree to 1e-12.
  #' KMO compares the correlations to the PARTIAL correlations (the anti-image): when the partials are small relative
  #' to the correlations, the measures share something a factor can carry.
  Rinv <- solve(R)
  dinv <- sqrt(diag(Rinv))
  Q <- -Rinv / outer(dinv, dinv); diag(Q) <- 1
  off <- !diag(k_meas)
  kmo_overall <- sum(R[off]^2) / (sum(R[off]^2) + sum(Q[off]^2))
  kmo_each <- vapply(seq_len(k_meas), function(i) {
    idx <- setdiff(seq_len(k_meas), i)
    r2 <- sum(R[i, idx]^2); q2 <- sum(Q[i, idx]^2)
    r2 / (r2 + q2)
  }, numeric(1))
  names(kmo_each) <- colnames(R)
  #' Bartlett's test of sphericity: is the correlation matrix distinguishable from the identity at all? A test that
  #' does NOT reject means the measures are unrelated and there is nothing for a factor to explain.
  bart_chi <- -((n - 1) - (2 * k_meas + 5) / 6) * log(det(R))
  bart_df  <- k_meas * (k_meas - 1) / 2
  bart_p   <- stats::pchisq(bart_chi, bart_df, lower.tail = FALSE)

  #' ## How many factors
  #' The most factors maximum likelihood can identify from k measures is bounded by the degrees of freedom; beyond it
  #' the model is not estimable and factanal refuses.
  max_factors <- max(1L, floor((2 * k_meas + 1 - sqrt(8 * k_meas + 1)) / 2))
  fit_k <- function(kf) tryCatch(stats::factanal(M, factors = kf, rotation = "none"), error = function(e) NULL)
  #' The smallest number the likelihood-ratio test does not reject at 0.05. This is a floor, NOT a selection: a larger
  #' number can fit just as well, so every sentence says k factors are ENOUGH and that k-1 is rejected.
  lr_table <- list()
  chosen <- NA_integer_
  for (kf in seq_len(max_factors)) {
    f <- fit_k(kf)
    if (is.null(f)) next
    pv <- if (is.null(f$PVAL) || is.na(f$PVAL)) NA_real_ else as.numeric(f$PVAL)
    lr_table[[length(lr_table) + 1]] <- data.frame(factors = kf, p_value = pv, dof = f$dof, stringsAsFactors = FALSE)
    if (is.na(chosen) && (is.na(pv) || pv >= 0.05)) chosen <- kf
  }
  if (!length(lr_table)) stop("No factor model could be fitted to these measures. With this many measures and rows, maximum likelihood has too few degrees of freedom.")
  lr_df <- do.call(rbind, lr_table)
  n_keep <- if (!is.null(factors_param)) {
    if (factors_param > max_factors) stop(sprintf("module_parameters$factors is %d, but %d measures support at most %d factors by maximum likelihood.", factors_param, k_meas, max_factors))
    as.integer(factors_param)
  } else if (!is.na(chosen)) as.integer(chosen) else as.integer(max_factors)

  rotated <- rotation_param == "varimax" && n_keep >= 2
  fa <- stats::factanal(M, factors = n_keep, rotation = if (rotated) "varimax" else "none")
  L <- unclass(fa$loadings)[, seq_len(n_keep), drop = FALSE]
  #' A factor's sign is arbitrary: fix it so the strongest-loading measure is positive, or a reader sees a reversal
  #' the data does not contain.
  for (j in seq_len(n_keep)) if (L[which.max(abs(L[, j])), j] < 0) L[, j] <- -L[, j]
  fac_names <- if (rotated) paste0("Factor ", seq_len(n_keep)) else paste0("Factor ", seq_len(n_keep))
  colnames(L) <- fac_names

  uniq  <- fa$uniquenesses[rownames(L)]
  commu <- 1 - uniq
  ss    <- colSums(L^2)
  #' Every figure a reader can add up must derive from the figures SHOWN: totalling the unrounded
  #' shares and rounding once printed a cumulative of 63.22 above parts reading 32.03 and 31.18,
  #' which sum to 63.21. The parts are rounded first and every total is built from them.
  var_pct <- round(100 * ss / k_meas, 2)
  cum_pct <- round(cumsum(var_pct), 2)
  var_total <- round(sum(var_pct), 2)

  #' Scree against a random baseline (Horn), on the CORRELATION matrix's eigenvalues, seed 42 so the line is stable.
  eig_obs <- eigen(R, symmetric = TRUE, only.values = TRUE)$values
  set.seed(42)
  rand <- replicate(100, {
    Rr <- stats::cor(matrix(stats::rnorm(n * k_meas), n, k_meas))
    eigen(Rr, symmetric = TRUE, only.values = TRUE)$values
  })
  rand95 <- apply(rand, 1, stats::quantile, probs = 0.95)

  #' ## Places
  lead <- which.max(var_pct)
  strongest_idx <- apply(abs(L), 1, which.max)
  measures_df <- data.frame(
    measure           = rownames(L),
    communality       = pct(commu),
    uniqueness        = pct(uniq),
    strongest_factor  = fac_names[strongest_idx],
    strongest_loading = tidy(L[cbind(seq_len(nrow(L)), strongest_idx)]),
    stringsAsFactors = FALSE)
  measures_df <- measures_df[order(-measures_df$communality), , drop = FALSE]

  factors_df <- data.frame(
    factor         = fac_names,
    ss_loadings    = tidy(ss),
    variance_pct   = var_pct,
    cumulative_pct = cum_pct,
    stringsAsFactors = FALSE)

  load_long <- data.frame(
    measure = rep(rownames(L), times = n_keep),
    factor  = rep(fac_names, each = nrow(L)),
    loading = tidy(as.vector(L)),
    stringsAsFactors = FALSE)

  scree_df <- data.frame(
    factor     = rep(seq_len(k_meas), 2),
    eigenvalue = tidy(c(eig_obs, rand95)),
    series     = rep(c("This data", "Random 95th"), each = k_meas),
    stringsAsFactors = FALSE)

  commu_df <- data.frame(measure = rownames(L), communality = pct(commu), stringsAsFactors = FALSE)
  commu_df <- commu_df[order(-commu_df$communality), , drop = FALSE]
  varexp_df <- data.frame(factor = fac_names, variance_pct = var_pct, stringsAsFactors = FALSE)

  #' ## Checks
  #' Every row states what it tested and what it found. Two things this tool must never say: that a test CHOSE the
  #' number of factors (a goodness-of-fit test that fails to reject has only failed to rule it out, and a larger
  #' number often fits equally well), and that sampling adequacy says a measure is well explained (KMO is about
  #' whether a measure belongs in the analysis, communality about how much of it the factors account for; the two
  #' routinely disagree).
  p_below <- lr_df$p_value[lr_df$factors == (n_keep - 1)]
  enough_p <- lr_df$p_value[lr_df$factors == n_keep]
  rows_per_measure <- n / k_meas
  checks_df <- data.frame(
    check = c(
      "Enough factors to fit the correlations",
      "The measures are correlated at all",
      "The measures belong in one analysis",
      "Rows behind each measure"),
    statistic = c(
      if (length(enough_p) && !is.na(enough_p[1]))
        sprintf("%d factors, fit p %s%s", n_keep, p_text(enough_p[1]),
                if (length(p_below) && !is.na(p_below[1])) sprintf("; %d rejected at p %s", n_keep - 1, p_text(p_below[1])) else "")
      else sprintf("%d factors; the fit test has no degrees of freedom left", n_keep),
      sprintf("chi-square %s on %d df", format(tidy(bart_chi)), bart_df),
      sprintf("KMO %s overall; lowest measure %s", format(tidy(kmo_overall)), format(tidy(min(kmo_each)))),
      sprintf("%s rows per measure (%d rows, %d measures)", format(tidy(rows_per_measure)), n, k_meas)),
    p_value = c(
      if (length(enough_p)) p_cell(enough_p[1]) else NA_real_,
      p_cell(bart_p), NA_real_, NA_real_),
    verdict = c(
      if (length(enough_p) && !is.na(enough_p[1]) && enough_p[1] >= 0.05) "holds" else if (length(enough_p) && !is.na(enough_p[1])) "strained" else "not testable",
      if (bart_p < 0.05) "holds" else "strained",
      if (kmo_overall >= 0.6) "holds" else "strained",
      if (rows_per_measure >= 10) "holds" else "strained"),
    note = c(
      sprintf("The model with %d factors reproduces the correlations; the test does not reject it. A larger number may fit as well, so this is the smallest number that is enough rather than a number the test selected.", n_keep),
      if (bart_p < 0.05) "The correlation matrix differs from one with no correlations at all, so there is shared variation for a factor to carry." else "The measures are close to uncorrelated, and a factor model has little to explain.",
      "Sampling adequacy asks whether a measure belongs in the analysis. It is not a statement about how much of that measure the factors explain, which is its communality and can be low while adequacy is high.",
      "Maximum likelihood is a large-sample method; ten rows per measure is the common floor."),
    stringsAsFactors = FALSE)

  method <- sprintf(
    "Maximum-likelihood exploratory factor analysis (stats::factanal) of %d measures on %d of %d rows%s%s. %s %d factors were fitted, and %s. Loadings are %s, and each factor's sign was set so its strongest measure loads positively. A measure's communality is the share of it the factors explain and its uniqueness the remainder. Sampling adequacy is Kaiser-Meyer-Olkin computed from the anti-image correlation matrix, and Bartlett's test of sphericity is a chi-square on the determinant of the correlation matrix; both are computed directly from the correlations. The p-values test model fit, not the size of any loading, and are not adjusted.",
    k_meas, n, n_in,
    if (n_incomplete > 0) sprintf(", %d rows excluded for a missing measure", n_incomplete) else "",
    if (length(excluded_cols)) sprintf(", excluding %s", paste(paste0(excluded_cols, " (", why, ")"), collapse = ", ")) else "",
    if (!is.null(factors_param)) "As requested," else "The number of factors is the smallest the likelihood-ratio test does not reject at 0.05, which is a floor rather than a selection, so a larger number may fit as well.",
    n_keep,
    if (length(enough_p) && !is.na(enough_p[1])) sprintf("the fit test for %d gives p %s", n_keep, p_text(enough_p[1])) else "the fit test has no degrees of freedom left",
    if (rotated) "varimax-rotated so each factor is driven by its own group of measures" else "unrotated")

  assumptions <- list(
    "A factor model explains only the variance the measures share; each measure's own variance is left out by design, so the factors are not expected to account for everything.",
    "A measure with a low communality is not explained by these factors. That is a statement about the model, not about the quality of the measure.",
    "Sampling adequacy and communality answer different questions and often disagree for the same measure.",
    "The number of factors is the smallest the fit test does not reject. A larger number may fit as well; the test rules numbers out, it does not choose one.",
    "Loadings are correlations between a measure and a factor, so they describe how things move together and not what causes what.")

  answer <- sprintf("%d factors account for %s%% of the variation across %d measures, led by %s at %s%%.",
                    n_keep, format(var_total), k_meas, fac_names[lead], format(var_pct[lead]))

  results <- list()
  results$summary_metrics <- place_metric(list(
      factors_kept = n_keep,
      variance_kept_pct = var_total,
      kmo = tidy(kmo_overall),
      measures = k_meas,
      rows = n),
    lead = "factors_kept",
    place = "summary_metrics")
  results$scree              <- place_trend(scree_df, x = "factor", y = "eigenvalue", series = "series", place = "scree")
  results$communalities      <- place_comparison(commu_df, category = "measure", value = "communality", place = "communalities")
  results$loadings_matrix    <- place_matrix(load_long, x = "factor", y = "measure", z = "loading", place = "loadings_matrix")
  results$variance_explained <- place_comparison(varexp_df, category = "factor", value = "variance_pct", place = "variance_explained")
  results$factors_table      <- place_table(factors_df, place = "factors_table")
  results$measures_table     <- place_table(measures_df, place = "measures_table")
  results$assumption_checks  <- place_table(checks_df, place = "assumption_checks")
  results$factor_method <- list(kind = "metric", values = list(
      method = method, n_in = n_in, n_used = n,
      excluded = as.list(unname(c(excluded_cols, ignored_cols))),
      assumptions = assumptions,
      x_column = paste(colnames(M), collapse = ", "),
      y_column = sprintf("%d factors", n_keep)),
    value_order = list("n_used", "n_in"))

  objects <- list()   # filled by the object layer, not here
  list(answer = answer, method = method, n = n, results = results, objects = objects,
       json_output = list(answer = answer, method = method, n = n))
}