Search code examples
rcluster-analysis

factoextra::fviz_gap_stat() versus factoextra::fviz_nbclust(df, method = "gap_stat")

I'm trying to figure out why these two functions from the factoextra package with seemingly similar parameters (e.g. kmeans, gap_stat,k.maxandB`) are yielding different results.

library(cluster)
library(cluster.datasets)
library(tidyverse)
library(factoextra)

# load data and scale it
data("all.mammals.milk.1956")
mammals <- all.mammals.milk.1956 %>% select(-name)
mammals_scaled <- scale(mammals)

The first method uses factoextra::clusGap() and factoextra::fviz_gap_stat()

gap_stat <- clusGap(mammals_scaled, FUN = kmeans, K.max = 24, B = 50)

fviz_gap_stat(gap_stat) + theme_minimal() + ggtitle("fviz_gap_stat: Gap Statistic")

enter image description here

The second method uses factoextra::fviz_nbclust() which 

fviz_nbclust(mammals_scaled, kmeans, method = "gap_stat", k.max = 24, nboot = 50) + theme_minimal() + ggtitle("fviz_nbClust_gap_stat: Gap Statistic")

enter image description here

I thought it could possibly be the nstart option from clusGap() but when I use the jimhester/lookup to read the source code of fviz_nbclust() with the following code I couldn't find what the issue was: 

devtools::install_github("jimhester/lookup")
lookup::lookup(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, ...) 
{
        set.seed(123)
        if (k.max < 2) 
                stop("k.max must bet > = 2")
        method = match.arg(method)
        if (!inherits(x, c("data.frame", "matrix")) & !("Best.nc" %in% 
                names(x))) 
                stop("x should be an object of class matrix/data.frame or ", 
                        "an object created by the function NbClust() [NbClust package].")
        if (inherits(x, "list") & "Best.nc" %in% names(x)) {
                best_nc <- x$Best.nc
                if (class(best_nc) == "numeric") 
                        print(best_nc)
                else if (class(best_nc) == "matrix") 
                        .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 (method %in% c("silhouette", "wss")) {
                if (is.data.frame(x)) 
                        x <- as.matrix(x)
                if (is.null(diss)) 
                        diss <- stats::dist(x)
                v <- rep(0, k.max)
                if (method == "silhouette") {
                        for (i in 2:k.max) {
                                clust <- FUNcluster(x, i, ...)
                                v[i] <- .get_ave_sil_width(diss, clust$cluster)
                        }
                }
                else if (method == "wss") {
                        for (i in 1:k.max) {
                                clust <- FUNcluster(x, i, ...)
                                v[i] <- .get_withinSS(diss, clust$cluster)
                        }
                }
                df <- data.frame(clusters = as.factor(1:k.max), y = v)
                ylab <- "Total Within Sum of Square"
                if (method == "silhouette") 
                        ylab <- "Average silhouette width"
                p <- ggpubr::ggline(df, x = "clusters", y = "y", group = 1, 
                        color = linecolor, ylab = ylab, xlab = "Number of clusters k", 
                        main = "Optimal number of clusters")
                if (method == "silhouette") 
                        p <- p + geom_vline(xintercept = which.max(v), linetype = 2, 
                                color = linecolor)
                return(p)
        }
        else if (method == "gap_stat") {
                extra_args <- list(...)
                gap_stat <- cluster::clusGap(x, FUNcluster, K.max = k.max, 
                        B = nboot, verbose = verbose, ...)
                if (!is.null(extra_args$maxSE)) 
                        maxSE <- extra_args$maxSE
                else maxSE <- list(method = "firstSEmax", SE.factor = 1)
                p <- fviz_gap_stat(gap_stat, linecolor = linecolor, 
                        maxSE = maxSE)
                return(p)
        }
}
Solution

  • The difference is right at the beginning of the fviz_nbclust function. In line 6 the random seed is set:
    set.seed(123)

    Because the kmeans algorithm uses a random start the results can be different in repeated runs. For example, I used your data with two different random seeds to arrive at slightly different results. 

    set.seed(123)  
gap_stat <- cluster::clusGap(mammals_scaled, FUN = kmeans, K.max = 24, B = 50)   
fviz_gap_stat(gap_stat) + theme_minimal() + ggtitle("fviz_gap_stat: Gap Statistic")

    seed 123 gap stat

    set.seed(42)  
gap_stat <- cluster::clusGap(mammals_scaled, FUN = kmeans, K.max = 24, B = 50)
fviz_gap_stat(gap_stat) + theme_minimal() + ggtitle("fviz_gap_stat: Gap Statistic")

    seed 42 gap stat

    I'm not entirely sure why the seed 123 results are not the same but I think it is connected to the fact that in my code it is executed right above the clusGap function and in Fviz_nbclust several other commands are evaluated in between.