R/fviz_nbclust.R

In factoextra: Extract and Visualize the Results of Multivariate Data Analyses

#' @include  hcut.R
 NULL
#' Determining and Visualizing the Optimal Number of Clusters
#' @description Partitioning methods, such as k-means clustering require the 
#'   users to specify the number of clusters to be generated. \itemize{
#'   \item fviz_nbclust(): Determines and visualizes the optimal number of
#'   clusters using different methods: \strong{within-cluster sum of squares},
#'   \strong{average silhouette}, and the \strong{gap statistic}. Silhouette values
#'   are evaluated only for \code{k = 2, ..., k.max} because average
#'   silhouette width is undefined for a one-cluster partition. The
#'   \code{"silhouette"} and \code{"gap_stat"} plots mark their optimum with a
#'   dashed guide line; set \code{mark_optimal = TRUE} to also mark the elbow on
#'   the \code{"wss"} plot, or \code{mark_optimal = FALSE} to omit the guide line
#'   for every method.
#'   \item fviz_gap_stat(): Visualizes the gap statistic generated by the
#'   function \code{\link[cluster]{clusGap}}() [in cluster package]. The optimal
#'   number of clusters is specified using the
#'   \code{\link[cluster]{maxSE}} method with
#'   \code{method = "firstSEmax"}. }
#'
#'   For \code{method = "wss"}, \code{factoextra} computes the \code{k = 1}
#'   baseline internally so helper functions such as \code{hcut()} and
#'   \code{hkmeans()} can keep rejecting direct \code{k = 1} inputs.
#'
#'   Read more:
#'   \href{https://www.datanovia.com/learn/machine-learning/clustering/optimal-clusters}{Determining the Optimal Number of Clusters in R}.
#'
#' @details When \code{mark_optimal = TRUE}, the \code{"wss"} elbow is located with
#'   a deterministic chord-distance heuristic: both axes are rescaled to the unit
#'   interval \code{[0, 1]} (so the choice does not depend on the magnitude of the
#'   within-cluster sum of squares) and the marked \code{k} is the one whose point
#'   lies farthest from the straight line joining the first and last points of the
#'   curve. This is related to the chord-distance idea used in knee detection,
#'   but it does not implement the full Kneedle sensitivity and local-maximum
#'   procedure of \enc{Satopää}{Satopaa} et al. (2011). The heuristic always returns a candidate, even for data
#'   with no clear cluster structure, which is why it is opt-in; the
#'   \code{"silhouette"} and
#'   \code{"gap_stat"} guide lines instead mark defined optima (the maximum average
#'   silhouette width and the \code{\link[cluster]{maxSE}} location).
#'
#' @references \enc{Satopää}{Satopaa}, V., Albrecht, J., Irwin, D. & Raghavan, B. (2011). Finding
#'   a "Kneedle" in a Haystack: Detecting Knee Points in System Behavior.
#'   \emph{2011 31st International Conference on Distributed Computing Systems
#'   Workshops}, 166-171. \doi{10.1109/ICDCSW.2011.20}
#'
#' @param x numeric matrix or data frame. In the function fviz_nbclust(), x can
#'   be the results of the function NbClust(). For \code{method = "silhouette"}
#'   or \code{"wss"}, x may also be a precomputed dissimilarity (an object of
#'   class \code{"dist"}); it is then passed directly to \code{FUNcluster},
#'   which must be a dissimilarity-capable method (e.g. \code{cluster::pam} or
#'   \code{factoextra::hcut}). \code{method = "gap_stat"} still needs the raw data.
#' @param method the method to be used for estimating the optimal number of 
#'   clusters. Possible values are "silhouette" (for average silhouette width), 
#'   "wss" (for total within-cluster sum of squares), and "gap_stat" (for the gap statistic).
#' @param FUNcluster a partitioning function which accepts as first argument a
#'   (data) matrix like \code{x}, second argument, say \code{k >= 2}, the
#'   number of clusters desired, and returns a list with a component named
#'   \code{cluster} which contains the grouping of observations. Allowed values
#'   include: \code{kmeans}, \code{cluster::pam}, \code{cluster::clara},
#'   \code{cluster::fanny}, \code{hcut}, etc. In \code{method = "wss"} mode,
#'   \code{fviz_nbclust()} computes the \code{k = 1} baseline internally instead
#'   of calling \code{FUNcluster(x, 1, ...)}. This argument is not required
#'   when \code{x} is an output of the function \code{NbClust::NbClust()}.
#' @param diss dist object as produced by dist(), i.e.: diss = dist(x, method = 
#'   "euclidean"). Used to compute the average silhouette width and
#'   within-cluster sum of squares. If NULL, dist(x) is
#'   computed with the default method = "euclidean"
#' @param k.max the maximum number of clusters to consider, must be at least two.
#' @param nboot integer, number of Monte Carlo ("bootstrap") samples. Used only for determining the number of clusters 
#' using the gap statistic.
#' @param verbose logical value. If TRUE, progress information is printed.
#' @param barfill,barcolor fill color and outline color for bars
#' @param linecolor color for lines
#' @param print.summary logical value. If TRUE, the optimal number of clusters
#'   is printed in \code{fviz_nbclust()}.
#' @param mark_optimal logical, or \code{NULL} (default). \code{NULL} keeps each
#'   method's standard marker: a dashed guide line at the estimated optimal number
#'   of clusters is drawn for \code{method = "silhouette"} (maximum average
#'   silhouette width) and \code{method = "gap_stat"} (the
#'   \code{\link[cluster]{maxSE}} location), and omitted for \code{method = "wss"}.
#'   Set \code{TRUE} to also mark the \code{"wss"} elbow (a maximum-distance
#'   heuristic that returns a candidate even when the data has no clear cluster
#'   structure; see Details), or \code{FALSE} to omit the guide line for every
#'   method.
#' @param ... optionally further arguments:
#'   arguments for FUNcluster() in "wss"/"silhouette" modes; arguments for
#'   \code{\link[cluster]{clusGap}}() in "gap_stat" mode. A \code{maxSE} list
#'   can also be supplied in "gap_stat" mode and is forwarded to
#'   \code{fviz_gap_stat()}.
#'   
#' @return Both \code{fviz_nbclust()} and \code{fviz_gap_stat()} return a
#'   ggplot2 object.
#' @seealso \code{\link{fviz_cluster}}, \code{\link{eclust}}.
#'   Online tutorial: \href{https://www.datanovia.com/learn/machine-learning/clustering/optimal-clusters}{Determining the Optimal Number of Clusters in R}.
#' @author Alboukadel Kassambara \email{alboukadel.kassambara@@gmail.com}
#'   
#' @examples 
#' set.seed(123)
#' 
#' # Data preparation
#' # +++++++++++++++
#' data("iris")
#' head(iris)
#' # Remove species column (5) and scale the data
#' iris.scaled <- scale(iris[, -5])
#' 
#' 
#' # Optimal number of clusters in the data
#' # ++++++++++++++++++++++++++++++++++++++
#' # Examples are provided only for kmeans, but
#' # you can also use cluster::pam (for pam) or
#' #  hcut (for hierarchical clustering)
#'  
#' ### Elbow method (look at the knee)
#' # Elbow method for kmeans
#' fviz_nbclust(iris.scaled, kmeans, method = "wss") +
#' geom_vline(xintercept = 3, linetype = 2)
#'
#' # Let factoextra mark the elbow automatically
#' fviz_nbclust(iris.scaled, kmeans, method = "wss", mark_optimal = TRUE)
#'
#' # WSS with hierarchical clustering keeps the internal k = 1 baseline
#' fviz_nbclust(iris.scaled, hcut, method = "wss", hc_method = "complete")
#' 
#' # Average silhouette for kmeans
#' fviz_nbclust(iris.scaled, kmeans, method = "silhouette")
#' 
#' ### Gap statistic
#' library(cluster)
#' set.seed(123)
#' # Compute gap statistic for kmeans
#' # we used B = 10 for demo. Recommended value is ~500
#' gap_stat <- clusGap(iris.scaled, FUN = kmeans, nstart = 25,
#'  K.max = 10, B = 10)
#'  print(gap_stat, method = "firstSEmax")
#' fviz_gap_stat(gap_stat)
#'  
#' # Gap statistic for hierarchical clustering
#' gap_stat <- clusGap(iris.scaled, FUN = hcut, K.max = 10, B = 10)
#' fviz_gap_stat(gap_stat)
#' 
#'  
#' @name fviz_nbclust
#' @rdname fviz_nbclust
#' @export
fviz_nbclust <- function (x, FUNcluster = NULL, method = c("silhouette", "wss", "gap_stat"),
                          diss = NULL, k.max = 10, nboot = 100, verbose = interactive(),
                          barfill="steelblue", barcolor="steelblue",
                          linecolor = "steelblue", print.summary = TRUE,  ...,
                          mark_optimal = NULL)
  {
  k.max <- .coerce_integerish(k.max, "k.max", lower = 2L,
                              value_label = "single integer value")
  method = match.arg(method)
  # mark_optimal = NULL keeps each method's historical default: the "silhouette"
  # and "gap_stat" plots mark their optimum (as they always have), while "wss"
  # is left unmarked. TRUE/FALSE force the guide line on/off for every method.
  if(is.null(mark_optimal)) mark_optimal <- (method != "wss")
  # x may also be a precomputed dissimilarity ("dist") for diss-capable methods (#90).
  is_dist <- inherits(x, "dist")
  if(!inherits(x, c("data.frame", "matrix", "dist")) && !("Best.nc" %in% names(x)))
    stop("x should be an object of class matrix/data.frame/dist or ",
         "an object created by the function NbClust() [NbClust package].")
  
  # x is an object created by the function NbClust() [NbClust package]
  if(inherits(x, "list") && "Best.nc" %in% names(x)){
      best_nc <- x$Best.nc
      # FIX: R 4.0.0+ deprecation - class() can return multiple values for matrices
      # Using inherits() or is.X() functions instead of class() == comparison
      # See: https://github.com/kassambara/factoextra/issues/171
      if(is.numeric(best_nc) && !is.matrix(best_nc)) print(best_nc)
      else if(is.matrix(best_nc))
        .viz_NbClust(x, print.summary, barfill, barcolor)
  }
  else if(is.null(FUNcluster)) stop("The argument FUNcluster is required. ",
                                    "Possible values are kmeans, pam, hcut, clara, ...")
  else if(!is.function(FUNcluster)){
    stop(
      "The argument FUNcluster should be a function. ",
      "Check if you're not overriding the specified function name somewhere."
      )
  }
  else if(method %in% c("silhouette", "wss")) {

      if (is.data.frame(x)) x <- as.matrix(x)
      # When x is a dissimilarity, use it directly (don't recompute) and pass it
      # to FUNcluster, which must be a dissimilarity-capable method (#90).
      # Detect capability via a 'diss'/'isdiss' formal (pam/fanny/hcut have one;
      # kmeans/clara do not and would silently mis-cluster a coerced dist).
      if(is_dist){
        fmls <- tryCatch(names(formals(FUNcluster)), error = function(e) character(0))
        if(!any(c("diss", "isdiss") %in% fmls))
          stop("FUNcluster does not accept a distance matrix (it has no 'diss'/'isdiss' ",
               "argument). Use a dissimilarity-capable method such as cluster::pam, ",
               "cluster::fanny, or factoextra::hcut, or pass the raw data instead of a 'dist'.",
               call. = FALSE)
      }
      if(is.null(diss)) diss <- if(is_dist) x else stats::dist(x)
      .cluster_at <- function(i) FUNcluster(x, i, ...)

      v <- rep(0, k.max)
      if(method == "silhouette"){
        for(i in 2:k.max){
          clust <- .cluster_at(i)
          v[i] <- .get_ave_sil_width(diss, clust$cluster)
        }
      }
      else if(method == "wss"){
        n_obs <- attr(diss, "Size")
        if(is.null(n_obs) || !is.numeric(n_obs))
          stop("Unable to determine the number of observations from diss")
        # Compute the one-cluster baseline internally so callers may reject k = 1.
        v[1] <- .get_withinSS(diss, rep(1L, n_obs))
        for(i in 2:k.max){
          clust <- .cluster_at(i)
          v[i] <- .get_withinSS(diss, clust$cluster)
        }
        
      }
      
      if(method == "silhouette"){
        silhouette_clusters <- seq.int(2, k.max)
        df <- data.frame(clusters = as.factor(silhouette_clusters), y = v[silhouette_clusters])
      }
      else{
        df <- data.frame(clusters = as.factor(1:k.max), y = v)
      }
      
      ylab <- "Total Within Sum of Square"
      main_title <- "Optimal number of clusters"
      if(method == "silhouette") {
        ylab <- "Average silhouette width"
        main_title <- "Optimal number of clusters (method = \"silhouette\")"
      }

      p <- ggpubr::ggline(df, x = "clusters", y = "y", group = 1,
                          color = linecolor, ylab = ylab,
                          xlab = "Number of clusters k",
                          main = main_title
                          )
      if(mark_optimal) {
        if(method == "silhouette") {
          valid_clusters <- which(!is.na(v[-1])) + 1L
          if(length(valid_clusters) > 0) {
            best_k <- valid_clusters[which.max(v[valid_clusters])]
            p <- p + geom_vline(
              xintercept = match(as.character(best_k), levels(df$clusters)),
              linetype=2, color = linecolor
            )
          }
        }
        else { # wss: mark the elbow (max distance to the first-last chord)
          elbow_k <- .wss_elbow(seq_len(k.max), v)
          if(!is.na(elbow_k)) {
            p <- p + geom_vline(
              xintercept = match(as.character(elbow_k), levels(df$clusters)),
              linetype=2, color = linecolor
            )
          }
        }
      }

      return(p)
  }
  
  else if(method == "gap_stat"){
    if(is_dist)
      stop("method = 'gap_stat' requires the original data, not a distance matrix ",
           "(the gap statistic simulates reference datasets). Use method = 'silhouette' or ",
           "'wss' with a 'dist', or pass the raw data.", call. = FALSE)
    extra_args <- list(...)
    if(!is.null(extra_args$maxSE)) {
      maxSE <- extra_args$maxSE
      extra_args$maxSE <- NULL
    } else {
      maxSE <- list(method = "firstSEmax", SE.factor = 1)
    }
    gap_stat <- do.call(cluster::clusGap, c(list(x = x, FUNcluster = FUNcluster, K.max = k.max, B = nboot,
                                                  verbose = verbose), extra_args))
    p <- fviz_gap_stat(gap_stat,  linecolor = linecolor, maxSE = maxSE,
                       mark_optimal = mark_optimal)
    return(p)
  }
  

}

#' @rdname fviz_nbclust
#' @param gap_stat an object of class "clusGap" returned by the function
#'   clusGap() [in cluster package]
#' @param maxSE a list containing the parameters \code{method} and
#'   \code{SE.factor} used by \code{\link[cluster]{maxSE}} to locate the gap
#'   statistic optimum. The default is
#'   \code{list(method = "firstSEmax", SE.factor = 1)}. Allowed methods include:
#'   \itemize{ \item "globalmax": simply corresponds to the global maximum,
#'   i.e., is which.max(gap) \item "firstmax": gives the location of the first
#'   local maximum \item "Tibs2001SEmax": uses the criterion, Tibshirani et al
#'   (2001) proposed: "the smallest k such that gap(k) >= gap(k+1) - s(k+1)".
#'   It's also possible to use "the smallest k such that gap(k) >= gap(k+1) -
#'   SE.factor*s(k+1)" where SE.factor is a numeric value which can be 1
#'   (default), 2, 3, etc. \item "firstSEmax": location of the first f() value
#'   which is not larger than the first local maximum minus SE.factor * SE.f,
#'   i.e, within an "f S.E." range of that maximum. \item
#'   see \code{\link[cluster]{maxSE}} for more options }
#'
#' @section Method selection for gap statistic:
#'   The default \code{"firstSEmax"} method returns the first value within
#'   \code{SE.factor} standard errors of the first local maximum. Other
#'   \code{\link[cluster]{maxSE}} rules can be selected explicitly through
#'   \code{maxSE}.
#'
#'
#' @export
fviz_gap_stat <- function(gap_stat,  linecolor = "steelblue",
                          maxSE = list(method = "firstSEmax", SE.factor = 1),
                          mark_optimal = NULL){
  if(!inherits(gap_stat, "clusGap"))
    stop("Only an object of class clusGap is allowed. (cluster package)")
  # NULL keeps the historical behavior of marking the optimum; FALSE omits it.
  if(is.null(mark_optimal)) mark_optimal <- TRUE
  if(is.list(maxSE)){
    if(is.null(maxSE$method)) maxSE$method = "firstSEmax"
    if(is.null(maxSE$SE.factor)) maxSE$SE.factor = 1
  }
  else stop("The argument maxSE must be a list containing the parameters method and SE.factor")
  
  # first local max
  gap <- gap_stat$Tab[, "gap"]
  se <- gap_stat$Tab[, "SE.sim"]
  k <- .maxSE(gap, se, method = maxSE$method, SE.factor = maxSE$SE.factor)

  df <- as.data.frame(gap_stat$Tab)
  df$clusters <- as.factor(seq_len(nrow(df)))
  df$ymin <- gap-se
  df$ymax <- gap + se
  p <- ggpubr::ggline(df, x = "clusters", y = "gap", group = 1, color = linecolor)+
    # FIX: ggplot2 3.0.0+ deprecation - aes_string() replaced with aes() + .data pronoun
    # See: https://github.com/kassambara/factoextra/issues/190
    ggplot2::geom_errorbar(aes(ymin = .data[["ymin"]], ymax = .data[["ymax"]]), width=.2, color = linecolor)+
    labs(y = "Gap statistic (k)", x = "Number of clusters k",
         title = "Optimal number of clusters (method = \"gap_stat\")")
  if(mark_optimal)
    p <- p + geom_vline(xintercept = k, linetype=2, color = linecolor)
  p
}




# Get the average silhouette width
# ++++++++++++++++++++++++++
# Cluster package required
# d: dist object
# cluster: cluster number of observation
# Returns NA if silhouette cannot be computed (e.g., k <= 1 or k >= n)
# Fixes GitHub issues #113 and #147
.get_ave_sil_width <- function(d, cluster){
  if (!requireNamespace("cluster", quietly = TRUE)) {
    stop("cluster package needed for this function to work. Please install it.")
  }
  ss <- cluster::silhouette(cluster, d)
  # Handle case where silhouette() returns NA (when k <= 1 or k >= n)
  if (length(ss) == 1 && is.na(ss)) {
    return(NA_real_)
  }
  mean(ss[, 3])
}

# Get total within sum of square
# +++++++++++++++++++++++++++++
# d: dist object
# cluster: cluster number of observation
#
# OPTIMIZATION: Replaced explicit loops with vapply() for better performance
# - Pre-compute cluster indices once using split()
# - Use vapply for type-safe vectorized computation
# - Avoid repeated subsetting of distance matrix
# - Maintains identical econometric results
.get_withinSS <- function(d, cluster){
  d <- stats::as.dist(d)
  dmat <- as.matrix(d)

  # Handle cluster renumbering if needed
  clusterf <- as.factor(cluster)
  clusterl <- levels(clusterf)
  cn <- length(clusterl)

  if (max(cluster) != cn) {
    warning("cluster renumbered because maximum != number of clusters")
    cluster <- as.integer(clusterf)
  }

  # OPTIMIZED: Pre-compute cluster membership indices
  # split() creates a list of indices for each cluster - computed once
  cluster_indices <- split(seq_along(cluster), cluster)

  # OPTIMIZED: Compute within-cluster SS using vapply (faster than for loop)
  # For each cluster: sum(d[i,j]^2) / cluster_size for all pairs in cluster
  within_ss <- vapply(cluster_indices, function(idx) {
    if (length(idx) <= 1) return(0)
    # Extract submatrix for this cluster and compute sum of squared distances
    di <- dmat[idx, idx, drop = FALSE]
    # Only use lower triangle (as.dist) to avoid double counting
    sum(di[lower.tri(di)]^2) / length(idx)
  }, FUN.VALUE = numeric(1))

  sum(within_ss)
}


# Locate the "elbow" of a within-cluster sum of squares curve.
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# k_values : integer cluster counts on the x-axis (e.g. 1:k.max)
# wss      : total within sum of squares (non-increasing in k)
# Deterministic max-distance-to-chord heuristic: after rescaling both axes to
# [0, 1] (so the choice does not depend on the units of wss), returns the k whose
# point lies farthest from the straight line joining the first and last points.
# Returns NA when the elbow is undefined: fewer than 3 points, non-finite values,
# a constant curve, or a straight-line curve (no interior point bends away from
# the chord, so there is no elbow to mark).
.wss_elbow <- function(k_values, wss){
  n <- length(k_values)
  if(n < 3 || any(!is.finite(wss))) return(NA_integer_)
  rx <- range(k_values); ry <- range(wss)
  if(diff(rx) == 0 || diff(ry) == 0) return(NA_integer_)
  x <- (k_values - rx[1]) / diff(rx)
  y <- (wss - ry[1]) / diff(ry)
  x1 <- x[1]; y1 <- y[1]; xn <- x[n]; yn <- y[n]
  # perpendicular distance from each point to the (first, last) chord
  num <- abs((yn - y1) * x - (xn - x1) * y + xn * y1 - yn * x1)
  den <- sqrt((yn - y1)^2 + (xn - x1)^2)
  d <- num / den
  interior <- 2:(n - 1)                       # endpoints are on the chord (d = 0)
  # A perfectly straight curve leaves every interior point on the chord (d ~ 0);
  # there is no elbow, so don't mark an arbitrary tie.
  if(max(d[interior]) < 1e-8) return(NA_integer_)
  k_values[interior][which.max(d[interior])]
}




# Visualization of the output returned by the function
# NbClust()
# x : an object generated by the function NbClust()
.viz_NbClust <- function(x, print.summary = TRUE,
                         barfill = "steelblue", barcolor = "steelblue")
  {
     best_nc <- x$Best.nc
    # FIX: R 4.0.0+ deprecation - class() can return multiple values for matrices
    # Using inherits() or is.X() functions instead of class() == comparison
    # See: https://github.com/kassambara/factoextra/issues/171
    if(is.numeric(best_nc) && !is.matrix(best_nc)) print(best_nc)
     else if(is.matrix(best_nc)){
    best_nc <- as.data.frame(t(best_nc))
    best_nc$Number_clusters <- as.factor(best_nc$Number_clusters)
    ss <- summary(best_nc$Number_clusters)
    
    # Summary
    if(print.summary){
      cat ("Among all indices: \n===================\n")
      for(i in seq_along(ss)){
        cat("*", ss[i], "proposed ", names(ss)[i], "as the best number of clusters\n" )
      }
      cat("\nConclusion\n=========================\n")
      cat("* According to the majority rule, the best number of clusters is ",
          names(which.max(ss)),  ".\n\n")
    }
    # Fix for Issue #131: ensure numeric ordering of clusters (not alphabetical)
    # When clusters > 9, alphabetical ordering would put "10" before "2"
    cluster_names <- names(ss)
    cluster_order <- order(as.numeric(cluster_names))
    df <- data.frame(
      Number_clusters = factor(cluster_names, levels = cluster_names[cluster_order]),
      freq = ss,
      stringsAsFactors = FALSE
    )
    p <- ggpubr::ggbarplot(df,  x = "Number_clusters", y = "freq", fill = barfill, color = barcolor)+
      labs(x = "Number of clusters k", y = "Frequency among all indices",
           title = paste0("Optimal number of clusters - k = ", names(which.max(ss)) ))
    
    return(p)
  }
}


#  Determines the location of the maximum; see cluster::maxSE.
# +++++++++++++++++++++++++++++++++++++++++++
# f: numeric vector containing the gap statistic
# SE.f : standard error of the gap statistic
# method : character string indicating how the "optimal" number of clusters, k, 
  # is computed from the gap statistics (and their standard deviations), 
  # or more generally how the location k^ of the maximum of f[k] should be determined.
# SE.factor:  Determining the optimal number of clusters, Tibshirani et al. proposed the "1 S.E."-rule.
.maxSE <- function (f, SE.f, method = c("firstSEmax", "Tibs2001SEmax", 
                                        "globalSEmax", "firstmax", "globalmax"), SE.factor = 1) 
{
  method <- match.arg(method)
  stopifnot((K <- length(f)) >= 1, K == length(SE.f), SE.f >= 
              0, SE.factor >= 0)
  fSE <- SE.factor * SE.f
  switch(method, firstmax = {
    decr <- diff(f) <= 0
    if (any(decr)) which.max(decr) else K
  }, globalmax = {
    which.max(f)
  }, Tibs2001SEmax = {
    g.s <- f - fSE
    if (any(mp <- f[-K] >= g.s[-1])) which.max(mp) else K
  }, firstSEmax = {
    decr <- diff(f) <= 0
    nc <- if (any(decr)) which.max(decr) else K
    if (any(mp <- f[seq_len(nc - 1)] >= f[nc] - fSE[nc])) which(mp)[1] else nc
  }, globalSEmax = {
    nc <- which.max(f)
    if (any(mp <- f[seq_len(nc - 1)] >= f[nc] - fSE[nc])) which(mp)[1] else nc
  })
}