R/eclust.R

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

#' @include utilities.R cluster_utilities.R dist.R fviz_cluster.R fviz_dend.R
NULL
#' Visual enhancement of clustering analysis
#' 
#' @description Provides a convenient workflow for clustering analyses and
#'   ggplot2-based data visualization. When \code{k = NULL}, the gap statistic
#'   selects the number of clusters. Hierarchical backends may validly return
#'   \code{k = 1}; in that case \code{eclust()} returns a one-cluster result
#'   without silhouette information. Read more:
#'  \href{https://www.datanovia.com/en/blog/cluster-analysis-in-r-simplified-and-enhanced/}{Visual enhancement of clustering analysis}.
#' @param x numeric vector, data matrix or data frame. For hierarchical
#'   clustering (\code{FUNcluster} = "hclust", "agnes" or "diana"), a
#'   precomputed dissimilarity matrix (an object of class \code{"dist"}) may be
#'   supplied directly; in that case \code{hc_metric} is ignored and \code{k}
#'   must be specified. This allows custom distances such as Bray-Curtis
#'   (e.g. \code{vegan::vegdist(df, "bray")}).
#' @param FUNcluster a clustering function including "kmeans", "pam", "clara",
#'   "fanny", "hkmeans", "hclust", "agnes" and "diana". Abbreviation is allowed.
#' @param k the number of clusters to be generated. If NULL, the gap statistic
#'   is used to estimate the appropriate number of clusters. For hierarchical
#'   clustering, this automatic selection may return \code{k = 1}. In the case
#'   of kmeans, \code{k} can be either the number of clusters, or a set of
#'   initial (distinct) cluster centers.
#' @param k.max the maximum number of clusters to consider, must be at least 
#'   two.
#' @param stand logical value; default is FALSE. If TRUE, then the data will be
#'   standardized using the function \code{scale()}. Measurements are
#'   standardized for each variable (column), by subtracting the variable's
#'   mean value and dividing by the variable's standard deviation. If scaling
#'   produces \code{NA} values, \code{eclust()} stops with a package-level
#'   error.
#' @param graph logical value. If TRUE, cluster plot is displayed.
#' @param hc_metric character string specifying the metric to be used for 
#'   calculating dissimilarities between observations. Allowed values are those 
#'   accepted by the function dist() [including "euclidean", "manhattan", 
#'   "maximum", "canberra", "binary", "minkowski"] and correlation based 
#'   distance measures ["pearson", "spearman" or "kendall"]. Used only when
#'   FUNcluster is a hierarchical clustering function such as one of "hclust",
#'   "agnes" or "diana". Ignored when \code{x} is already a \code{"dist"} object.
#' @param hc_method the agglomeration method to be used (?hclust): "ward.D", 
#'   "ward.D2", "single", "complete", "average", ...
#' @param gap_maxSE a list containing the parameters (method and SE.factor) for 
#'   determining the location of the maximum of the gap statistic (Read the 
#'   documentation ?cluster::maxSE).
#' @param nboot integer, number of Monte Carlo ("bootstrap") samples. Used only 
#'   for determining the number of clusters using gap statistic.
#' @param verbose logical value. If TRUE, the result of progress is printed.
#' @param seed integer used for seeding the random number generator.
#' @param ... other arguments to be passed to FUNcluster.
#' @return Returns an object of class "eclust" containing the result of the 
#'   standard function used (e.g., kmeans, pam, hclust, agnes, diana, etc.).
#'   
#'   It also includes: \itemize{ \item cluster: the cluster assignment of
#'   observations after cutting the tree \item nbclust: the number of clusters
#'   \item silinfo: the silhouette information of observations, when available
#'   for solutions with at least two clusters, including $widths (silhouette
#'   width values of each observation), $clus.avg.widths (average silhouette
#'   width of each cluster) and $avg.width (average width of all clusters)
#'   \item size: the size of clusters \item data: a matrix containing the
#'   original or the standardized data (if stand = TRUE) } The "eclust" class
#'   has method for fviz_silhouette(), fviz_dend(), fviz_cluster().
#' @seealso \code{\link{fviz_silhouette}}, \code{\link{fviz_dend}}, 
#'   \code{\link{fviz_cluster}}
#' @author Alboukadel Kassambara \email{alboukadel.kassambara@@gmail.com}
#'   
#' @examples 
#' # Load and scale data
#' data("USArrests")
#' df <- scale(USArrests)
#' 
#' # Enhanced k-means clustering
#' # nboot >= 500 is recommended
#' res.km <- eclust(df, "kmeans", nboot = 2)
#' # Silhouette plot
#' fviz_silhouette(res.km)
#' # Optimal number of clusters using gap statistics
#' res.km$nbclust
#' # Print result
#'  res.km
#'  
#' \dontrun{
#' # Enhanced hierarchical clustering
#' res.hc <- eclust(df, "hclust", nboot = 2) # compute hclust
#' fviz_dend(res.hc) # dendrogram
#' if (res.hc$nbclust > 1) fviz_silhouette(res.hc) # silhouette plot
#' }
#'  
#' @name eclust
#' @rdname eclust
#' @export
eclust <- function(x, FUNcluster = c("kmeans", "pam", "clara", "fanny", "hclust", "agnes", "diana", "hkmeans"),
                   k = NULL, k.max = 10, stand = FALSE,
                   graph = TRUE,
                   hc_metric = "euclidean", hc_method = "ward.D2",
                   gap_maxSE = list(method = "firstSEmax", SE.factor = 1), 
                   nboot = 100, verbose = interactive(),
                   seed = 123,  ...)
  {
  has_seed <- exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)
  if(has_seed) old_seed <- get(".Random.seed", envir = .GlobalEnv)
  on.exit({
    if(!has_seed && exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)){
      rm(".Random.seed", envir = .GlobalEnv)
    } else if(has_seed){
      assign(".Random.seed", old_seed, envir = .GlobalEnv)
    }
  }, add = TRUE)
  set.seed(seed)
  data <- x
  # A precomputed distance matrix (class "dist") may be passed as x for
  # hierarchical clustering (mirrors hcut()); hc_metric is then ignored (#182).
  is_dist <- inherits(x, "dist")
  if(stand) {
    if(is_dist)
      stop("'stand = TRUE' is not supported when 'x' is a distance matrix (class 'dist').")
    x <- scale(x)
    if(anyNA(x))
      stop("Scaling produced NA values. Check for constant columns or non-finite values.")
  }
  # Define the type of clustering
  FUNcluster <- match.arg(FUNcluster)
  if(is_dist && FUNcluster %in% c("kmeans", "pam", "clara", "fanny", "hkmeans"))
    stop("A distance matrix (class 'dist') is supported only for hierarchical clustering ",
         "(FUNcluster = 'hclust', 'agnes', or 'diana'). For '", FUNcluster,
         "', supply the raw data, or use hcut() for precomputed distances.")
  fun_clust <- switch(FUNcluster,
                 kmeans = stats::kmeans,
                 pam = cluster::pam,
                 clara = cluster::clara,
                 fanny = cluster::fanny,
                 hkmeans = hkmeans,
                 diana = hcut,
                 agnes = hcut,
                 hclust = hcut
                 )
  
  if(!inherits(data, c("matrix", "data.frame")) ) graph = FALSE
  else if(ncol(data)< 2) graph = FALSE
      
  gap_stat <- NULL
  auto_k <- is.null(k)
  # Partitioning clustering
  # ++++++++++++++++++++++++++++++
  clust <- list()
  if(FUNcluster %in% c("kmeans", "pam", "clara", "fanny", "hkmeans")){
    # Number of cluster
    if(is.null(k)) {
      gap <- .gap_stat(x, fun_clust, k.max = k.max, nboot = nboot,
                       gap_maxSE = gap_maxSE, verbose = verbose, ...)
      k <- gap$k
      gap_stat <- gap$stat
    }
    clust <- fun_clust(x, k, ...)
    
    if(inherits(k, c("matrix", "data.frame"))) k <- nrow(k) # cluster centers are provided as k
    
    # Plot
    if(graph) {
      clust$clust_plot <- fviz_cluster(clust, x)
      print(clust$clust_plot + labs(title = paste0(toupper(FUNcluster), " Clustering")))
    }
    if(k > 1) clust$silinfo <-.get_silinfo(clust$cluster, stats::dist(x))
  }
  
  
  # Hierarchical clustering
  # ++++++++++++++++++++++++++++++++
  else if(FUNcluster %in% c("hclust", "agnes", "diana")){

    # Use a precomputed distance directly; otherwise compute it from the data.
    res.dist <- if(is_dist) x else get_dist(x, method = hc_metric)
    if(auto_k && is_dist)
      stop("When 'x' is a distance matrix, the number of clusters 'k' must be specified: ",
           "the gap statistic requires the original data and cannot be computed from a 'dist'.")
    # Number of cluster
    if(auto_k) {
      gap <- .gap_stat(x, fun_clust, k.max = k.max, nboot = nboot,
                       gap_maxSE = gap_maxSE, verbose = verbose, diss = res.dist)
      k <- gap$k
      gap_stat <- gap$stat
    }
    if(auto_k && k == 1) {
      res.hc <- .single_cluster_hcut(res.dist, hc_func = FUNcluster, hc_method = hc_method)
    } else {
      res.hc <- hcut(res.dist, k, hc_func = FUNcluster, hc_method = hc_method )
    }
    clust <- res.hc
    if(graph) {
      if(k > 1) fviz_dend(clust, k)
      else fviz_dend(clust)
    }
  }
  
  clust$nbclust <- k
  clust$data <- x
  clust$gap_stat <- gap_stat
  class(clust) <- c(class(clust), "eclust")
  clust
}


# Compute the hierarchical tree for the requested backend
.compute_hc_tree <- function(diss, hc_func, hc_method){
  if(hc_func == "hclust") stats::hclust(diss, method = hc_method)
  else if(hc_func == "agnes") {
    if(hc_method %in% c("ward.D", "ward.D2")) hc_method <- "ward"
    cluster::agnes(diss, method = hc_method)
  }
  else if(hc_func == "diana") cluster::diana(diss)
  else stop("Don't support the function ", hc_func)
}

# Build a one-cluster hierarchical result without relaxing direct hcut validation.
.single_cluster_hcut <- function(diss, hc_func, hc_method){
  n_obs <- attr(diss, "Size")
  if(is.null(n_obs) || !is.numeric(n_obs))
    stop("Unable to determine number of observations from distance data")

  hc <- .compute_hc_tree(diss, hc_func = hc_func, hc_method = hc_method)
  hc$cluster <- rep(1L, n_obs)
  cluster_labels <- attr(diss, "Labels")
  if(!is.null(cluster_labels) && length(cluster_labels) == n_obs)
    names(hc$cluster) <- cluster_labels
  hc$nbclust <- 1L
  hc$size <- as.vector(table(hc$cluster))
  hc$data <- diss
  class(hc) <- c(class(hc), "hcut")
  hc
}


# Compute gap stat and get k
.gap_stat <- function(x, fun_clust, k.max = 10, nboot = 100,
                   gap_maxSE = list(method = "firstmax", SE.factor = 1),
                   verbose = interactive(), ...)
  {
  gap_stat <- cluster::clusGap(x, fun_clust, K.max = k.max,  B = nboot, 
                               verbose = verbose, ...)
  gap <- gap_stat$Tab[, "gap"]
  se <- gap_stat$Tab[, "SE.sim"]
  k <- .maxSE(gap, se, method = gap_maxSE$method, SE.factor = gap_maxSE$SE.factor)
  list(stat = gap_stat, k = k)
}