Elbow Method

This method involves plotting the within-cluster sum of squares (WSS) as a function of the number of clusters. The WSS is a measure of how well the data points are clustered together. The elbow method works by looking for the point in the WSS curve where the rate of decrease starts to slow down. This point is often considered to be the optimal number of clusters.

library(factoextra)
fviz_nbclust(df_scaled, kmeans, method = "wss") +  
  geom_vline(xintercept = 2, linetype = 2) + # add line for better viz
  labs(subtitle = "Elbow method") # add subtitle

Silhouette Score/Width

The silhouette score is a measure of how well each data point is clustered with its own cluster compared to other clusters. It ranges from -1 to 1, with a value of 1 indicating that the data point is perfectly clustered with its own cluster, 0 indicating the data point is on the border between the two clusters and a value of -1 indicating that the data point is equally well-clustered with two or more clusters. The optimal number of clusters is often the number that maximizes the average silhouette score.

fviz_nbclust(df_scaled, kmeans, method = "silhouette") +  
  labs(subtitle = "Silhouette method")

To compute Silhouette information for a cluster:

library(cluster)
km_res <- kmeans(df_scaled, centers = 2)  # defining cluster

sil <- silhouette(km_res$cluster, dist(df_scaled))
plot(sil)

fviz_silhouette(sil) + theme_minimal()

#   cluster size ave.sil.width
# 1       1   12          0.26
# 2       2    8          0.48

A positive silhouette coefficient indicates that an observation is well-matched to its own cluster.

Gap Statistic Method

The gap statistic helps determine the optimal number of clusters by comparing the observed within-cluster dispersion with a reference distribution. The idea is to select the number of clusters that maximizes the gap statistic. A larger gap indicates a more significant deviation from randomness, suggesting a better-defined clustering structure. By comparing the gaps for different numbers of clusters, we can identify the number of clusters that provides the most distinct and meaningful clustering pattern.

fviz_nbclust(df_scaled, kmeans, method = "gap_stat") +
  labs(subtitle = "Gap statistic method")

Combination of All

But the clever way is to try out different methods and see what majority returns. See here

The parameter package comes to help in this case:

library(parameters)
n_clust <- n_clusters(df_scaled,
                      package = c("easystats", "NbClust", "mclust"),
                      standardize = TRUE
                      )
print(n_clust)

# # Method Agreement Procedure:
# 
# The choice of 2 clusters is supported by 7 (31.82%) methods out of 22 (Elbow, Silhouette, kl, Duda, Pseudot2, Ratkowsky, Mcclain).

plot(n_clust)

Step 5.1: Distance

df_dist_euc <- dist(df_scaled, method = 'euclidian')
# df_dist_max <- dist(df_scaled, method = 'maximum')
# df_dist_man <- dist(df_scaled, method = 'manhattan')
# df_dist_can <- dist(df_scaled, method = 'canberra')
# df_dist_mink <- dist(df_scaled, method = 'minkowski')

# dist(dummified_data, method = "binary")

df <- dplyr::mutate_if(df, is.character, as.factor) # convert character to factor

df_dist_gower <- cluster::daisy(df, metric = "gower") # Dissimilarity Matrix Calculation

Step 5.2: Linkage Criteria

Linkage refers to the way in which clusters are merged together. There are several different linkage criteria that can be used, and the choice of linkage criterion can have a significant impact on the results of the clustering -

• Complete linkage: Maximum distance between two sets
  
• Single linkage: Minimum distance between two sets
  
• Average linkage: Average distance between two sets
  
• Ward’s linkage: Minimize the total within cluster variance

Using the Euclidean distance on the numeric data the following code finds the hierarchical clusters by different linkage methods:

hclust_single <- hclust(df_dist_euc, method = "single")
hclust_complete <- hclust(df_dist_euc, method = "complete")
hclust_average <- hclust(df_dist_euc, method = "average")
hclust_centroid <- hclust(df_dist_euc, method = "centroid")
hclust_ward <- hclust(df_dist_euc, method = "ward.D2")

Get the cluster information:

cutree(hclust_ward, k = 2) # specify number of clusters

# 202 102 112 151 198   4  55  70  98 135   7 183  43 189 140  51 146  25   2 179 
#   1   1   1   1   1   2   2   2   2   1   2   1   2   1   1   2   1   2   2   1

cutree(hclust_ward, h = 0.7) # specify height to cut

# 202 102 112 151 198   4  55  70  98 135   7 183  43 189 140  51 146  25   2 179 
#   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20

Summary by clusters:

library(dplyr)
df_num %>%
  mutate(cluster = cutree(hclust_ward, k = 2)) %>% 
  group_by(cluster) %>%
  summarise(across(everything(), list(Average = mean))) 

# # A tibble: 2 × 12
#   cluster Aroma_Average Flavor_Average Aftertaste_Average Acidity_Average
#     <int>         <dbl>          <dbl>              <dbl>           <dbl>
# 1       1          7.55           7.55               7.37            7.48
# 2       2          8.00           8.03               7.85            7.91
# # ℹ 7 more variables: Body_Average <dbl>, Balance_Average <dbl>,
# #   Overall_Average <dbl>, Total.Cup.Points_Average <dbl>,
# #   Moisture.Percentage_Average <dbl>, Category.One.Defects_Average <dbl>,
# #   Quakers_Average <dbl>

Dendrogram will help to visualize as well as to determine the number of clusters to use or the height where to cut.

Step 5.3: Dendrogram

Pass the hclust object to the plot function:

plot(hclust_ward, hang = -1)
rect.hclust(hclust_ward, k = 2, border = "blue")

Using dendextend package:

# install.packages("dendextend")
library(dendextend)

plot(color_branches(dend = hclust_ward, k = 2)) # specify k, can also use h to specify height

plot(color_branches(dend = hclust_ward, k = 2), 
     # parameters for nodes
     nodePar = list(col = 3:2, 
                    cex = c(1.5, 0.75), 
                    pch =  21:22,
                    bg =  c("light blue", "pink"),
                    lab.cex = 0.75, 
                    lab.col = "tomato"))

# function to get color labels; Change color and cut based on number of clusters
colLab <- function(n, labelColors = c("#036564", "#EB6841"), cut = 2) {
  clusMember = cutree(hclust_ward, cut)
  if (is.leaf(n)) {
    a <- attributes(n)
    labCol <- labelColors[clusMember[which(names(clusMember) == a$label)]]
    attr(n, "nodePar") <- c(a$nodePar, lab.col = labCol)
  }
  n
}

hclust_ward %>% 
  color_branches(k = 2, col = c("#036564", "#EB6841")) %>% 
  dendrapply(colLab) %>%  # applying to all nodes
  plot()  # finally plotting

library("ape")
plot(as.phylo(hclust_ward), type = "unrooted", cex = 0.7, no.margin = T)