R/efa_compare.R

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

#' Compare two vectors or matrices (communalities or loadings)
#'
#' The function takes two objects of the same dimensions containing numeric
#' information (loadings or communalities) and returns a list of class
#' `efa_compare` containing summary information of the differences of the objects.
#'
#' @details `digits`, `m_red`, `range_red`, `round_red`, `print_diff`, and `plot_red`
#'  only control how the result is displayed; each is stored in the returned object's
#'  `settings` and can be overridden later without recomputing the comparison --
#'  `digits`, `m_red`, `range_red`, `round_red`, and `print_diff` in a call to
#'  [print.efa_compare()], and `plot_red` in a call to [plot.efa_compare()].
#'
#' @param x matrix, or vector. Loadings or communalities of a factor
#'  analysis output.
#' @param y matrix, or vector. Loadings or communalities of another
#'  factor analysis output to compare to x.
#' @param reorder character. Whether and how elements / columns should be
#' reordered. If "congruence" (default), the columns of `y` are matched to those of
#' `x` by a joint one-to-one assignment that maximizes the total Tucker's congruence
#' coefficient (a standard measure of similarity between two loading vectors) across
#' all columns at once, and each matched column's sign is flipped if needed. This way,
#' mismatched factor order or sign between two solutions does not distort the
#' comparison. It applies to matrices only, and warns when `x`
#' and `y` are vectors. If "names", the columns of a matrix -- or the elements of
#' a vector -- are put in alphabetical order of their names; the rows of a matrix
#' are assumed to be aligned already and are left untouched. If "none", no
#' reordering is done.
#' @param corres logical. Whether factor correspondences should be compared if a
#'  matrix is entered. Default is TRUE.
#' @param thresh numeric. The threshold at or above which a loading is classified as substantial. Default is .3.
#' @param digits numeric. Number of decimals to print in the output. Default is 4.
#' @param m_red numeric. Number above which the mean and median should be printed
#'  in red (i.e., if .001 is used, the mean will be in red if it is larger than
#'  .001, otherwise it will be displayed in green.) Default is .001.
#' @param range_red numeric. Number above which the min and max should be printed
#'  in red (i.e., if .001 is used, min and max will be in red if the max is larger
#'   than .001, otherwise it will be displayed in green). Default is .001. Note that
#'   the color of min also depends on max, that is min will be displayed in the
#'   same color as max.
#' @param round_red numeric. The number of agreeing decimals below which the
#'  report highlights the agreement in red (i.e., if 3 is used, the value is
#'  shown in red when the compared numbers agree to fewer than 3 decimals,
#'  otherwise in green). Default is 3.
#' @param print_diff logical. Whether the difference vector or matrix should be
#'  printed or not. Default is TRUE.
#' @param na.rm logical. Whether NAs should be removed from the difference
#'  summaries and factor-correspondence classifications. With `FALSE`, a missing
#'  loading makes the correspondence counts undefined (`NA`). Default is FALSE.
#' @param x_labels character. A vector of length two containing identifying
#'  labels for the two objects x and y that will be compared. These will be used
#'  as labels on the x-axis of the plot, and to name the direction of the signed
#'  elementwise differences in the printed report (see [print.efa_compare()]).
#'  Default is "x" and "y".
#' @param plot `r lifecycle::badge("superseded")` Accepted and validated, but
#'  without effect; retained for backwards compatibility. The difference plot is
#'  drawn by [plot.efa_compare()]. Default is TRUE.
#' @param plot_red numeric. Threshold above which to plot the absolute differences
#'  in red. Default is .01.
#'
#' @return A list of class `efa_compare` with the following components:
#'
#' \item{diff}{The vector or matrix containing the differences between x and y.}
#' \item{mean_abs_diff}{The mean absolute difference between x and y.}
#' \item{median_abs_diff}{The median absolute difference between x and y.}
#' \item{min_abs_diff}{The minimum absolute difference between x and y.}
#' \item{max_abs_diff}{The maximum absolute difference between x and y.}
#' \item{max_dec}{The maximum number of decimals to which a comparison makes sense.
#'  For example, if x contains only values up to the third decimals, and y is a
#'  normal double, max_dec will be three.}
#' \item{are_equal}{The maximal number of decimals to which all elements of x and y
#'  agree in absolute value. The comparison is on magnitudes, so two elements that
#'  are equal in size but opposite in sign count as agreeing; signed disagreements
#'  are reflected in `diff` and the mean / median / min / max absolute differences.
#'  `0` means the two agree in their integer parts but in no decimal place. `NA`
#'  means there is no agreement at all: either they already differ in their integer
#'  parts, or `na.rm = FALSE` and an element is missing.}
#' \item{diff_corres}{The number of differing variable-to-factor correspondences
#'  between x and y, when only the highest loading is considered. `NA` whenever the
#'  correspondences were not compared: for vector input, for a matrix with a single
#'  column, with `corres = FALSE`, and when a loading is missing under
#'  `na.rm = FALSE`.}
#' \item{diff_corres_cross}{The number of differing variable-to-factor correspondences
#'  between x and y when all loadings `>= thresh` are considered. `NA` under the same
#'  conditions as `diff_corres`.}
#' \item{g}{The root mean squared distance (RMSE) between x and y.}
#' \item{settings}{List of the settings used.}
#'
#' @seealso [efa_fit()] for the solutions being compared, and [efa_procrustes()] to rotate
#'  one solution onto another before comparing.
#'
#' @family factor comparison
#'
#' @export
#'
#' @examples
#' # A type SPSS EFA to mimick the SPSS implementation
#' EFA_SPSS_6 <- efa_fit(test_models$case_11b$cormat, n_factors = 6,
#'                       estimate_control = estimate_control(type = "SPSS"),
#'                       rotate_control = rotate_control(type = "SPSS"))
#'
#' # A type psych EFA to mimick the psych::fa() implementation
#' EFA_psych_6 <- efa_fit(test_models$case_11b$cormat, n_factors = 6,
#'                        estimate_control = estimate_control(type = "psych"),
#'                        rotate_control = rotate_control(type = "psych"))
#'
#' # compare the two
#' efa_compare(EFA_SPSS_6$unrot_loadings, EFA_psych_6$unrot_loadings,
#'             x_labels = c("SPSS", "psych"))
efa_compare <- function(x,
                        y,
                        reorder = c("congruence", "names", "none"),
                        corres = TRUE,
                        thresh = .3,
                        digits = 4,
                        m_red = .001,
                        range_red = .001,
                        round_red = 3,
                        print_diff = TRUE,
                        na.rm = FALSE,
                        x_labels = c("x", "y"),
                        plot = TRUE,
                        plot_red = .01)  {

  reorder <- .match_arg_ci(reorder)
  checkmate::assert_flag(corres)
  checkmate::assert_number(thresh, lower = 0)
  checkmate::assert_count(digits)
  checkmate::assert_number(m_red)
  checkmate::assert_number(range_red)
  checkmate::assert_number(round_red)
  checkmate::assert_flag(print_diff)
  checkmate::assert_flag(na.rm)
  checkmate::assert_character(x_labels, len = 2)
  checkmate::assert_flag(plot)
  checkmate::assert_number(plot_red)

  # reclass data.frames and tibbles to matrices so the stats functions afterwards
  # work
  if ((inherits(x, c("loadings", "LOADINGS", "SLLOADINGS", "matrix"))) &&
      (inherits(y, c("loadings", "LOADINGS", "SLLOADINGS", "matrix")))) {

    if (inherits(x, c("loadings", "LOADINGS", "SLLOADINGS"))) {
      x <- unclass(x)
    }

    if (inherits(y, c("loadings", "LOADINGS", "SLLOADINGS"))) {
      y <- unclass(y)
    }

    if (!is.numeric(x) || !is.numeric(y)) {
      cli::cli_abort(
        "{.arg x} and {.arg y} must be numeric vectors or matrices.",
        class = "efa_compare_bad_input"
      )
    }

    # check if dimensions match:
    if (any(dim(x) != dim(y))) {

      cli::cli_abort("{.arg x} and {.arg y} have different dimensions; {.fun efa_compare} only works with identical dimensions.",
                     class = "efa_compare_dim_mismatch")

    }

  } else if (inherits(x, c("numeric", "integer")) &&
             inherits(y, c("numeric", "integer"))) {

    if (length(x) != length(y)) {

      cli::cli_abort("{.arg x} and {.arg y} have different lengths; {.fun efa_compare} only works with identical dimensions.",
                     class = "efa_compare_dim_mismatch")

    }

  } else {

    # A fitted solution is the commonest thing to hand over here, so name the component that
    # holds its loadings rather than leaving the reader to find it.
    cli::cli_abort(
      c("{.arg x} ({.cls {class(x)}}) and {.arg y} ({.cls {class(y)}}) must be numeric vectors or matrices.",
        "i" = "From an {.cls efa} solution, use {.code $rot_loadings}, or {.code $unrot_loadings}
               if it was fitted without a rotation."),
      class = "efa_compare_bad_input")

  }

  if (length(x) == 0L) {
    cli::cli_abort(
      "{.arg x} and {.arg y} must contain at least one value.",
      class = "efa_compare_empty"
    )
  }

  # Matrices headed for congruence reordering are exempt: that branch reports a
  # non-finite input as a reordering failure, which names the step that actually
  # cannot proceed. The exemption has to cover every matrix the branch handles,
  # single-column ones included, or the same input reports differently depending on
  # how many columns it has.
  if ((any(is.infinite(x)) || any(is.infinite(y))) &&
      !(is.matrix(x) && identical(reorder, "congruence"))) {
    cli::cli_abort(
      "{.arg x} and {.arg y} must not contain infinite values.",
      class = "efa_compare_nonfinite"
    )
  }

  if (inherits(x, "matrix")) {

    n_factors <- ncol(x)

    if (reorder == "congruence") {
      # Match y's columns to those of x with the shared congruence alignment: an
      # optimal one-to-one assignment that maximises the matched absolute Tucker
      # congruences and reflects signs accordingly. The linear assignment
      # guarantees a permutation (a greedy row-wise which.max could send two
      # x-factors to the same y column, duplicating one column and dropping
      # another), and carries x's dimnames onto the realigned y. A single column
      # has nothing to permute but the same arbitrary sign as any other, so it
      # goes through the alignment too: two one-factor solutions that differ only
      # in the sign of their column are the same solution.
      #
      # Tucker's congruence is undefined for missing or non-finite inputs and for
      # zero/near-zero or non-finite congruences. Reject the former up front; map
      # the alignment's classed aborts for the latter to the same reorder error.
      if (!all(is.finite(x)) || !all(is.finite(y))) {
        cli::cli_abort(
          c("{.arg x} or {.arg y} contains missing or non-finite values; cannot reorder columns by congruence.",
            "i" = "Try another reordering method."),
          class = "efa_compare_congruence_na"
        )
      }

      congruence_undefined <- function(e) {
        cli::cli_abort(
          c("Tucker's congruence is undefined for a zero, near-zero, or non-finite column; cannot reorder columns by congruence.",
            "i" = "Try another reordering method."),
          class = "efa_compare_congruence_na"
        )
      }
      y <- tryCatch(
        .align_solution(L_target = x, L = y)$loadings,
        efa_zero_column = congruence_undefined,
        efa_undefined_congruence = congruence_undefined
      )
    } else if (reorder == "names" && n_factors > 1) {

      if (!is.null(colnames(x)) && !is.null(colnames(y))) {
        x <- x[, order(colnames(x))]
        y <- y[, order(colnames(y))]

        if(!all(colnames(x) == colnames(y))) {
          cli::cli_warn("{.code reorder = \"names\"} was used but the colnames of {.arg x} and {.arg y} differ; results might be inaccurate.",
                        class = "efa_compare_reorder_mismatch")
        }
      } else if (is.null(colnames(x)) || is.null(colnames(y))) {
        cli::cli_warn("{.arg reorder} was {.val names} but at least one of {.arg x} and {.arg y} is unnamed; proceeding without reordering.",
                      class = "efa_compare_reorder_unnamed")
      }

    }

  } else if (inherits(x, c("numeric", "integer"))) {

    if (reorder == "congruence"){

      # Congruence reordering needs columns to permute, so it does nothing here.
      # That is true of any vector, named or not, and the default is
      # `"congruence"` -- so the warning does not depend on the names. Only the
      # pointer to the alternative does: `reorder = "names"` is no use to someone
      # whose vectors carry no names.
      msg <- "{.arg reorder} was {.val congruence}, which only works for matrices; proceeding without reordering."
      if (!is.null(names(x)) && !is.null(names(y))) {
        msg <- c(msg, "i" = "To reorder vectors, set {.code reorder = \"names\"}.")
      }
      cli::cli_warn(msg, class = "efa_compare_reorder_vectors")

    } else if (reorder == "names") {

      if (!is.null(names(x)) && !is.null(names(y))) {

        x <- x[order(names(x))]
        y <- y[order(names(y))]

        if (!all(names(x) == names(y))) {
          cli::cli_warn("{.code reorder = \"names\"} was used but the names of {.arg x} and {.arg y} differ; results might be inaccurate.",
                        class = "efa_compare_reorder_mismatch")
        }



      } else if (is.null(names(x)) || is.null(names(y))) {
        cli::cli_warn("{.arg reorder} was {.val names} but at least one of {.arg x} and {.arg y} is unnamed; proceeding without reordering.",
                      class = "efa_compare_reorder_unnamed")
      }

      }

  }

  # delegate the difference statistics and factor correspondences to the
  # printless core; x and y are already coerced and aligned at this point, so the
  # core does not (and must not) reorder them.
  core <- .compare_loadings(x, y, thresh = thresh, na.rm = na.rm, corres = corres)

  settings <- list(
    reorder = reorder,
    corres = corres,
    digits = digits,
    thresh = thresh,
    m_red = m_red,
    range_red = range_red,
    round_red = round_red,
    print_diff = print_diff,
    na.rm = na.rm,
    x_labels = x_labels,
    plot_red = plot_red
  )



  # create output list: the core statistics followed by the recorded settings
  out <- c(core, list(settings = settings))

  # the superseded COMPARE() name keeps its old class string alongside the new one,
  # so `inherits(x, "COMPARE")` in existing code still resolves
  class(out) <- c("efa_compare", "COMPARE")

  out
}


# Difference statistics and factor correspondences for two already coerced and
# aligned loading/communality objects. Kept free of coercion, reordering, and
# printing so a caller that has aligned x and y to a shared target can reuse it
# without triggering a second alignment.
.compare_loadings <- function(x, y, thresh = 0.3, na.rm = FALSE, corres = TRUE,
                              decimals = TRUE) {

  # Factor correspondences need a matrix with more than one column, an enabled
  # comparison, and no missing loading to work around. Wherever one of those is
  # absent the correspondences were not compared, and the result must say so with
  # NA: a count of 0 would read as "they were compared and agreed everywhere".
  if (inherits(x, "matrix")) {
    if (ncol(x) > 1 && isTRUE(corres)) {
      if (!na.rm && (anyNA(x) || anyNA(y))) {
        diff_corres <- NA_integer_
        diff_corres_cross <- NA_integer_
      } else {
        corres_list <- .factor_corres(x, y, thresh = thresh)
        diff_corres <- corres_list$diff_corres
        diff_corres_cross <- corres_list$diff_corres_cross
      }
    } else {
      # a single factor leaves nothing to disagree on, and `corres = FALSE` asked
      # for no comparison at all
      diff_corres <- NA_integer_
      diff_corres_cross <- NA_integer_
    }
  } else {
    # a vector has no factors to correspond to
    diff_corres <- NA_integer_
    diff_corres_cross <- NA_integer_
  }

  # compute differences and statistics
  diff <- x - y

  # RMSE of the elementwise differences. na.rm drops missing entries from both
  # the numerator and the denominator so the statistic stays a proper mean
  # square; the unified form also covers the matrix case (length(diff) is the
  # full element count) without forming t(diff) %*% diff.
  sq <- diff ^ 2
  n_ok <- if (na.rm) sum(!is.na(diff)) else length(diff)
  if (na.rm && n_ok == 0L) {
    g <- mean_abs_diff <- median_abs_diff <- min_abs_diff <- max_abs_diff <- NA_real_
  } else {
    g <- sqrt(sum(sq, na.rm = na.rm) / n_ok)
    mean_abs_diff <- mean(abs(diff), na.rm = na.rm)
    median_abs_diff <- stats::median(abs(diff), na.rm = na.rm)
    min_abs_diff <- min(abs(diff), na.rm = na.rm)
    max_abs_diff <- max(abs(diff), na.rm = na.rm)
  }

  # max_dec: the most decimal places at which a comparison is meaningful (the
  # fewest decimals carried by either input). are_equal: the most decimal places
  # to which every corresponding pair of x and y still agrees. Comparing the
  # truncated values place by place is correct whatever the magnitude of the
  # integer part; concatenating the digit strings would let extra integer digits
  # inflate the count. signif() strips the floating-point representation noise
  # before truncating, so a value such as 0.57 (held as 0.5699999999999999) is
  # not mis-truncated to 0.56 at the second decimal.
  #
  # The two are the only non-trivial work in this function -- a decimal-place scan
  # of every element, then a truncation loop over the decimal places -- and only
  # the printed comparison reports them, so `decimals = FALSE` skips both (leaving
  # them NA) for a caller that needs the difference summaries alone.
  max_dec <- NA_real_
  are_equal <- NA_real_
  if (isTRUE(decimals)) {
    max_dec <- min(c(.decimals(x), .decimals(y)))

    ax <- abs(x)
    ay <- abs(y)
    if (anyNA(diff) && (!na.rm || all(is.na(diff)))) {
      # are_equal is undefined when there is nothing to compare: under na.rm = FALSE
      # any missing element poisons the comparison, and under na.rm = TRUE a fully
      # missing overlap leaves no comparable pair (the loop's all(logical(0)) would
      # otherwise spuriously report agreement to every decimal place). Report NA (the
      # printed "none"), mirroring how mean_abs_diff and g propagate NA, and matching
      # max_dec, which ignores missings.
      are_equal <- NA_real_
    } else {
      # Walk the decimal places from the integer part outwards, stopping at the
      # first one that differs: once a place disagrees no deeper place can count, so
      # there is nothing to gain from testing the rest. The counter starts at NA and
      # is only ever set to a place that was reached, so a comparison that already
      # fails at d = 0 (the integer parts differ) stays NA and is distinguishable
      # from one that succeeds at d = 0 but at no decimal place.
      for (d in 0:max_dec) {
        if (!isTRUE(all(trunc(signif(ax * 10^d, 13)) == trunc(signif(ay * 10^d, 13)),
                        na.rm = na.rm))) {
          break
        }
        are_equal <- as.double(d)
      }
    }
  }

  list(
    diff = diff,
    mean_abs_diff = mean_abs_diff,
    median_abs_diff = median_abs_diff,
    min_abs_diff = min_abs_diff,
    max_abs_diff = max_abs_diff,
    max_dec = max_dec,
    are_equal = are_equal,
    diff_corres = diff_corres,
    diff_corres_cross = diff_corres_cross,
    g = g
  )
}