Determining the optimal number of clusters in your dataset

1. ELBOW METHOD

In short, the elbow method maps the within-cluster sum of squares onto the number of possible clusters. As a rule of thumb, you pick the number for which you see a significant decrease in the within-cluster dissimilarity, or so called ‘elbow’. You can find more details here or here.

Let’s see what happens when we apply the Elbow Method to our dataset:

set.seed(13)

wss <- (nrow(scaled_wine2)-1)*sum(apply(scaled_wine2,2,var))
for (i in 2:15) wss[i] <- sum(kmeans(scaled_wine2,
                                     centers=i)$withinss)
plot(1:15, wss, type="b", xlab="Number of Clusters",
     ylab="Within groups sum of squares")

In this example, we see an ‘elbow’ at 3 clusters. By adding more clusters than that we get relatively smaller gains. So… 3 clusters!

2. SILHOUETTE ANALYSIS

The silhouette plots display a measure of how close each point in one cluster is to points in the neighboring clusters. This measure ranges from -1 to 1, where 1 means that points are very close to their own cluster and far from other clusters, whereas -1 indicates that points are close to the neighbouring clusters. You’ll find more details about silhouette analysis here and here.

When we run the silhouette analysis on our wine dataset, we get the following graphs with indication on the optimal number of clusters:

library(fpc)

set.seed(13)

pamk.best2 <- pamk(scaled_wine2)
cat("number of clusters estimated by optimum average silhouette width:", pamk.best2$nc, "\n")

## number of clusters estimated by optimum average silhouette width: 3

plot(pam(scaled_wine2, pamk.best2$nc))

Thus, it’s 3 clusters again.

3. GAP STATISTICS

Gap statistic is a goodness of clustering measure, where for each hypothetical number of clusters k, it compares two functions: log of within-cluster sum of squares (wss) with its expectation under the null reference distribution of the data. In essence, it standardizes wss. It chooses the value where the log(wss) is the farthest below the reference curve, ergo the gap statistic is maximum. Pretty neat!

If you want to know more, check out this presentation for a very good explanation and visualization of what gap statistic is.

Also, for even more details, see the original paper where the gap statistic was first described, it’s like a good novel!

In the case of our wine dataset, gap analysis picks (again!) 3 cluster as an optimum:

library(cluster)

set.seed(13)

clusGap(scaled_wine2, kmeans, 10, B = 100, verbose = interactive())

## Clustering Gap statistic ["clusGap"] from call:
## clusGap(x = scaled_wine2, FUNcluster = kmeans, K.max = 10, B = 100,     verbose = interactive())
## B=100 simulated reference sets, k = 1..10; spaceH0="scaledPCA"
##  --> Number of clusters (method 'firstSEmax', SE.factor=1): 3
##           logW   E.logW       gap     SE.sim
##  [1,] 5.377557 5.864072 0.4865153 0.01261541
##  [2,] 5.203502 5.757547 0.5540456 0.01356018
##  [3,] 5.066921 5.697040 0.6301193 0.01404664
##  [4,] 5.033404 5.651970 0.6185658 0.01410140
##  [5,] 4.992739 5.616114 0.6233751 0.01342023
##  [6,] 4.955442 5.584673 0.6292313 0.01304916
##  [7,] 4.933420 5.556602 0.6231820 0.01250233
##  [8,] 4.908585 5.531523 0.6229379 0.01299415
##  [9,] 4.884336 5.508464 0.6241281 0.01209969
## [10,] 4.858595 5.488182 0.6295870 0.01336910

4. HIERARCHICAL CLUSTERING

Hierarchical clustering is usually used to better understand the structure and relationships in your data and based on them you decide what number of clusters seems appropriate for your purpose. The way the algorithm works is very well described here and here.

How to choose the optimal number of clusters based on the output of this analysis, the dendogram? As a rule of thumb, look for the clusters with the longest ‘branches’, the shorter they are, the more similar they are to following ‘twigs’ and ‘leaves’. But keep in mind that, as always, the optimal number will also depend on the context, expert knowledge, application, etc.

In our case, the solution is rather blurred:

set.seed(13)

wine_dist2 <- dist(as.matrix(scaled_wine2))   # find distance matrix
plot(hclust(wine_dist2))

Personally, I would go for something between 2 and 4 clusters. What do you think? :-)

5. CALINSKY CRITERION

yet another method evaluating optimal number of clusters based on within- and between-cluster distances. For more details, see here.

In the example below we set the number of clusters to test between 1 and 10 and we iterated the test 1000 times.

library(vegan)

set.seed(13)

cal_fit2 <- cascadeKM(scaled_wine2, 1, 10, iter = 1000)
plot(cal_fit2, sortg = TRUE, grpmts.plot = TRUE)

calinski.best2 <- as.numeric(which.max(cal_fit2$results[2,]))
cat("Calinski criterion optimal number of clusters:", calinski.best2, "\n")

## Calinski criterion optimal number of clusters: 3

Again, 3 clusters!

6. BAYESIAN INFORMATION CRITERION FOR EXPECTATION - MAXIMIZATION

This time in order to determine the optimal number of cluters we simultaneously run and compare several probabilistic models and pick the best one based on their BIC (Bayesian Information Criterion) score. For more details, make sure to read this and this.

In the example below, we set the search of optimal number of clusters to be between 1 and 20.

library(mclust)

set.seed(13)

d_clust2 <- Mclust(as.matrix(scaled_wine2), G=1:20)
m.best2 <- dim(d_clust2$z)[2]

cat("model-based optimal number of clusters:", m.best2, "\n")

## model-based optimal number of clusters: 3

plot(d_clust2)

And this time, the optimal number of clusters seems to be 4.

7. AFFINITY PROPAGATION CLUSTERING

I’ll be quite honest with you: I’ve struggled to get my head around what affinity propagation is about! I understood that it clusters data by identifying a subset of representative examples, so called exemplars. And it does so by considering all data points as potential exemplars, but the rest is magic to me..

\[and\] exchanging real-valued messages between data points until a high-quality set of exemplars and corresponding clusters gradually emerges.

(according to Delbert Dueck)

So if you have a better explanation or better materials.. send them over!

In our particular case, this algorithm identified

library(apcluster)

set.seed(13)

d.apclus2 <- apcluster(negDistMat(r=2), scaled_wine2)
cat("affinity propogation optimal number of clusters:", length(d.apclus2@clusters), "\n")

## affinity propogation optimal number of clusters: 15

heatmap(d.apclus2)

plot(d.apclus2, scaled_wine2)