Determining and Visualizing the Optimal Number of Clusters

Partitioning methods, such as k-means clustering require the users to specify the number of clusters to be generated.

fviz_nbclust(): Determines and visualizes the optimal number of clusters using different methods: within-cluster sum of squares, average silhouette, and the gap statistic. Silhouette values are evaluated only for k = 2, ..., k.max because average silhouette width is undefined for a one-cluster partition. The "silhouette" and "gap_stat" plots mark their optimum with a dashed guide line; set mark_optimal = TRUE to also mark the elbow on the "wss" plot, or mark_optimal = FALSE to omit the guide line for every method.
fviz_gap_stat(): Visualizes the gap statistic generated by the function clusGap() [in cluster package]. The optimal number of clusters is specified using the maxSE method with method = "firstSEmax".

For method = "wss", factoextra computes the k = 1 baseline internally so helper functions such as hcut() and hkmeans() can keep rejecting direct k = 1 inputs.

Usage

fviz_nbclust(
  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
)

fviz_gap_stat(
  gap_stat,
  linecolor = "steelblue",
  maxSE = list(method = "firstSEmax", SE.factor = 1),
  mark_optimal = NULL
)

Arguments

x

numeric matrix or data frame. In the function fviz_nbclust(), x can be the results of the function NbClust(). For method = "silhouette" or "wss", x may also be a precomputed dissimilarity (an object of class "dist"); it is then passed directly to FUNcluster, which must be a dissimilarity-capable method (e.g. cluster::pam or factoextra::hcut). method = "gap_stat" still needs the raw data.

FUNcluster

a partitioning function which accepts as first argument a (data) matrix like x, second argument, say k >= 2, the number of clusters desired, and returns a list with a component named cluster which contains the grouping of observations. Allowed values include: kmeans, cluster::pam, cluster::clara, cluster::fanny, hcut, etc. In method = "wss" mode, fviz_nbclust() computes the k = 1 baseline internally instead of calling FUNcluster(x, 1, ...). This argument is not required when x is an output of the function NbClust::NbClust().

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).

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"

k.max

the maximum number of clusters to consider, must be at least two.

nboot

integer, number of Monte Carlo ("bootstrap") samples. Used only for determining the number of clusters using the gap statistic.

verbose

logical value. If TRUE, progress information is printed.

barfill, barcolor

fill color and outline color for bars

linecolor

color for lines

print.summary

logical value. If TRUE, the optimal number of clusters is printed in fviz_nbclust().

...

optionally further arguments: arguments for FUNcluster() in "wss"/"silhouette" modes; arguments for clusGap() in "gap_stat" mode. A maxSE list can also be supplied in "gap_stat" mode and is forwarded to fviz_gap_stat().

mark_optimal

logical, or NULL (default). NULL keeps each method's standard marker: a dashed guide line at the estimated optimal number of clusters is drawn for method = "silhouette" (maximum average silhouette width) and method = "gap_stat" (the maxSE location), and omitted for method = "wss". Set TRUE to also mark the "wss" elbow (a maximum-distance heuristic that returns a candidate even when the data has no clear cluster structure; see Details), or FALSE to omit the guide line for every method.

gap_stat

an object of class "clusGap" returned by the function clusGap() [in cluster package]

maxSE

a list containing the parameters method and SE.factor used by maxSE to locate the gap statistic optimum. The default is list(method = "firstSEmax", SE.factor = 1). Allowed methods include:

"globalmax": simply corresponds to the global maximum, i.e., is which.max(gap)
"firstmax": gives the location of the first local maximum
"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.
"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.
see maxSE for more options

Value

Both fviz_nbclust() and fviz_gap_stat() return a ggplot2 object.

Details

When mark_optimal = TRUE, the "wss" elbow is located with a deterministic chord-distance heuristic: both axes are rescaled to the unit interval [0, 1] (so the choice does not depend on the magnitude of the within-cluster sum of squares) and the marked 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 Satopää 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 "silhouette" and "gap_stat" guide lines instead mark defined optima (the maximum average silhouette width and the maxSE location).

Method selection for gap statistic

The default "firstSEmax" method returns the first value within SE.factor standard errors of the first local maximum. Other maxSE rules can be selected explicitly through maxSE.

References

Satopää, V., Albrecht, J., Irwin, D. & Raghavan, B. (2011). Finding a "Kneedle" in a Haystack: Detecting Knee Points in System Behavior. 2011 31st International Conference on Distributed Computing Systems Workshops, 166-171. doi:10.1109/ICDCSW.2011.20

Author

Alboukadel Kassambara alboukadel.kassambara@gmail.com

Examples

set.seed(123)

# Data preparation
# +++++++++++++++
data("iris")
head(iris)
#>   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
#> 1          5.1         3.5          1.4         0.2  setosa
#> 2          4.9         3.0          1.4         0.2  setosa
#> 3          4.7         3.2          1.3         0.2  setosa
#> 4          4.6         3.1          1.5         0.2  setosa
#> 5          5.0         3.6          1.4         0.2  setosa
#> 6          5.4         3.9          1.7         0.4  setosa
# 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)
#> Warning: package ‘cluster’ was built under R version 4.5.2
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")
#> Clustering Gap statistic ["clusGap"] from call:
#> clusGap(x = iris.scaled, FUNcluster = kmeans, K.max = 10, B = 10, nstart = 25)
#> B=10 simulated reference sets, k = 1..10; spaceH0="scaledPCA"
#>  --> Number of clusters (method 'firstSEmax', SE.factor=1): 2
#>           logW   E.logW       gap      SE.sim
#>  [1,] 4.534565 4.745781 0.2112157 0.024271468
#>  [2,] 4.021316 4.481045 0.4597287 0.023247363
#>  [3,] 3.806577 4.287108 0.4805310 0.022005009
#>  [4,] 3.699263 4.138042 0.4387785 0.022976513
#>  [5,] 3.589284 4.046911 0.4576270 0.018998839
#>  [6,] 3.522810 3.971789 0.4489795 0.013118375
#>  [7,] 3.448288 3.906691 0.4584031 0.011975368
#>  [8,] 3.379870 3.851128 0.4712584 0.008759733
#>  [9,] 3.310088 3.801559 0.4914709 0.007393207
#> [10,] 3.278659 3.757545 0.4788863 0.009164325
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)