I am using K-means to cluster my data and was looking for a way to suggest an "optimal" cluster number. Gap statistics seems to be a common way to find a good cluster number.

For some reason it returns 1 as optimal cluster number, but when I look at the data it's obvious that there are 2 clusters:


This is how I call gap in R:

gap <- clusGap(data, FUN=kmeans, K.max=10, B=500)
with(gap, maxSE(Tab[,"gap"], Tab[,"SE.sim"], method="firstSEmax"))

The result set:

> Number of clusters (method 'firstSEmax', SE.factor=1): 1
          logW   E.logW           gap    SE.sim
[1,]  5.185578 5.085414 -0.1001632148 0.1102734
[2,]  4.438812 4.342562 -0.0962498606 0.1141643
[3,]  3.924028 3.884438 -0.0395891064 0.1231152
[4,]  3.564816 3.563931 -0.0008853886 0.1387907
[5,]  3.356504 3.327964 -0.0285393917 0.1486991
[6,]  3.245393 3.119016 -0.1263766015 0.1544081
[7,]  3.015978 2.914607 -0.1013708665 0.1815997
[8,]  2.812211 2.734495 -0.0777154881 0.1741944
[9,]  2.672545 2.561590 -0.1109558011 0.1775476
[10,] 2.656857 2.403220 -0.2536369287 0.1945162

Am I doing something wrong or does someone know a better way to get a good cluster number?


3 Answers 3


Clustering depends on scale, among other things. For discussions of this issue see (inter alia) When should you center and standardize data? and PCA on covariance or correlation?.

Here are your data drawn with a 1:1 aspect ratio, revealing how much the scales of the two variables differ:

Figure 1

To its right, the plot of the gap stats shows the statistics by number of clusters ($k$) with standard errors drawn with vertical segments and the optimal value of $k$ marked with a vertical dashed blue line. According to the clusGap help,

The default method "firstSEmax" looks for the smallest $k$ such that its value $f(k)$ is not more than 1 standard error away from the first local maximum.

Other methods behave similarly. This criterion does not cause any of the gap statistics to stand out, resulting in an estimate of $k=1$.

Choice of scale depends on the application, but a reasonable default starting point is a measure of dispersion of the data, such as the MAD or standard deviation. This plot repeats the analysis after recentering to zero and rescaling to make a unit standard deviation for each component $a$ and $b$:

Figure 2

The $k=2$ K-means solution is indicated by varying symbol type and color in the scatterplot of the data at left. Among the set $k\in\{1,2,3,4,5\}$, $k=2$ is clearly favored in the gap statistics plot at right: it is the first local maximum and the stats for smaller $k$ (that is, $k=1$) are significantly lower. Larger values of $k$ are likely overfit for such a small dataset, and none are significantly better than $k=2$. They are shown here only to illustrate the general method.

Here is R code to produce these figures. The data approximately match those shown in the question.

xy <- matrix(c(29,391, 31,402, 31,380, 32.5,391, 32.5,360, 33,382, 33,371,
        34,405, 34,400, 34.5,404, 36,343, 36,320, 36,303, 37,344,
        38,358, 38,356, 38,351, 39,318, 40,322, 40, 341), ncol=2, byrow=TRUE)
colnames(xy) <- c("a", "b")
title <- "Raw data"
for (i in 1:2) {
  # Estimate optimal cluster count and perform K-means with it.
  gap <- clusGap(xy, kmeans, K.max=10, B=500)
  k <- maxSE(gap$Tab[, "gap"], gap$Tab[, "SE.sim"], method="Tibs2001SEmax")
  fit <- kmeans(xy, k)
  # Plot the results.
  pch <- ifelse(fit$cluster==1,24,16); col <- ifelse(fit$cluster==1,"Red", "Black")
  plot(xy, asp=1, main=title, pch=pch, col=col)
  plot(gap, main=paste("Gap stats,", title))
  abline(v=k, lty=3, lwd=2, col="Blue")
  # Prepare for the next step.
  xy <- apply(xy, 2, scale)
  title <- "Standardized data"
  • $\begingroup$ alright thanks for the explanation. Btw: Do you know any other cluster metric like gap statistics? I found some, but I don't know which one are usually used with k-means? $\endgroup$
    – MikeHuber
    Mar 6, 2015 at 23:22
  • $\begingroup$ +1. Very nice demonstration, and it's impressive that you seem to have digitized OP's figure to get the same data. $\endgroup$
    – amoeba
    Mar 6, 2015 at 23:37
  • 3
    $\begingroup$ @amoeba I eyeballed the scatterplot and typed in the coordinates exactly as you see here (that is, the digits involved were my own fingers :-)). Sometimes the simplest approach is efficient. $\endgroup$
    – whuber
    Mar 7, 2015 at 1:55
  • $\begingroup$ Can the gap statistic be used for finding the number of clusters in a single 1-d array of numeric values? $\endgroup$ Apr 27, 2016 at 9:41
  • 1
    $\begingroup$ @majom Negative standardized values don't mean what you seem to think they do. As the first figure shows, they still correspond to positive gaps. If you don't like negative numbers, add a sufficiently large positive constant to all the standardized values: that will not change the result of this procedure in any way. $\endgroup$
    – whuber
    Jun 10, 2021 at 20:16

I think you do not understand anything wrong in your use of the GAP statistic. I believe though you are partially mislead by the scale of the data in the visualization. You see two clusters but actually the x direction is rather small compared to the y direction. Based on that you would expect two enlonged clusters. Nevertheless it looks like your one mode of variance dominates the other. As the GAP statistics assumes a null model with a single component ($K=1$) and then tries to reject this model for an alternative one with $K>1$; what you observe is the inability to reject the null. Please note that the inability to reject the null hypothesis does not make it true. The methodological paper describing the GAP statistic it is available online if you want to check the technical particulars more.

I run your model using a Gaussian Mixture Model (GMM - a generalization of $k$-means, see this thread for more on that matter). True enough in that case too the GAP statistic suggested a single cluster. The BIC also suggested a single cluster. AIC suggests 4 clusters (!), this being a clear sign we start to overfit. The sample used is not extremely big; you have 21 points where one mode of variance dominates over the other. It is a bit of a stretch to have two 2-D clusters (ie, fit two 2-D means and two $2 \times 2$ covariances matrices) with just 21 2-D points. :) (In the case of $k$-means your covariance matrix is more structured (you do not look at covariances) but I would not focus on that matter here.)

EDIT: Just for completeness: @whuber showed that two clusters would appear as optimal in $k$-means if one standardised his data; the GAP criterion applied on the GMM fit will also give $K=2$ as the optimal number of clusters if one standardises the data.

  • 1
    $\begingroup$ +1 You saw the potential problem by carefully reading the plot: well done! The link to Hastie's paper is welcome support to your answer, too. $\endgroup$
    – whuber
    Mar 6, 2015 at 23:06
  • $\begingroup$ @whuber: We had this discussion about scales, didn't we? :) $\endgroup$
    – usεr11852
    Mar 6, 2015 at 23:08
  • $\begingroup$ It was such a different context I didn't make the connection ... . $\endgroup$
    – whuber
    Mar 6, 2015 at 23:11
  • $\begingroup$ It was a different context indeed; I mentioned it to you just because it was "scales" there, and "scales" here. $\endgroup$
    – usεr11852
    Mar 6, 2015 at 23:14

I had the same problem as the original poster. R documentation currently says that original and default setting of d.power = 1 was incorrect and should be replaced by d.power: "The default, d.power = 1, corresponds to the “historical” R implementation, whereas d.power = 2 corresponds to what Tibshirani et al had proposed. This was found by Juan Gonzalez, in 2016-02."

Consequently, changing d.power = 2 solved the problem for me.



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.