R/efa_hull.R

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

#' Hull method for determining the number of factors to retain
#'
#' Implementation of the Hull method suggested by Lorenzo-Seva, Timmerman,
#' and Kiers (2011), with an extension to principal axis factoring. See details for
#' parallelization.
#'
#' @param x matrix or data.frame. Dataframe or matrix of raw data or matrix with
#' correlations.
#' @param N numeric. Number of cases in the data. This is passed to [efa_parallel].
#'  Only has to be specified if x is a correlation matrix, otherwise it is determined
#'  based on the dimensions of x.
#' @param n_fac_theor numeric. Theoretical number of factors to retain. One plus
#'   the larger of this number and the number of factors suggested by [efa_parallel]
#'   is used as the upper bound *J* of factors to extract in the Hull method.
#' @param estimator character. The estimator to use. One of  `"PAF"`,
#'    `"ULS"`, or  `"ML"`, for principal axis factoring, unweighted
#'    least squares, and maximum likelihood, respectively. Default is `"PAF"`.
#' @param gof character. The goodness of fit index to use. Either `"CAF"`,
#'   `"CFI"`, or `"RMSEA"`, or any combination of them.
#'   With the `"PAF"` estimator, only
#'   the CAF can be used as goodness of fit index. For details on the CAF, see
#'   Lorenzo-Seva, Timmerman, and Kiers (2011).
#' @param eigen_type character. On what the eigenvalues should be found in the
#'  parallel analysis. Can be one of `"SMC"`, `"PCA"`, or `"EFA"`.
#'   If using  `"SMC"` (default), the diagonal of the correlation matrices is
#'    replaced by the squared multiple correlations (SMCs) of the indicators. If
#'     using  `"PCA"`, the diagonal values of the correlation
#'  matrices are left to be 1. If using  `"EFA"`, eigenvalues are found on the
#'  correlation  matrices with the final communalities of an EFA solution as
#'  diagonal. This is passed to  [efa_parallel()].
#' @param use character. Passed to [stats::cor()] if raw data
#' is given as input. Default is `"pairwise.complete.obs"`.
#' @param cor_method character. One of `"pearson"`, `"spearman"`, or `"kendall"`,
#'   passed to [stats::cor()]. `"poly"` and `"tetra"` are not supported because
#'   `HULL` derives its factor-search bound from an internal parallel analysis
#'   against continuous reference data.
#'  Default is  `"pearson"`.
#' @param n_datasets numeric. The number of datasets to simulate. Must be at
#'   least 1. Default is 1000. This is passed to [efa_parallel()].
#' @param percent numeric. The percentile to take from the simulated eigenvalues.
#'  Default is 95. This is passed to [efa_parallel()].
#' @param decision_rule character. Which rule to use to determine the number of
#' factors to retain. Default is `"means"`, which will use the average
#' simulated eigenvalues. `"percentile"`, uses the percentiles specified
#' in percent. `"crawford"` uses the 95th percentile for the first factor
#' and the mean afterwards (based on Crawford et al, 2010). This is passed to [efa_parallel()].
#' @param n_factors numeric. Number of factors to extract if  `"EFA"` is
#' included in `eigen_type`. Default is 1. This is passed to
#' [efa_parallel()].
#' @param estimate_control an [estimate_control()] object with the estimation settings for the
#'  [efa_fit()] fits of the 0 to *J* factor solutions, and for the fit inside the internal
#'  [efa_parallel()] call. `NULL` (default) uses the [efa_fit()] defaults. This object carries
#'  estimation settings only; the fits are always unrotated, which the hull statistics (CFI,
#'  RMSEA, CAF) do not depend on.
#' @param ... Further arguments passed to [efa_fit()], also in
#' [efa_parallel()]. The estimation tuning knobs are not passed here; they live in
#' `estimate_control`, a rotation setting is not accepted because the fits are unrotated, and
#' neither are the standard-error arguments (`se`, `b_boot`, `ci`, `seed`), because the fits
#' are internal steps whose standard errors are not reported.
#'
#' @details The Hull method aims to find a model with an optimal balance between
#'  model fit and number of parameters, retaining only major factors
#'  (Lorenzo-Seva, Timmerman, & Kiers, 2011). It fits 0 to *J* factors -- where
#'  *J* is the number of factors suggested by parallel analysis (or `n_fac_theor`,
#'  if that is larger), plus one -- keeps the solutions on the upper boundary of
#'  the convex hull of goodness-of-fit against degrees of freedom, and selects the
#'  one at the sharpest elbow, i.e. with the highest *st* value.
#'
#'  Because it trades fit against parsimony instead of testing against a null model
#'  of uncorrelated variables, the Hull method does not lose accuracy for the
#'  correlated-factor structures where parallel analysis ([efa_parallel()]) tends to
#'  under-extract; the CAF variant in particular was among the more accurate criteria
#'  in Auerswald and Moshagen (2019). It needs at least six indicators and fits a
#'  model at every candidate factor count, so it is comparatively slow and is not an
#'  option for very short scales.
#'
#' The [efa_parallel] function and the principal axis factoring of the
#'   different number of factors can be parallelized using the future framework,
#'   by calling the [future::plan()] function. The examples
#'    provide example code on how to enable parallel processing.
#'
#'   The upper bound *J* comes from [efa_parallel()], which compares against simulated
#'   data, so the suggested number of factors varies slightly from run to run; a
#'   criterion-based rotation passed through `...` adds its own random starts. Call
#'   [base::set.seed()] beforehand to make a run reproducible; the result is then also
#'   independent of the parallel plan.
#'
#'   Note that if `gof = "RMSEA"` is used, 1 - RMSEA is actually used to
#'   compare the different solutions. This is necessary due to how the heuristic to
#'   locate the elbow of the hull works.
#'
#'   The solutions are fitted without inequality constraints, so a solution can be
#'   inadmissible (a Heywood case, or a fit that did not converge). Only the selected
#'   solution is checked for this; if it is inadmissible a warning is raised and the
#'   retained number of factors should be interpreted with caution.
#'
#'   The ML estimation method uses the [psych::fa()]
#'    starting values. See also the [efa_fit] documentation.
#'
#' @returns An object of class `efa_retention` (see [print.efa_retention()] and
#'   [plot.efa_retention()] for the print and plot methods). Its main fields are:
#' \item{n_factors}{A named numeric vector with the suggested number of factors
#'   for each requested goodness-of-fit index (`"CAF"`, `"CFI"`, and/or
#'   `"RMSEA"`).}
#' \item{results}{A list with one record per goodness-of-fit index, each holding
#'   the goodness-of-fit values, the degrees of freedom, the hull membership, and
#'   the retained solution used for printing and plotting. Each record also carries
#'   `st`, the elbow sharpness of every solution (see details): the retained
#'   solution has the largest value, and the runner-up shows how close the
#'   selection was. `st` is `NA` for the solutions where it is undefined, that is
#'   for those not on the hull and for the two hull endpoints, which have no
#'   neighbouring hull solution on one side. When fewer than three solutions
#'   remain on the hull, `st` is undefined throughout and the whole vector is `NA`;
#'   the retained solution is then the one with the highest goodness of fit, and a
#'   warning says so.}
#' \item{settings}{A list of the settings used, including `n_fac_max`, the upper
#'   bound *J* of the number of factors to extract (see details). For backwards
#'   compatibility the estimator is also repeated in `settings$method`.}
#'
#' For backwards compatibility the per-index suggestions are additionally available
#' as the top-level fields `n_fac_CAF`, `n_fac_CFI` and `n_fac_RMSEA`, each `NA` if
#' that index was not requested in `gof`. New code should read them from
#' `n_factors` instead.
#'
#' @source Auerswald, M., & Moshagen, M. (2019). How to determine the number of
#' factors to retain in exploratory factor analysis: A comparison of extraction
#' methods under realistic conditions. Psychological Methods, 24(4), 468–491.
#' https://doi.org/10.1037/met0000200
#'
#' @source Lorenzo-Seva, U., Timmerman, M. E., & Kiers, H. A. (2011).
#' The Hull method for selecting the number of common factors. Multivariate
#' Behavioral Research, 46(2), 340-364.
#'
#' @family factor retention criteria
#'
#' @seealso [efa_retain()] as a wrapper function for this and the other factor
#'   retention criteria.
#'
#' @export
#'
#' @examples
#' \donttest{
#' # using PAF (this will print a message if gof is not specified manually
#' # and CAF will be used automatically)
#' efa_hull(test_models$baseline$cormat, N = 500, gof = "CAF", n_datasets = 100)
#'
#' # using ML with all available fit indices (CAF, CFI, and RMSEA)
#' efa_hull(test_models$baseline$cormat, N = 500, estimator = "ML", n_datasets = 100)
#'
#' # using ULS with only RMSEA
#' efa_hull(test_models$baseline$cormat, N = 500, estimator = "ULS", gof = "RMSEA",
#'          n_datasets = 100)
#'}
#'
#'\dontrun{
#' # using parallel processing (Note: plans can be adapted, see the future
#' # package for details). future::plan() returns the plan it replaces, so
#' # on.exit() puts the session back as it was -- also if the call fails.
#' local({
#'   old_plan <- future::plan(future::multisession, workers = 2)
#'   on.exit(future::plan(old_plan), add = TRUE)
#'   efa_hull(test_models$baseline$cormat, N = 500, gof = "CAF")
#' })
#' }
efa_hull <- function(x, N = NA, n_fac_theor = NA,
                 estimator = c("PAF", "ULS", "ML"), gof = c("CAF", "CFI", "RMSEA"),
                 eigen_type = c("SMC", "PCA", "EFA"),
                 use = c("pairwise.complete.obs", "all.obs", "complete.obs",
                         "everything", "na.or.complete"),
                 cor_method = c("pearson", "spearman", "kendall", "poly", "tetra"),
                 n_datasets = 1000, percent = 95,
                 decision_rule = c("means", "percentile", "crawford"),
                 n_factors = 1, estimate_control = NULL, ...) {
  # Perform hull method following Lorenzo-Seva, Timmerman, and Kiers (2011)

  .reject_flat_knobs(...names(), fn = "efa_hull")
  .reject_unknown_fit_dots(...names(), fn = "efa_hull", unrotated = TRUE)
  .reject_rotation_dots(list(...), fn = "efa_hull")
  .assert_cor_input(x)

  estimator <- .match_arg_ci(estimator)
  use <- .match_arg_ci(use)
  cor_method <- .match_arg_ci(cor_method)
  # The Hull method derives its factor-search bound from an internal parallel
  # analysis, whose reference data are continuous; poly/tetra are therefore not
  # supported, consistent with PARALLEL/NEST/CD.
  .reject_poly_reference(cor_method, "efa_hull")
  gof <- .match_arg_ci(gof, several.ok = TRUE)
  eigen_type <- .match_arg_ci(eigen_type)
  .assert_args({
    checkmate::assert_count(n_fac_theor, na.ok = TRUE)
    checkmate::assert_count(N, na.ok = TRUE)
    decision_rule <- .match_arg_ci(decision_rule)
    .assert_estimate_control(estimate_control)
    checkmate::assert_count(n_factors)
    checkmate::assert_count(n_datasets, positive = TRUE)
    checkmate::assert_number(percent, lower = 0, upper = 100)
  })

  if (ncol(x) < 6) {
    cli::cli_abort(
      c("The data has fewer than 6 indicators.",
        "i" = "The Hull method needs at least 6."),
      class = "efa_hull_min_indicators"
    )
  }

  if (estimator == "PAF" && !all(gof == "CAF")) {
    cli::cli_inform(
      c("i" = 'Only CAF can be used as gof if estimator "PAF" is used. Setting gof to "CAF"'),
      class = "efa_hull_gof_caf"
    )
    gof <- "CAF"
  }

  # Detect or compute the correlation matrix, check it, and smooth it if needed
  prep <- .prepare_cor_input(
    x, N = N, use = use, cor_method = cor_method, N_policy = "required",
    singular_tail = "the Hull method cannot be executed",
    N_required_msg = "{.arg N} is not specified but is needed to compute some fit indices.")
  R <- prep$R
  N <- prep$N

  m <- ncol(R)

  # 1) perform parallel analysis to find J as n_fac_theor + 1
  par_res <- efa_parallel(R, N = N, eigen_type = eigen_type, estimator = estimator,
                      n_datasets = n_datasets, percent = percent,
                      decision_rule = decision_rule, n_factors = n_factors,
                      estimate_control = estimate_control, ...)

  n_fac_PA <- unname(par_res$n_factors[eigen_type])

  if (is.na(n_fac_PA)) {

    if (!is.na(n_fac_theor)) {
      J <- n_fac_theor + 1
    } else {
      J <- .det_max_factors(ncol(R))
    }

  } else {

    J <- max(c(n_fac_PA, n_fac_theor), na.rm = TRUE) + 1

  }

  if (J > .det_max_factors(ncol(R))) {
    J <- .det_max_factors(ncol(R))
    cli::cli_warn("Setting the maximum number of factors to {J} to ensure overidentified models.",
                  class = "efa_hull_max_factors")
  }

  # The Hull method needs at least three solutions (for 0 to J factors) to form a
  # hull, i.e. J of at least 2. With few indicators the largest over-identified
  # model can itself be below the usual floor of three factors, so cap the floor at
  # that maximum to avoid forcing an under-identified (df = 0) top model.
  hull_floor <- min(3L, .det_max_factors(ncol(R)))
  if (J < hull_floor) {
    cli::cli_warn(
      c("The suggested maximum number of factors was {J}, but the Hull method needs at least {hull_floor}.",
        "i" = "Setting it to {hull_floor}."),
      class = "efa_hull_min_factors"
    )
    J <- hull_floor

  }

  # 2) perform factor analysis for the range of dimensions 1:J and compute f and
  #    df for every solution
  s <- matrix(0, ncol = 4, nrow = J + 1)
  s[, 1] <- 0:J

  # first for 0 factors
  if ("CAF" %in% gof) {
    s_CAF <- s
    colnames(s_CAF) <- c("nfactors", "CAF", "df", "st")
    s_CAF[1, 2] <- 1 - .compute_kmo(R)$KMO
    s_CAF[1, 3] <- (m**2 - m) / 2
  }

  if ("CFI" %in% gof) {
    s_CFI <- s
    colnames(s_CFI) <- c("nfactors", "CFI", "df", "st")
    s_CFI[1, 2] <- 0
    s_CFI[1, 3] <- (m**2 - m) / 2
  }

  if ("RMSEA" %in% gof) {
    s_RMSEA <- s
    colnames(s_RMSEA) <- c("nfactors", "RMSEA", "df", "st")
    # 0-factor (independence model) reference, on the same uncorrected (N - 1) discrepancy scale
    # the RMSEA of the 1:J solutions uses in .chi_fit_indices(). When N is too small for the
    # Bartlett correction (the corrected null chi-square is NA) the chi-square asymptotics break
    # down, so the reference is left undefined too -- matching how .gof()/SMT() drop the RMSEA at
    # such N (where every fitted 1:J solution is likewise NA, the factor-count term only lowering
    # the multiplier further). The log-determinant is computed once and reused for both calls.
    ld_R <- determinant(R, logarithm = TRUE)
    df <- (m**2 - m) / 2
    # compute 1 - RMSEA
    s_RMSEA[1, 2] <- if (is.na(.null_chisq(R, N, ld = ld_R))) {
      NA_real_
    } else {
      1 - .rmsea_point(.null_chisq(R, N, ld = ld_R, corrected = FALSE), df, N)
    }
    s_RMSEA[1, 3] <- (m**2 - m) / 2

  }

  # Calculate loadings with the EFA function. future.seed = TRUE because a
  # criterion-based rotation passed through `...` (e.g. oblimin) draws random starts, so
  # the workers need a managed random-number stream; without one an oblique HULL run
  # would not be reproducible from set.seed() under a parallel plan.
  loadings <- suppressWarnings(future.apply::future_lapply(seq_len(J), efa_fit,
                                                           x = R,
                                                           estimator = estimator,
                                                           N = N,
                                                           estimate_control = estimate_control,
                                                           ...,
                                                           future.seed = TRUE))

  # then for 1 to J factors. estimator == "PAF" forces gof to "CAF" above, so the
  # gof-keyed blocks already cover it; the df is the same for every index (Eq 4
  # gives the free parameters, and the difference in df equals the difference in
  # free parameters).
  for (i in seq_len(J)) {
    if ("CAF" %in% gof) {
      # compute goodness of fit "f" as CAF (common part accounted for; Eq 3)
      s_CAF[i + 1, 2] <- loadings[[i]]$fit_indices$CAF
      s_CAF[i + 1, 3] <- loadings[[i]]$fit_indices$df
    }

    if ("CFI" %in% gof) {
      s_CFI[i + 1, 2] <- loadings[[i]]$fit_indices$CFI
      s_CFI[i + 1, 3] <- loadings[[i]]$fit_indices$df
    }

    if ("RMSEA" %in% gof) {
      # compute 1 - RMSEA
      s_RMSEA[i + 1, 2] <- 1 - loadings[[i]]$fit_indices$RMSEA
      s_RMSEA[i + 1, 3] <- loadings[[i]]$fit_indices$df
    }

  }

  out_CAF <- list(s_complete = NA, retain = NA)
  out_CFI <- list(s_complete = NA, retain = NA)
  out_RMSEA <- list(s_complete = NA, retain = NA)

  if("CAF" %in% gof) {
    out_CAF <- .hull_calc(s = s_CAF, J = J, gof_t = "CAF")
  }
  if("CFI" %in% gof) {
    out_CFI <- .hull_calc(s = s_CFI, J = J, gof_t = "CFI")
  }
  if("RMSEA" %in% gof) {
    out_RMSEA <- .hull_calc(s = s_RMSEA, J = J, gof_t = "RMSEA")
  }

  gof_results <- list(CAF = out_CAF, CFI = out_CFI, RMSEA = out_RMSEA)

  # one record per requested goodness-of-fit index (df vs. fit, with the hull
  # membership and the retained solution used for printing and plotting)
  results <- list()
  inadmissible <- character(0)
  for (g in c("CAF", "CFI", "RMSEA")) {
    if (!(g %in% gof)) next
    sol <- gof_results[[g]]$s_complete
    retain <- gof_results[[g]]$retain

    # The right-most hull solutions over-extract and may be inadmissible, so flag
    # only the selected solution if it has a Heywood case or did not converge.
    if (!is.na(retain) && retain >= 1) {
      fit_sel <- loadings[[retain]]
      issues <- c(if (length(fit_sel$heywood) > 0) "Heywood case",
                  if (isTRUE(fit_sel$convergence != 0)) "non-convergence")
      if (length(issues) > 0) {
        inadmissible <- c(
          inadmissible,
          paste0(g, ": ", .retention_count(retain), " factor",
                 if (retain != 1) "s" else "",
                 " (", paste(issues, collapse = ", "), ")")
        )
      }
    }

    results[[g]] <- list(
      name = g,
      label = g,
      n_factors = retain,
      plot_type = "hull",
      x = unname(sol[, "df"]),
      y = unname(sol[, g]),
      reference = NULL,
      threshold = NULL,
      highlight = retain,
      point_labels = unname(sol[, "nfactors"]),
      on_hull = gof_results[[g]]$on_hull,
      # elbow sharpness of each solution (Eq 5); the selected one has the largest
      # value, so the runner-up shows how close the decision was
      st = unname(sol[, "st"])
    )
  }

  if (length(inadmissible) > 0) {
    cli::cli_warn(
      c("The Hull method selected an inadmissible solution: {inadmissible}.",
        "i" = "The selected solution has a Heywood case or did not converge, so the retained number of factors may be unreliable; interpret it with caution and cross-check with other criteria."),
      class = "efa_hull_inadmissible"
    )
  }

  out <- .new_efa_retention(
    "HULL",
    results = unname(results),
    settings = list(N = N,
                    estimator = estimator,
                    # back-compat alias, as for the frozen P_type key
                    method = estimator,
                    gof = gof,
                    n_fac_theor = n_fac_theor,
                    eigen_type = eigen_type,
                    use = use,
                    cor_method = cor_method,
                    n_fac_max = J),
    subtitle = paste0("Estimator: ", estimator)
  )

  # back-compat aliases of the per-index suggestions, as for the frozen method key
  out$n_fac_CAF <- unname(out$n_factors["CAF"])
  out$n_fac_CFI <- unname(out$n_factors["CFI"])
  out$n_fac_RMSEA <- unname(out$n_factors["RMSEA"])

  return(out)

}


.hull_calc <- function(s, J, gof_t){

  # 3) sort n solutions by their df values and denoted by s (already done)

  # 4) all solutions s are excluded for which a solution sj (j<i) exists such
  #    that fj > fi (eliminate solutions not on the boundary of the convex hull)

  s_complete <- s

  # A non-finite goodness-of-fit value (e.g. an undefined CFI/RMSEA for a Heywood
  # or near-singular model) cannot lie on the convex hull; drop those solutions
  # with a classed warning rather than failing in the comparisons below.
  na_rows <- !is.finite(s[, 2])
  if (any(na_rows)) {
    cli::cli_warn(
      c("{sum(na_rows)} solution{?s} had a non-finite {gof_t} value and {?was/were} excluded from the hull.",
        "i" = "Inspect the affected models, or try a different goodness-of-fit index or estimation method."),
      class = "efa_hull_na_fit"
    )
    s <- s[!na_rows, , drop = FALSE]
  }

  if (nrow(s) < 1) {
    cli::cli_abort(
      c("No solution had a finite {gof_t} value, so the Hull method cannot proceed.",
        "i" = "Try a different goodness-of-fit index or estimation method."),
      class = "efa_hull_no_fit"
    )
  }

  d_s <- diff(s[, 2])
  while (any(d_s < 0)) {
    s <- s[c(1, d_s) > 0, , drop = FALSE]
    if(nrow(s) == 1){
      break
    }
    d_s <- diff(s[, 2])
  }

  # 5) all triplets of adjacent solutions are considered consecutively.
  #    The middle solution is excluded if its point is below or on the line
  #    connecting its neighbors in GOF vs df.

  # 6) repeat 5) until no solution can be excluded

  # The `i <= nr_s - 1` bound ensures the final interior triplet is also tested;
  # the loop is a no-op while fewer than three boundary solutions remain.
  nr_s <- nrow(s)
  i <- 2

  while(i <= nr_s - 1) {

    f1 <- s[i - 1, 2]
    f2 <- s[i, 2]
    f3 <- s[i + 1, 2]
    df1 <- s[i - 1, 3]
    df2 <- s[i, 3]
    df3 <- s[i + 1, 3]

    # compute f2 if it were on the line between f1 and f3
    p_f2 <- f1 + (f3 - f1) / (df3 - df1) * (df2 - df1)

    # check if f2 is below or on the predicted line and if so, remove it
    if (f2 <= p_f2) {
      s <- s[-i, , drop = FALSE]
      nr_s <- nr_s -1
      i <- 1
    }
    i <- i + 1
  }

  if (nrow(s) < 3) {
    cli::cli_warn(
      c("Fewer than three solutions were located on the hull using {gof_t} as goodness-of-fit index.",
        "i" = "Proceeding with the maximum-{gof_t} value as a heuristic; consider additional indices or methods as a robustness check."),
      class = "efa_hull_few_solutions"
    )

    # the st values are undefined with fewer than three hull solutions
    s_complete[, 4] <- NA_real_

    # 8) select solution with highest gof value
    retain <- s[which.max(s[, 2]), 1]


  } else {

    # 7) the st values of the hull solutions are determined (Eq 5)
    for (i in 2:(nrow(s) - 1)) {

      f_i <- s[i, 2]
      f_p <- s[i - 1, 2]
      f_n <- s[i + 1, 2]
      df_i <- s[i, 3]
      df_p <- s[i - 1, 3]
      df_n <- s[i + 1, 3]

      s[i, 4] <- ((f_i - f_p) / (df_i - df_p)) / ((f_n - f_i) / (df_n - df_i))

    }

    # Combine values. Only the interior hull solutions have an st value: the ratio
    # needs both a predecessor and a successor on the hull, so the two hull
    # endpoints are as undefined as the solutions that are not on the hull at all,
    # and are reported as NA. The local `s` still carries the zero its endpoints
    # were initialised with, which is why the reported column is rebuilt here
    # rather than copied wholesale from it.
    interior <- s[-c(1, nrow(s)), , drop = FALSE]
    s_complete[, 4] <- NA_real_
    s_complete[match(interior[, 1], s_complete[, 1]), 4] <- interior[, 4]

    # 8) select solution with highest st value
    retain <- s[which.max(s[, 4]), 1]

  }


  # hull membership: the solutions that survived the boundary/triplet elimination.
  # Derived from the surviving solutions rather than from `st`, so the retained
  # solution is still flagged as on the hull in the < 3 fallback where `st` is
  # undefined.
  on_hull <- s_complete[, 1] %in% s[, 1]

  out <- list(s_complete = s_complete,
              retain = unname(retain),
              on_hull = on_hull)

  return(out)

}