Exploratory Factor Analysis — Latent Structure Finder

suppressPackageStartupMessages(library(htmltools))
suppressPackageStartupMessages(library(jsonlite))
suppressPackageStartupMessages(library(plotly))
suppressPackageStartupMessages(library(DT))
suppressPackageStartupMessages(library(htmlwidgets))
suppressPackageStartupMessages(library(arrow))
suppressPackageStartupMessages(library(knitr))
suppressPackageStartupMessages(library(rmarkdown))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(tidyr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(stringr))
suppressPackageStartupMessages(library(lubridate))
suppressPackageStartupMessages(library(broom))
suppressPackageStartupMessages(library(Matrix))
suppressPackageStartupMessages(library(cluster))
suppressPackageStartupMessages(library(data.table))

initial_rows <- nrow(df)
  feat_cols <- grep("^feature_[0-9]+$", names(df), value = TRUE)
  feat_cols <- feat_cols[order(as.integer(sub("^feature_", "", feat_cols)))]
  if (length(feat_cols) < 3) {
    stop("column_mapping must map at least three feature columns(feature_1, feature_2, feature_3) — factor analysis needs 3 or more variables")
  }
  feature_names <- setNames(humanize_semantic(feat_cols, col_map), feat_cols)

dropped_features <- character(0)
  for (fc in feat_cols) {
    v <- df[[fc]]
    if (!is.numeric(v)) {
      conv <- suppressWarnings(as.numeric(as.character(v)))
      n_orig <- sum(!is.na(v) & as.character(v) != "")
      if (n_orig > 0 && sum(!is.na(conv)) >= 0.95 * n_orig) {
        df[[fc]] <- conv
      } else {
        dropped_features <- c(dropped_features, fc)
        next
      }
    }
    v <- df[[fc]]
    med <- median(v, na.rm = TRUE)
    if (is.na(med)) { dropped_features <- c(dropped_features, fc); next }
    v[is.na(v)] <- med
    df[[fc]] <- v
    if (isTRUE(var(v) == 0) || is.na(var(v))) {
      dropped_features <- c(dropped_features, fc)
    }
  }
  used_features <- setdiff(feat_cols, dropped_features)
  if (length(used_features) < 3) {
    stop(paste0(
      "Factor analysis needs at least three usable numeric columns; only ",
      length(used_features), " remained after cleaning. Excluded: ",
      paste(feature_names[dropped_features], collapse = ", ")
    ))
  }

  X <- as.matrix(df[, used_features, drop = FALSE])
  final_rows <- nrow(X)
  rows_removed <- initial_rows - final_rows
  if (final_rows < 20) {
    stop(sprintf("Only %d usable rows — factor analysis needs at least 20.", final_rows))
  }
  p <- length(used_features)
  n <- final_rows

Xs <- scale(X)

k_try_max <- max(1L, min(5L, as.integer(floor((p - sqrt(p)) / 2))))
  fits <- vector("list", k_try_max)
  for (k in seq_len(k_try_max)) {
    fits[[k]] <- tryCatch(
      suppressWarnings(factanal(x = Xs, factors = k, rotation = "varimax",
                                scores = "none")),
      error = function(e) NULL
    )
    if (!is.null(fits[[k]]) && isFALSE(fits[[k]]$converged)) fits[[k]] <- NULL
  }
  ok_ks <- which(!vapply(fits, is.null, logical(1)))
  if (length(ok_ks) == 0) {
    stop(paste0(
      "Factor analysis could not be fit to these columns(",
      paste(feature_names[used_features], collapse = ", "),
      ") — they may be near-duplicates of each other or too few rows for ",
      "the model. Try removing redundant columns or adding rows."
    ))
  }
  retried <- length(ok_ks) < k_try_max

  get_pval <- function(f) {
    pv <- suppressWarnings(as.numeric(f$PVAL %||% NA_real_))
    if (length(pv) == 0) NA_real_ else pv[1]
  }
  get_bic <- function(f) {
    st <- suppressWarnings(as.numeric(f$STATISTIC %||% NA_real_))
    dof <- suppressWarnings(as.numeric(f$dof %||% NA_real_))
    if (length(st) == 0 || length(dof) == 0 || is.na(st) || is.na(dof) ||
        dof <= 0) return(NA_real_)
    st - dof * log(n)
  }
  pvals <- vapply(fits[ok_ks], get_pval, numeric(1))
  bics  <- vapply(fits[ok_ks], get_bic, numeric(1))

  ev <- eigen(cor(Xs), symmetric = TRUE, only.values = TRUE)$values
  kaiser_k <- max(1L, sum(ev > 1))

  fit_ok <- !is.na(pvals) & pvals > 0.05
  if (any(fit_ok)) {
    n_factors <- ok_ks[fit_ok][1]                      # smallest adequate k
  } else if (any(!is.na(bics))) {
    valid <- which(!is.na(bics))                       # filter NA BEFORE which.min
    n_factors <- ok_ks[valid][which.min(bics[valid])]
  } else {
    below <- ok_ks[ok_ks <= kaiser_k]                  # Kaiser fallback, capped
    n_factors <- if (length(below) > 0) max(below) else min(ok_ks)
  }

is_heywood <- function(f) any(f$uniquenesses <= 0.005 + 1e-6)
  heywood <- FALSE
  while (is_heywood(fits[[n_factors]]) && any(ok_ks < n_factors)) {
    retried <- TRUE
    n_factors <- max(ok_ks[ok_ks < n_factors])
  }
  if (is_heywood(fits[[n_factors]])) heywood <- TRUE   # k = 1 minimum: accept + flag
  fa <- fits[[n_factors]]
  fit_pval <- get_pval(fa)

L <- unclass(fa$loadings)[, seq_len(n_factors), drop = FALSE]
  ss <- colSums(L^2)
  ord <- order(-ss)                                    # factors by variance, desc
  L <- L[, ord, drop = FALSE]
  ss <- ss[ord]
  var_pct <- 100 * ss / p
  cumulative_pct <- sum(var_pct)
  hum <- unname(feature_names[used_features])

  factor_names_short <- paste0("Factor ", seq_len(n_factors))
  top2_join <- character(n_factors)
  top_vars_list <- character(n_factors)
  loading_groups <- list()
  for (j in seq_len(n_factors)) {
    aload <- abs(L[, j])
    idx <- order(-aload)
    top2_join[j] <- paste(hum[idx[seq_len(min(2, p))]], collapse = " + ")
    strong <- which(aload >= 0.3)
    strong <- strong[order(-aload[strong])]
    if (length(strong) == 0) strong <- idx[1]
    top_vars_list[j] <- paste(hum[head(strong, 5)], collapse = ", ")
    loading_groups[[factor_names_short[j]]] <- hum[strong]
  }
  factor_names_full <- paste0(factor_names_short, ": ", top2_join)

  variance_by_factor_df <- data.frame(
    factor_name  = factor_names_full,
    variance_pct = round(var_pct, 1),
    stringsAsFactors = FALSE
  )
  # var_pct carries Factor1/Factor2 names; reset so no "Row" column renders.
  rownames(variance_by_factor_df) <- NULL

  loadings_df <- data.frame(variable = hum, stringsAsFactors = FALSE)
  for (j in seq_len(n_factors)) {
    loadings_df[[paste0("f", j)]] <- round(L[, j], 3)
  }
  loadings_df$uniqueness <- round(unname(fa$uniquenesses[used_features]), 2)
  rownames(loadings_df) <- NULL

  interpretation_df <- data.frame(
    factor_name       = factor_names_short,
    top_variables     = top_vars_list,
    variance_pct      = round(var_pct, 1),
    suggested_meaning = top2_join,
    stringsAsFactors  = FALSE
  )
  # var_pct carries Factor1/Factor2 names; reset so no "Row" column renders.
  rownames(interpretation_df) <- NULL

  uniqueness_df <- data.frame(
    variable       = hum,
    uniqueness_pct = round(100 * unname(fa$uniquenesses[used_features]), 1),
    stringsAsFactors = FALSE
  )
  uniqueness_df <- uniqueness_df[order(-uniqueness_df$uniqueness_pct), , drop = FALSE]
  rownames(uniqueness_df) <- NULL

  metrics <- list(
    `Observations`               = final_rows,
    `Variables Analysed`         = p,
    `Factors Found`              = as.integer(n_factors),
    `Total Variance Explained %` = round(cumulative_pct, 1),
    `Fit Test P Value`           = fmt_p(fit_pval),
    `Top Factor`                 = factor_names_full[1]
  )

  json_output <- list(
    answer = paste0(
      "Maximum-likelihood factor analysis on ", p, " variables across ",
      format(final_rows, big.mark = ","), " rows found ", n_factors,
      " latent ", plural_word(n_factors, "factor"), ": ",
      paste0(factor_names_full, " (", round(var_pct, 1), "% of variance)",
             collapse = "; "),
      ". Together the shared ",
      plural_word(n_factors, "factor explains", "factors explain"), " ",
      round(cumulative_pct, 1), "% of total variation(likelihood-ratio fit ",
      "test p ", fmt_p(fit_pval), ")."
    ),
    cards = lapply(
      c("tldr", "overview", "preprocessing", "variance_by_factor",
        "loadings_table", "factor_names", "uniqueness_chart"),
      function(cid) list(id = cid, metrics = metrics)
    )
  )

  list(
    initial_rows          = initial_rows,
    final_rows            = final_rows,
    rows_removed          = rows_removed,
    feature_names         = feature_names,
    used_features         = used_features,
    dropped_features      = dropped_features,
    n_factors             = as.integer(n_factors),
    k_try_max             = k_try_max,
    fit_pval              = fit_pval,
    heywood               = heywood,
    retried               = retried,
    factor_names_short    = factor_names_short,
    factor_names_full     = factor_names_full,
    variance_by_factor_df = variance_by_factor_df,
    var_pct               = var_pct,
    cumulative_pct        = cumulative_pct,
    loadings_df           = loadings_df,
    interpretation_df     = interpretation_df,
    uniqueness_df         = uniqueness_df,
    loading_groups        = loading_groups,
    kaiser_k              = kaiser_k,
    metrics               = metrics,
    json_output           = json_output
  )
}

shared <- compute_shared(df, params, col_map)