Introduction

Large amounts of data are collected every day from satellite images, bio-medical, security, marketing, web search, geo-spatial or other automatic equipment. Mining knowledge from these big data far exceeds human’s abilities.

Clustering is one of the important data mining methods for discovering knowledge in multidimensional data. The goal of clustering is to identify pattern or groups of similar objects within a data set of interest.

In the litterature, it is referred as “pattern recognition” or “unsupervised machine learning” - “unsupervised” because we are not guided by a priori ideas of which variables or samples belong in which clusters. “Learning” because the machine algorithm “learns” how to cluster.

Cluster analysis is popular in many fields, including:

• In cancer research for classifying patients into subgroups according their gene expression profile. This can be useful for identifying the molecular profile of patients with good or bad prognostic, as well as for understanding the disease.

• In marketing for market segmentation by identifying subgroups of customers with similar profiles and who might be receptive to a particular form of advertising.

• In City-planning for identifying groups of houses according to their type, value and location.

Preview of the first 38 pages of the book: Practical Guide to Cluster Analysis in R (preview).

Download the ebook through payhip:

Order a physical copy from amazon:

Key features of this book

Although there are several good books on unsupervised machine learning/clustering and related topics, we felt that many of them are either too high-level, theoretical or too advanced. Our goal was to write a practical guide to cluster analysis, elegant visualization and interpretation.

The main parts of the book include:

• distance measures,
• partitioning clustering,
• hierarchical clustering,
• cluster validation methods, as well as,
• advanced clustering methods such as fuzzy clustering, density-based clustering and model-based clustering.

The book presents the basic principles of these tasks and provide many examples in R. This book offers solid guidance in data mining for students and researchers.

Key features:

• Covers clustering algorithm and implementation
• Key mathematical concepts are presented
• Short, self-contained chapters with practical examples. This means that, you don’t need to read the different chapters in sequence.

How this book is organized?

This book contains 5 parts. Part I (Chapter 1 - 3) provides a quick introduction to R (chapter 1) and presents required R packages and data format (Chapter 2) for clustering analysis and visualization.

The classification of objects, into clusters, requires some methods for measuring the distance or the (dis)similarity between the objects. Chapter 3 covers the common distance measures used for assessing similarity between observations.

Part II starts with partitioning clustering methods, which include:

• K-means clustering (Chapter 4),
• K-Medoids or PAM (partitioning around medoids) algorithm (Chapter 5) and
• CLARA algorithms (Chapter 6).

Partitioning clustering approaches subdivide the data sets into a set of k groups, where k is the number of groups pre-specified by the analyst.

cluster analysis in R

In Part III, we consider agglomerative hierarchical clustering method, which is an alternative approach to partitionning clustering for identifying groups in a data set. It does not require to pre-specify the number of clusters to be generated. The result of hierarchical clustering is a tree-based representation of the objects, which is also known as dendrogram (see the figure below).

In this part, we describe how to compute, visualize, interpret and compare dendrograms:

• Agglomerative clustering (Chapter 7)
  • Algorithm and steps
  • Verify the cluster tree
  • Cut the dendrogram into different groups
• Compare dendrograms (Chapter 8)
  • Visual comparison of two dendrograms
  • Correlation matrix between a list of dendrograms
• Visualize dendrograms (Chapter 9)
  • Case of small data sets
  • Case of dendrogram with large data sets: zoom, sub-tree, PDF
  • Customize dendrograms using dendextend
• Heatmap: static and interactive (Chapter 10)
  • R base heat maps
  • Pretty heat maps
  • Interactive heat maps
  • Complex heatmap
  • Real application: gene expression data

In this section, you will learn how to generate and interpret the following plots.

• Standard dendrogram with filled rectangle around clusters:

cluster analysis in R

• Compare two dendrograms:

cluster analysis in R

• Heatmap:

cluster analysis in R

Part IV describes clustering validation and evaluation strategies, which consists of measuring the goodness of clustering results. Before applying any clustering algorithm to a data set, the first thing to do is to assess the clustering tendency. That is, whether applying clustering is suitable for the data. If yes, then how many clusters are there. Next, you can perform hierarchical clustering or partitioning clustering (with a pre-specified number of clusters). Finally, you can use a number of measures, described in this chapter, to evaluate the goodness of the clustering results.

The different chapters included in part IV are organized as follow:

• Assessing clustering tendency (Chapter 11)

• Determining the optimal number of clusters (Chapter 12)

• Cluster validation statistics (Chapter 13)

• Choosing the best clustering algorithms (Chapter 14)

• Computing p-value for hierarchical clustering (Chapter 15)

In this section, you’ll learn how to create and interpret the plots hereafter.

• Visual assessment of clustering tendency (left panel): Clustering tendency is detected in a visual form by counting the number of square shaped dark blocks along the diagonal in the image.
• Determine the optimal number of clusters (right panel) in a data set using the gap statistics.

cluster analysis in R

• Cluster validation using the silhouette coefficient (Si): A value of Si close to 1 indicates that the object is well clustered. A value of Si close to -1 indicates that the object is poorly clustered. The figure below shows the silhouette plot of a k-means clustering.

cluster analysis in R

Part V presents advanced clustering methods, including:

• Hierarchical k-means clustering (Chapter 16)
• Fuzzy clustering (Chapter 17)
• Model-based clustering (Chapter 18)
• DBSCAN: Density-Based Clustering (Chapter 19)

The hierarchical k-means clustering is an hybrid approach for improving k-means results.

In Fuzzy clustering, items can be a member of more than one cluster. Each item has a set of membership coefficients corresponding to the degree of being in a given cluster.

In model-based clustering, the data are viewed as coming from a distribution that is mixture of two ore more clusters. It finds best fit of models to data and estimates the number of clusters.

The density-based clustering (DBSCAN is a partitioning method that has been introduced in Ester et al. (1996). It can find out clusters of different shapes and sizes from data containing noise and outliers.

cluster analysis in R

1 How this article is organized

We’ll start by demonstrating why we should combine k-means and hierarcical clustering. An application is provided using R software.

Finally, we’ll provide an easy to use R function (in factoextra package) for computing hybrid hierachical k-means clustering.

2 Required R packages

We’ll use the R package factoextra which is very helpful for simplifying clustering workflows and for visualizing clusters using ggplot2 plotting system

Install factoextra package as follow:

if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/factoextra")

Load the package:

library(factoextra)

3 Data preparation

We’ll use USArrest dataset and we start by scaling the data:

# Load the data
data(USArrests)
# Scale the data
df <- scale(USArrests)
head(df)

##                Murder   Assault   UrbanPop         Rape
## Alabama    1.24256408 0.7828393 -0.5209066 -0.003416473
## Alaska     0.50786248 1.1068225 -1.2117642  2.484202941
## Arizona    0.07163341 1.4788032  0.9989801  1.042878388
## Arkansas   0.23234938 0.2308680 -1.0735927 -0.184916602
## California 0.27826823 1.2628144  1.7589234  2.067820292
## Colorado   0.02571456 0.3988593  0.8608085  1.864967207

If you want to understand why the data are scaled before the analysis, then you should read this section: Distances and scaling.

4 R function for clustering analyses

We’ll use the function eclust() [in factoextra] which provides several advantages as described in the previous chapter: Visual Enhancement of Clustering Analysis.

eclust() stands for enhanced clustering. It simplifies the workflow of clustering analysis and, it can be used for computing hierarchical clustering and partititioning clustering in a single line function call.

library("factoextra")
# K-means clustering
km.res <- eclust(df, "kmeans", k = 4,
                 nstart = 25, graph = FALSE)
# k-means group number of each observation
head(km.res$cluster, 15)

##     Alabama      Alaska     Arizona    Arkansas  California    Colorado 
##           4           3           3           4           3           3 
## Connecticut    Delaware     Florida     Georgia      Hawaii       Idaho 
##           2           2           3           4           2           1 
##    Illinois     Indiana        Iowa 
##           3           2           1

# Visualize k-means clusters
fviz_cluster(km.res,  frame.type = "norm", frame.level = 0.68)

# Visualize the silhouette of clusters
fviz_silhouette(km.res)

##   cluster size ave.sil.width
## 1       1   13          0.37
## 2       2   16          0.34
## 3       3   13          0.27
## 4       4    8          0.39

# Silhouette width of observation
sil <- km.res$silinfo$widths[, 1:3]
# Objects with negative silhouette
neg_sil_index <- which(sil[, &#39;sil_width&#39;] < 0)
sil[neg_sil_index, , drop = FALSE]

##          cluster neighbor   sil_width
## Missouri       3        2 -0.07318144

# Enhanced hierarchical clustering
res.hc <- eclust(df, "hclust", k = 4,
                method = "ward.D2", graph = FALSE) 
head(res.hc$cluster, 15)

##     Alabama      Alaska     Arizona    Arkansas  California    Colorado 
##           1           2           2           3           2           2 
## Connecticut    Delaware     Florida     Georgia      Hawaii       Idaho 
##           3           3           2           1           3           4 
##    Illinois     Indiana        Iowa 
##           2           3           4

# Dendrogram
fviz_dend(res.hc, rect = TRUE, show_labels = TRUE, cex = 0.5) 

# Visualize the silhouette of clusters
fviz_silhouette(res.hc)

##   cluster size ave.sil.width
## 1       1    7          0.46
## 2       2   12          0.29
## 3       3   19          0.26
## 4       4   12          0.43

# Silhouette width of observation
sil <- res.hc$silinfo$widths[, 1:3]
# Objects with negative silhouette
neg_sil_index <- which(sil[, &#39;sil_width&#39;] < 0)
sil[neg_sil_index, , drop = FALSE]

##          cluster neighbor   sil_width
## Kentucky       3        4 -0.06459230
## Arkansas       3        1 -0.08467352

5 Combining hierarchical clustering and k-means

res.hc <- eclust(df, "hclust", k = 4,
                method = "ward.D2", graph = FALSE) 
grp <- res.hc$cluster

# Compute cluster centers
clus.centers <- aggregate(df, list(grp), mean)
clus.centers

##   Group.1     Murder    Assault   UrbanPop        Rape
## 1       1  1.5803956  0.9662584 -0.7775109  0.04844071
## 2       2  0.7298036  1.1188219  0.7571799  1.32135653
## 3       3 -0.3621789 -0.3444705  0.3953887 -0.21863180
## 4       4 -1.0782511 -1.1370610 -0.9296640 -1.00344660

# Remove the first column
clus.centers <- clus.centers[, -1]
clus.centers

##       Murder    Assault   UrbanPop        Rape
## 1  1.5803956  0.9662584 -0.7775109  0.04844071
## 2  0.7298036  1.1188219  0.7571799  1.32135653
## 3 -0.3621789 -0.3444705  0.3953887 -0.21863180
## 4 -1.0782511 -1.1370610 -0.9296640 -1.00344660

km.res2 <- eclust(df, "kmeans", k = clus.centers, graph = FALSE)
fviz_silhouette(km.res2)

##   cluster size ave.sil.width
## 1       1    8          0.39
## 2       2   13          0.27
## 3       3   16          0.34
## 4       4   13          0.37

# res.hc$cluster: Initial clusters defined using hierarchical clustering
# km.res2$cluster: Final clusters defined using k-means
table(km.res2$cluster, res.hc$cluster)

##    
##      1  2  3  4
##   1  7  0  1  0
##   2  0 12  1  0
##   3  0  0 16  0
##   4  0  0  1 12

fviz_dend(res.hc, k = 4, 
          k_colors = c("black", "red",  "blue", "green3"),
          label_cols =  km.res2$cluster[res.hc$order], cex = 0.6)

# Final clusters defined using hierarchical k-means clustering
km.clust <- km.res$cluster

# Standard k-means clustering
set.seed(123)
res.km <- kmeans(df, centers = 4, iter.max = 100)


# comparison
table(km.clust, res.km$cluster)

##         
## km.clust  1  2  3  4
##        1 13  0  0  0
##        2  0 16  0  0
##        3  0  0 13  0
##        4  0  0  0  8

# Compute hierarchical k-means clustering
res.hk <-hkmeans(df, 4)
# Elements returned by hkmeans()
names(res.hk)

##  [1] "cluster"      "centers"      "totss"        "withinss"    
##  [5] "tot.withinss" "betweenss"    "size"         "iter"        
##  [9] "ifault"       "data"         "hclust"

# Print the results
res.hk

## Hierarchical K-means clustering with 4 clusters of sizes 8, 13, 16, 13
## 
## Cluster means:
##       Murder    Assault   UrbanPop        Rape
## 1  1.4118898  0.8743346 -0.8145211  0.01927104
## 2  0.6950701  1.0394414  0.7226370  1.27693964
## 3 -0.4894375 -0.3826001  0.5758298 -0.26165379
## 4 -0.9615407 -1.1066010 -0.9301069 -0.96676331
## 
## Clustering vector:
##        Alabama         Alaska        Arizona       Arkansas     California 
##              1              2              2              1              2 
##       Colorado    Connecticut       Delaware        Florida        Georgia 
##              2              3              3              2              1 
##         Hawaii          Idaho       Illinois        Indiana           Iowa 
##              3              4              2              3              4 
##         Kansas       Kentucky      Louisiana          Maine       Maryland 
##              3              4              1              4              2 
##  Massachusetts       Michigan      Minnesota    Mississippi       Missouri 
##              3              2              4              1              2 
##        Montana       Nebraska         Nevada  New Hampshire     New Jersey 
##              4              4              2              4              3 
##     New Mexico       New York North Carolina   North Dakota           Ohio 
##              2              2              1              4              3 
##       Oklahoma         Oregon   Pennsylvania   Rhode Island South Carolina 
##              3              3              3              3              1 
##   South Dakota      Tennessee          Texas           Utah        Vermont 
##              4              1              2              3              4 
##       Virginia     Washington  West Virginia      Wisconsin        Wyoming 
##              3              3              4              4              3 
## 
## Within cluster sum of squares by cluster:
## [1]  8.316061 19.922437 16.212213 11.952463
##  (between_SS / total_SS =  71.2 %)
## 
## Available components:
## 
##  [1] "cluster"      "centers"      "totss"        "withinss"    
##  [5] "tot.withinss" "betweenss"    "size"         "iter"        
##  [9] "ifault"       "data"         "hclust"

# Visualize the tree
fviz_dend(res.hk, cex = 0.6, rect = TRUE)

# Visualize the hkmeans final clusters
fviz_cluster(res.hk, frame.type = "norm", frame.level = 0.68)

6 Infos

This analysis has been performed using R software (ver. 3.2.4)

1 Required packages

The following R packages are required in this chapter:

• factoextra for data visualization
• clustertend for assessing clustering tendency
• seriation for visually assessment of cluster tendency

1. factoextra can be installed as follow:

if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/factoextra")

2. Install clustertend and seriation:

install.packages("clustertend")
install.packages("seriation")

3. Load required packages:

library(factoextra)
library(clustertend)
library(seriation)

2 Data preparation

We’ll use two datasets: the built-in R dataset faithful and a simulated dataset.

# Load the data
data("faithful")
df <- faithful
head(df)

##   eruptions waiting
## 1     3.600      79
## 2     1.800      54
## 3     3.333      74
## 4     2.283      62
## 5     4.533      85
## 6     2.883      55

library("ggplot2")
ggplot(df, aes(x=eruptions, y=waiting)) +
  geom_point() +  # Scatter plot
  geom_density_2d() # Add 2d density estimation

# Generate random dataset
set.seed(123)
n <- nrow(df)

random_df <- data.frame(
  x = runif(nrow(df), min(df$eruptions), max(df$eruptions)),
  y = runif(nrow(df), min(df$waiting), max(df$waiting)))

# Plot the data
ggplot(random_df, aes(x, y)) + geom_point()

random_df <- apply(df, 2, 
                function(x, n){runif(n, min(x), (max(x)))}, n)

3 Why assessing clustering tendency?

As shown above, we know that faithful dataset contains 2 real clusters. However the randomly generated dataset doesn’t contain any meaningful clusters.

The R code below computes k-means clustering and/or hierarchical clustering on the two datasets. The function fviz_cluster() and fviz_dend() [in factoextra] will be used to visualize the results.

library(factoextra)
set.seed(123)
# K-means on faithful dataset
km.res1 <- kmeans(df, 2)
fviz_cluster(list(data = df, cluster = km.res1$cluster),
             frame.type = "norm", geom = "point", stand = FALSE)

# K-means on the random dataset
km.res2 <- kmeans(random_df, 2)
fviz_cluster(list(data = random_df, cluster = km.res2$cluster),
             frame.type = "norm", geom = "point", stand = FALSE)

# Hierarchical clustering on the random dataset
fviz_dend(hclust(dist(random_df)), k = 2,  cex = 0.5)

It can be seen that, k-means algorithm and hierarchical clustering impose a classification on the random uniformly distributed dataset even if there are no meaningful clusters present in it.

Clustering tendency assessment methods are used to avoid this issue.

4 Methods for assessing clustering tendency

Clustering tendency assessment determines whether a given dataset contains meaningful clusters (i.e., non-random structure).

In this section, we’ll describe two methods for determining the clustering tendency: i) a statistical (Hopkins statistic) and ii) a visual methods (Visual Assessment of cluster Tendency (VAT) algorithm).

hopkins(data, n, byrow = F, header = F)

library(clustertend)
# Compute Hopkins statistic for faithful dataset
set.seed(123)
hopkins(faithful, n = nrow(faithful)-1)

## $H
## [1] 0.1588201

# Compute Hopkins statistic for a random dataset
set.seed(123)
hopkins(random_df, n = nrow(random_df)-1)

## $H
## [1] 0.5388899

library("seriation")
# faithful data: ordered dissimilarity image
df_scaled <- scale(faithful)
df_dist <- dist(df_scaled) 
dissplot(df_dist)

# faithful data: ordered dissimilarity image
random_df_scaled <- scale(random_df)
random_df_dist <- dist(random_df_scaled) 
dissplot(random_df_dist)

set.seed(123)
km.res <- kmeans(scale(faithful), 2)
dissplot(df_dist, labels = km.res$cluster)

5 A single function for Hopkins statistic and VAT

The function get_clust_tendency() [in factoextra package] can be used to compute Hopkins statistic and provides also an ordered dissimilarity image using ggplot2, in a single function call. The ordering of dissimilarity matrix is done using hierarchical clustering.

# Cluster tendency
clustend <- get_clust_tendency(scale(faithful), 100)
# Hopkins statistic
clustend$hopkins_stat

## [1] 0.1482683

# Customize the plot
clustend$plot + 
  scale_fill_gradient(low = "steelblue", high = "white")

6 Infos

This analysis has been performed using R software (ver. 3.2.4)

1 Introduction

2 How this document is organized?

Here,

• we start by describing the two standard clustering strategies [partitioning methods (k-MEANS, PAM, CLARA) and hierarchical clustering] as well as how to assess the quality of clustering analysis.
• next, we provide a step-by-step guide for clustering analysis and an R package, named factoextra, for ggplot2-based elegant clustering visualization.
• finally, we describe advanced clustering approaches to find pattern of any shape in large data sets with noise and outliers.

To be published late in 2016. Subscribe to our mailing list at: STHDA mailing list. You will be notified about this book.

This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivs 3.0 United States License.

3 Data preparation

The built-in R dataset USArrest is used as demo data.

• Remove missing data
• Scale variables to make them comparable

# Load data
data("USArrests")
my_data <- USArrests

# Remove any missing value (i.e, NA values for not available)
my_data <- na.omit(my_data)

# Scale variables
my_data <- scale(my_data)

# View the firt 3 rows
head(my_data, n = 3)

##             Murder   Assault   UrbanPop         Rape
## Alabama 1.24256408 0.7828393 -0.5209066 -0.003416473
## Alaska  0.50786248 1.1068225 -1.2117642  2.484202941
## Arizona 0.07163341 1.4788032  0.9989801  1.042878388

4 Installing and loading required R packages

1. Install required packages

# Install factoextra
install.packages("factoextra")

# Install cluster package
install.packages("cluster")

2. Loading required packages

library("cluster")
library("factoextra")

5 Clarifying distance measures

In this chapter, we covered the common distance measures used for assessing similarity between observations. Some R codes, for computing and visualizing pairwise-distances between observations, are also provided.

How this chapter is organized?

• Methods for measuring distances
• Distances and scaling
• Data preparation
• R functions for computing distances
  • The standard dist() function
  • Correlation based distance measures
  • The function daisy() in cluster package
• Visualizing distance matrices

Read more: Clarifying distance measures.

It’s simple to compute and visualize distance matrix using the functions get_dist() and fviz_dist() in factoextra R package:

• get_dist(): for computing a distance matrix between the rows of a data matrix. Compared to the standard dist() function, it supports correlation-based distance measures including “pearson”, “kendall” and “spearman” methods.

• fviz_dist(): for visualizing a distance matrix

res.dist <- get_dist(USArrests, stand = TRUE, method = "pearson")

fviz_dist(res.dist, 
   gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07"))

Clustering - Unsupervised Machine Learning

Read more: Clarifying distance measures.

6 Basic clustering methods

6.1 Partitioning clustering

Partitioning algorithms are clustering approaches that split the data sets, containing n observations, into a set of k groups (i.e. clusters). The algorithms require the analyst to specify the number of clusters to be generated.

This chapter describes the most commonly used partitioning algorithms including:

• K-means clustering (MacQueen, 1967), in which, each cluster is represented by the center or means of the data points belonging to the cluster.
• K-medoids clustering or PAM (Partitioning Around Medoids, Kaufman & Rousseeuw, 1990), in which, each cluster is represented by one of the objects in the cluster. It’s a “non-parametric” alternative of k-means clustering. We’ll describe also a variant of PAM named CLARA (Clustering Large Applications) which is used for analyzing large data sets.

For each of these methods, we provide:

• the basic idea and the key mathematical concepts
• the clustering algorithm and implementation in R software
• R lab sections with many examples for computing clustering methods and visualizing the outputs

Clustering - Unsupervised Machine Learning

How this chapter is organized?

1. Required packages: cluster (for computing clustering algorithm) and factoextra (for elegant visualization)
2. K-means clustering
   • Concept
   • Algorithm
   • R function for k-means clustering: stats::kmeans()
   • Data format
   • Compute k-means clustering
   • Application of K-means clustering on real data
     • Data preparation and descriptive statistics
     • Determine the number of optimal clusters in the data: factoextra::fviz_nbclust()
     • Compute k-means clustering
     • Plot the result: factoextra::fviz_cluster()
3. PAM: Partitioning Around Medoids
   • Concept
   • Algorithm
   • R function for computing PAM: cluster::pam() or fpc::pamk()
   • Compute PAM
4. CLARA: Clustering Large Applications
   • Concept
   • Algorithm
   • R function for computing CLARA: cluster::clara()
5. R packages and functions for visualizing partitioning clusters
   • cluster::clusplot() function
   • factoextra::fviz_cluster() function

Read more: Partitioning cluster analysis. If you are in hurry, read the following quick-start guide.

• K-means clustering: split the data into a set of k groups (i.e., cluster), where k must be specified by the analyst. Each cluster is represented by means of points belonging to the cluster.

Determine the optimal number of clusters: use factoextra::fviz_nbclust()

library("factoextra")
fviz_nbclust(my_data, kmeans, method = "gap_stat")

Clustering - Unsupervised Machine Learning

Compute and visualize k-means clustering

km.res <- kmeans(my_data, 4, nstart = 25)

# Visualize
library("factoextra")
fviz_cluster(km.res, data = my_data, frame.type = "convex")+
  theme_minimal()

Clustering - Unsupervised Machine Learning

• PAM clustering: Partitioning Around Medoids. Robust alternative to k-means clustering, less sensitive to outliers.

# Compute PAM
library("cluster")
pam.res <- pam(my_data, 4)

# Visualize
fviz_cluster(pam.res)

Read more: Partitioning cluster analysis.

6.2 Hierarchical clustering

Hierarchical clustering is an alternative approach to k-means clustering for identifying groups in the dataset. It does not require to pre-specify the number of clusters to be generated.

Hierarchical clustering can be subdivided into two types:

• Agglomerative hierarchical clustering (AHC) in which, each observation is initially considered as a cluster of its own (leaf. Then, the most similar clusters are iteratively merged until there is just one single big cluster (root).
• Divise hierarchical clustering which is an inverse of AHC. It begins with the root, in witch all objects are included in one cluster. Then the most heterogeneous clusters are iteratively divided until all observation are in their own cluster.

The result of hierarchical clustering is a tree-based representation of the observations which is called a dendrogram. Observations can be subdivided into groups by cutting the dendogram at a desired similarity level.

This chapter provides:

• The description of the different types of hierarchical clustering algorithms
• R lab sections with many examples for computing hierarchical clustering, visualizing and comparing dendrogram
• The interpretation of dendrogram
• R codes for cutting the dendrograms into groups

How this chapter is organized?

1. Required R packages
2. Algorithm
3. Data preparation and descriptive statistics
4. R functions for hierarchical clustering
   • hclust() function
   • agnes() and diana() functions
5. Interpretation of the dendrogram
6. Cut the dendrogram into different groups
7. Hierarchical clustering and correlation based distance
8. What type of distance measures should we choose?
9. Comparing two dendrograms
   • Tanglegram
   • Correlation matrix between a list of dendrogram

Read more: Hierarchical clustering essentials. If you are in hurry, read the following quick-start guide.

1. Install and load required packages (cluster, factoextra) as previously described

2. Compute and visualize hierarchical clustering using R base functions

# 1. Loading and preparing data
data("USArrests")
my_data <- scale(USArrests)

# 2. Compute dissimilarity matrix
d <- dist(my_data, method = "euclidean")

# Hierarchical clustering using Ward&#39;s method
res.hc <- hclust(d, method = "ward.D2" )

# Cut tree into 4 groups
grp <- cutree(res.hc, k = 4)

# Visualize
plot(res.hc, cex = 0.6) # plot tree
rect.hclust(res.hc, k = 4, border = 2:5) # add rectangle

Clustering - Unsupervised Machine Learning

3. Elegant visualization using factoextra functions: factoextra::hcut(), factoextra::fviz_dend()

library("factoextra")
# Compute hierarchical clustering and cut into 4 clusters
res <- hcut(USArrests, k = 4, stand = TRUE)

# Visualize
fviz_dend(res, rect = TRUE, cex = 0.5,
          k_colors = c("#00AFBB","#2E9FDF", "#E7B800", "#FC4E07"))

Clustering - Unsupervised Machine Learning

We’ll see also, how to customize the dendrogram:

Clustering - Unsupervised Machine Learning

Read more: Hierarchical clustering essentials.

7 Clustering validation

The aim of this part is to:

• describe the different methods for clustering validation
• compare the quality of clustering results obtained with different clustering algorithms
• provide R lab section for validating clustering results

7.1 Assessing clustering tendency

Assessing clustering tendency consists of examining whether the data is clusterable, that is, whether the data contains any inherent grouping structure. This should be checked before applying clustering analysis.

In this chapter:

• We describe why we should evaluate the clustering tendency before applying any cluster analysis on a dataset.
• We describe statistical and visual methods for assessing the clustering tendency
• R lab sections containing many examples are also provided for computing clustering tendency and visualizing clusters

How this chapter is organized?

1. Required packages
2. Data preparation
3. Why assessing clustering tendency?
4. Methods for assessing clustering tendency
   • Hopkins statistic
     • Algorithm
     • R function for computing Hopkins statistic: clustertend::hopkins()
   • VAT: Visual Assessment of cluster Tendency: seriation::dissplot()
     • VAT Algorithm
     • R functions for VAT
5. A single function for Hopkins statistic and VAT: factoextra::get_clust_tendency()

Read more: Assessing clustering tendency. If you are in hurry, read the following quick-start guide.

1. Install and load factoextra as previously described

2. Assessing clustering tendency: use factoextra::get_clust_tendency(). Assess clustering tendency using Hopkins’ statistic and a visual approach. An ordered dissimilarity image (ODI) is shown.

• Hopkins statistic: If the value of Hopkins statistic is close to zero (far below 0.5), then we can conclude that the dataset is significantly clusterable.

• VAT (Visual Assessment of cluster Tendency): The VAT detects the clustering tendency in a visual form by counting the number of square shaped dark (or colored) blocks along the diagonal in a VAT image.

library("factoextra")
my_data <- scale(iris[, -5])
get_clust_tendency(my_data, n = 50,
                   gradient = list(low = "steelblue",  high = "white"))

## $hopkins_stat
## [1] 0.2002686
## 
## $plot

Clustering - Unsupervised Machine Learning

Read more: Assessing clustering tendency.

7.2 Determining the optimal number of clusters

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

In this chapter, we’ll describe different methods to determine the optimal number of clusters for k-means, PAM and hierarchical clustering.

How this chapter is organized?

1. Required packages
2. Data preparation
3. Example of partitioning method results
4. Example of hierarchical clustering results
5. Three popular methods for determining the optimal number of clusters
   • Elbow method
     • Concept
     • Algorithm
     • R codes
   • Average silhouette method
     • Concept
     • Algorithm
     • R codes
   • Conclusions about elbow and silhouette methods
   • Gap statistic method
     • Concept
     • Algorithm
     • R codes
6. NbClust: A Package providing 30 indices for determining the best number of clusters
   • Overview of NbClust package
   • NbClust R function
   • Examples of usage
     • Compute only an index of interest
     • Compute all the 30 indices

Read more: Determining the optimal number of clusters. If you are in hurry, read the following quick-start guide.

• Estimate the number of clusters in the data using gap statistics : factoextra::fviz_nbclust()

my_data <- scale(USArrests)
library("factoextra")
fviz_nbclust(my_data, kmeans, method = "gap_stat")

Clustering - Unsupervised Machine Learning

• NbClust: A Package providing 30 indices for determining the best number of clusters

library("NbClust")
set.seed(123)
res.nbclust <- NbClust(my_data, distance = "euclidean",
                  min.nc = 2, max.nc = 10, 
                  method = "complete", index ="all") 

Visualize using factoextra:

factoextra::fviz_nbclust(res.nbclust) + theme_minimal()

## Among all indices: 
## ===================
## * 2 proposed  0 as the best number of clusters
## * 1 proposed  1 as the best number of clusters
## * 9 proposed  2 as the best number of clusters
## * 4 proposed  3 as the best number of clusters
## * 6 proposed  4 as the best number of clusters
## * 2 proposed  5 as the best number of clusters
## * 1 proposed  8 as the best number of clusters
## * 1 proposed  10 as the best number of clusters
## 
## Conclusion
## =========================
## * According to the majority rule, the best number of clusters is  2 .

Clustering - Unsupervised Machine Learning

Read more: Determining the optimal number of clusters.

7.3 Clustering validation statistics

A variety of measures has been proposed in the literature for evaluating clustering results. The term clustering validation is used to design the procedure of evaluating the results of a clustering algorithm.

The aim of this chapter is to:

• describe the different methods for clustering validation
• compare the quality of clustering results obtained with different clustering algorithms
• provide R lab section for validating clustering results

How this chapter is organized?

1. Required packages: cluster, factoextra, NbClust, fpc
2. Data preparation
3. Relative measures - Determine the optimal number of clusters: NbClust::NbClust()
4. Clustering analysis
   • Example of partitioning method results
   • Example of hierarchical clustering results
5. Internal clustering validation measures
   • Silhouette analysis
     • Concept and algorithm
     • Interpretation of silhouette width
     • R functions for silhouette analysis: cluster::silhouette(), factoextra::fviz_silhouette()
   • Dunn index
     • Concept and algorithm
     • R function for computing Dunn index: fpc::cluster.stats(), NbClust::NbClust()
   • Clustering validation statistics: fpc::cluster.stats()
6. External clustering validation

Read more: Clustering Validation Statistics. If you are in hurry, read the following quick-start guide.

1. Compute and visualize hierarchical clustering

• Compute: factoextra::eclust()
• Elegant visualization: factoextra::fviz_dend()

my_data <- scale(iris[, -5])

# Enhanced hierarchical clustering, cut in 3 groups
library("factoextra")
res.hc <- eclust(my_data, "hclust", k = 3, graph = FALSE) 

# Visualize
fviz_dend(res.hc, rect = TRUE, show_labels = FALSE)

Clustering - Unsupervised Machine Learning

2. Validate clustering results by inspection the cluster silhouette plot

Recall that the silhouette (\(S_i\)) measures how similar an object \(i\) is to the the other objects in its own cluster versus those in the neighbor cluster. \(S_i\) values range from 1 to - 1:

• A value of \(S_i\) close to 1 indicates that the object is well clustered. In the other words, the object \(i\) is similar to the other objects in its group.
• A value of \(S_i\) close to -1 indicates that the object is poorly clustered, and that assignment to some other cluster would probably improve the overall results.

# Visualize the silhouette plot
fviz_silhouette(res.hc)

##   cluster size ave.sil.width
## 1       1   49          0.63
## 2       2   30          0.44
## 3       3   71          0.32

Clustering - Unsupervised Machine Learning

Which samples have negative silhouette? To what cluster are they closer?

# Silhouette width of observations
sil <- res.hc$silinfo$widths[, 1:3]

# Objects with negative silhouette
neg_sil_index <- which(sil[, &#39;sil_width&#39;] < 0)
sil[neg_sil_index, , drop = FALSE]

##     cluster neighbor   sil_width
## 84        3        2 -0.01269799
## 122       3        2 -0.01789603
## 62        3        2 -0.04756835
## 135       3        2 -0.05302402
## 73        3        2 -0.10091884
## 74        3        2 -0.14761137
## 114       3        2 -0.16107155
## 72        3        2 -0.23036371

Read more: Clustering Validation Statistics.

my_data <- scale(USArrests)

# Compute clValid
library("clValid")
intern <- clValid(my_data, nClust = 2:6, 
              clMethods = c("hierarchical","kmeans","pam"),
              validation = "internal")
# Summary
summary(intern)

## 
## Clustering Methods:
##  hierarchical kmeans pam 
## 
## Cluster sizes:
##  2 3 4 5 6 
## 
## Validation Measures:
##                                  2       3       4       5       6
##                                                                   
## hierarchical Connectivity   6.6437  9.5615 13.9563 22.5782 31.2873
##              Dunn           0.2214  0.2214  0.2224  0.2046  0.2126
##              Silhouette     0.4085  0.3486  0.3637  0.3213  0.2720
## kmeans       Connectivity   6.6437 13.6484 16.2413 24.6639 33.7194
##              Dunn           0.2214  0.2224  0.2224  0.1983  0.2231
##              Silhouette     0.4085  0.3668  0.3573  0.3377  0.3079
## pam          Connectivity   6.6437 13.8302 20.4421 29.5726 38.2643
##              Dunn           0.2214  0.1376  0.1849  0.1849  0.2019
##              Silhouette     0.4085  0.3144  0.3390  0.3105  0.2630
## 
## Optimal Scores:
## 
##              Score  Method       Clusters
## Connectivity 6.6437 hierarchical 2       
## Dunn         0.2231 kmeans       6       
## Silhouette   0.4085 hierarchical 2

7.5 How to compute p-value for hierarchical clustering in R?

This chapter describes the R package pvclust (Suzuki et al., 2004) which uses bootstrap resampling techniques to compute p-value for each clusters.

How this chapter is organized?

1. Concept
2. Algorithm
3. Required R packages
4. Data preparation
5. Compute p-value for hierarchical clustering
   • Description of pvclust() function
   • Usage of pvclust() function

Read more: How to compute p-value for hierarchical clustering in R?. If you are in hurry, read the following quick-start guide.

Note that, pvclust() performs clustering on the columns of the dataset, which correspond to samples in our case.

library(pvclust)
# Data preparation
set.seed(123)
data("lung")
ss <- sample(1:73, 30) # extract 20 samples out of
my_data <- lung[, ss]

# Compute pvclust
res.pv <- pvclust(my_data, method.dist="cor", 
                  method.hclust="average", nboot = 10)

## Bootstrap (r = 0.5)... Done.
## Bootstrap (r = 0.6)... Done.
## Bootstrap (r = 0.7)... Done.
## Bootstrap (r = 0.8)... Done.
## Bootstrap (r = 0.9)... Done.
## Bootstrap (r = 1.0)... Done.
## Bootstrap (r = 1.1)... Done.
## Bootstrap (r = 1.2)... Done.
## Bootstrap (r = 1.3)... Done.
## Bootstrap (r = 1.4)... Done.

# Default plot
plot(res.pv, hang = -1, cex = 0.5)
pvrect(res.pv)

Clustering - Unsupervised Machine Learning

Clusters with AU > = 95% are indicated by the rectangles and are considered to be strongly supported by data.

Read more: How to compute p-value for hierarchical clustering in R?.

8 The guide for clustering analysis on a real data: 4 steps you should know

Read more: The guide for clustering analysis on a real data: 4 steps you should know.

9 Visualization of clustering results

In this chapter, we’ll describe how to visualize the result of clustering using dendrograms as well as static and interactive heatmap.

Heat map is a false color image with a dendrogram added to the left side and to the top. It’s used to visualize a hidden pattern in a data matrix in order to reveal some associations between rows or columns.

9.2 Beautiful dendrogram visualizations

Read more: Beautiful dendrogram visualizations in R: 5+ must known methods

Clustering - Unsupervised Machine Learning

9.3 Static and Interactive Heatmap

Read more: Static and Interactive Heatmap in R

Clustering - Unsupervised Machine Learning

10 Advanced clustering methods

10.2 Model-based clustering

In model-based clustering, the data are viewed as coming from a distribution that is mixture of two ore more clusters. It finds best fit of models to data and estimates the number of clusters. Read more: Model-based clustering.

Clustering - Unsupervised Machine Learning

10.3 DBSCAN: Density-based clustering

DBSCAN is a partitioning method that has been introduced in Ester et al. (1996). It can find out clusters of different shapes and sizes from data containing noise and outliers.The basic idea behind density-based clustering approach is derived from a human intuitive clustering method.

The description and implementation of DBSCAN in R are provided in this chapter : DBSCAN.

Clustering - Unsupervised Machine Learning

10.4 Hybrid clustering methods

• Hybrid hierarchical k-means clustering for optimizing clustering outputs - Hybrid approach (1/1)
• HCPC: Hierarchical clustering on principal components - Hybrid approach (2/2)

Clustering - Unsupervised Machine Learning

11 Infos

This analysis has been performed using R software (ver. 3.2.4)

• Martin Ester, Hans-Peter Kriegel, Joerg Sander, Xiaowei Xu (1996). A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise. Institute for Computer Science, University of Munich. Proceedings of 2nd International Conference on Knowledge Discovery and Data Mining (KDD-96). pdf

1 Required packages

Three R packages are required for this chapter:

1. cluster and e1071 for computing fuzzy clustering
2. factoextra for visualizing clusters

install.packages("cluster")
install.packages("e1071")
install.packages("factoextra")

2 Concept of fuzzy clustering

In K-means or PAM clustering, the data is divided into distinct clusters, where each element is affected exactly to one cluster. This type of clustering is also known as hard clustering or non-fuzzy clustering. Unlike K-means, Fuzzy clustering is considered as a soft clustering, in which each element has a probability of belonging to each cluster. In other words, each element has a set of membership coefficients corresponding to the degree of being in a given cluster.

Points close to the center of a cluster, may be in the cluster to a higher degree than points in the edge of a cluster. The degree, to which an element belongs to a given cluster, is a numerical value in [0, 1].

Fuzzy c-means (FCM) algorithm is one of the most widely used fuzzy clustering algorithms. It was developed by Dunn in 1973 and improved by Bezdek in 1981. It’s frequently used in pattern recognition.

3 Algorithm of fuzzy clustering

FCM algorithm is very similar to the k-means algorithm and the aim is to minimize the objective function defined as follow:

\[ \sum\limits_{j=1}^k \sum\limits_{x_i \in C_j} u_{ij}^m (x_i - \mu_j)^2 \]

Where,

• \(u_{ij}\) is the degree to which an observation \(x_i\) belongs to a cluster \(c_j\)
• \(\mu_j\) is the center of the cluster j
• \(u_{ij}\) is the degree to which an observation \(x_i\) belongs to a cluster \(c_j\)
• \(m\) is the fuzzifier.

It can be seen that, FCM differs from k-means by using the membership values \(u_{ij}\) and the fuzzifier \(m\).

The variable \(u_{ij}^m\) is defined as follow:

\[ u_{ij}^m = \frac{1}{\sum\limits_{l=1}^k \left( \frac{| x_i - c_j |}{| x_i - c_k |}\right)^{\frac{2}{m-1}}} \]

The degree of belonging, \(u_{ij}\), is linked inversely to the distance from x to the cluster center.

The parameter \(m\) is a real number greater than 1 (\(1.0 < m < \infty\)) and it defines the level of cluster fuzziness. Note that, a value of \(m\) close to 1 gives a cluster solution which becomes increasingly similar to the solution of hard clustering such as k-means; whereas a value of \(m\) close to infinite leads to complete fuzzyness.

Note that, a good choice is to use m = 2.0 (Hathaway and Bezdek 2001).

In fuzzy clustering the centroid of a cluster is he mean of all points, weighted by their degree of belonging to the cluster:

\[ C_j = \frac{\sum\limits_{x \in C_j} u_{ij}^m x}{\sum\limits_{x \in C_j} u_{ij}^m} \]

Where,

• \(C_j\) is the centroid of the cluster j
• \(u_{ij}\) is the degree to which an observation \(x_i\) belongs to a cluster \(c_j\)

The algorithm of fuzzy clustering can be summarize as follow:

1. Specify a number of clusters k (by the analyst)
2. Assign randomly to each point coefficients for being in the clusters.
3. Repeat until the maximum number of iterations (given by “maxit”) is reached, or when the algorithm has converged (that is, the coefficients’ change between two iterations is no more than \(\epsilon\), the given sensitivity threshold):
   • Compute the centroid for each cluster, using the formula above.
   • For each point, compute its coefficients of being in the clusters, using the formula above.

The algorithm minimizes intra-cluster variance as well, but has the same problems as k-means; the minimum is a local minimum, and the results depend on the initial choice of weights. Hence, different initializations may lead to different results.

Using a mixture of Gaussians along with the expectation-maximization algorithm is a more statistically formalized method which includes some of these ideas: partial membership in classes.

fanny(x, k, memb.exp = 2, metric = "euclidean", 
      stand = FALSE, maxit = 500)

library(cluster)
set.seed(123)
# Load the data
data("USArrests")

# Subset of USArrests
ss <- sample(1:50, 20)
df <- scale(USArrests[ss,])

# Compute fuzzy clustering
res.fanny <- fanny(df, 4)

# Cluster plot using fviz_cluster()
# You can use also : clusplot(res.fanny)
library(factoextra)
fviz_cluster(res.fanny, frame.type = "norm",
             frame.level = 0.68)

# Silhouette plot
fviz_silhouette(res.fanny, label = TRUE)

##   cluster size ave.sil.width
## 1       1    4          0.52
## 2       2    6          0.10
## 3       3    6          0.41
## 4       4    4          0.04

print(res.fanny)

## Fuzzy Clustering object of class 'fanny' :                      
## m.ship.expon.        2
## objective     6.052789
## tolerance        1e-15
## iterations         215
## converged            1
## maxit              500
## n                   20
## Membership coefficients (in %, rounded):
##              [,1] [,2] [,3] [,4]
## Iowa           75   11    7    7
## Rhode Island   26   32   21   21
## Maryland        8   19   37   37
## Tennessee      10   24   33   33
## Utah           23   36   20   20
## Arizona        10   23   34   34
## Mississippi    16   25   29   29
## Wisconsin      65   15   10   10
## Virginia       17   37   23   23
## Maine          63   15   11   11
## Texas           8   25   33   33
## Louisiana       9   22   35   35
## Montana        41   26   17   17
## Michigan        8   20   36   36
## Arkansas       19   30   25   25
## New York        9   24   34   34
## Florida        10   21   35   35
## Alaska         15   24   31   31
## Hawaii         27   34   20   20
## New Jersey     16   37   23   23
## Fuzzyness coefficients:
## dunn_coeff normalized 
## 0.31337355 0.08449807 
## Closest hard clustering:
##         Iowa Rhode Island     Maryland    Tennessee         Utah 
##            1            2            3            4            2 
##      Arizona  Mississippi    Wisconsin     Virginia        Maine 
##            3            4            1            2            1 
##        Texas    Louisiana      Montana     Michigan     Arkansas 
##            3            4            1            3            2 
##     New York      Florida       Alaska       Hawaii   New Jersey 
##            3            3            4            2            2 
## 
## Available components:
##  [1] "membership"  "coeff"       "memb.exp"    "clustering"  "k.crisp"    
##  [6] "objective"   "convergence" "diss"        "call"        "silinfo"    
## [11] "data"

# Membership coefficient
res.fanny$membership

##                    [,1]      [,2]       [,3]       [,4]
## Iowa         0.75234997 0.1056742 0.07098791 0.07098791
## Rhode Island 0.26129280 0.3198982 0.20940449 0.20940449
## Maryland     0.07559096 0.1906031 0.36690296 0.36690296
## Tennessee    0.10351700 0.2444743 0.32600436 0.32600436
## Utah         0.23177048 0.3631831 0.20252321 0.20252321
## Arizona      0.09505979 0.2329621 0.33598906 0.33598906
## Mississippi  0.15957721 0.2511123 0.29465525 0.29465525
## Wisconsin    0.65274007 0.1530047 0.09712764 0.09712764
## Virginia     0.16856415 0.3654879 0.23297397 0.23297397
## Maine        0.62818484 0.1532966 0.10925930 0.10925930
## Texas        0.08407125 0.2465250 0.33470188 0.33470188
## Louisiana    0.09152177 0.2159634 0.34625741 0.34625741
## Montana      0.40788012 0.2556886 0.16821562 0.16821562
## Michigan     0.07811792 0.1957270 0.36307753 0.36307753
## Arkansas     0.19473888 0.2992279 0.25301662 0.25301662
## New York     0.08723572 0.2392572 0.33675356 0.33675356
## Florida      0.09725070 0.2073927 0.34767830 0.34767830
## Alaska       0.14688036 0.2428630 0.30512830 0.30512830
## Hawaii       0.26945561 0.3356724 0.19743602 0.19743602
## New Jersey   0.16160093 0.3720897 0.23315470 0.23315470

# Visualize using corrplot
library(corrplot)
corrplot(res.fanny$membership, is.corr = FALSE)

# Dunn's partition coefficient
res.fanny$coeff

## dunn_coeff normalized 
## 0.31337355 0.08449807

# Observation groups
res.fanny$clustering

##         Iowa Rhode Island     Maryland    Tennessee         Utah 
##            1            2            3            4            2 
##      Arizona  Mississippi    Wisconsin     Virginia        Maine 
##            3            4            1            2            1 
##        Texas    Louisiana      Montana     Michigan     Arkansas 
##            3            4            1            3            2 
##     New York      Florida       Alaska       Hawaii   New Jersey 
##            3            3            4            2            2

cmeans(x, centers, iter.max = 100, dist = "euclidean", m = 2)

set.seed(123)
library(e1071)
cm <- cmeans(df, 4)
cm

## Fuzzy c-means clustering with 4 clusters
## 
## Cluster centers:
##       Murder    Assault   UrbanPop       Rape
## 1  0.6290005  0.9705484  0.5006389  0.8647698
## 2  0.8560350  0.3375298 -0.7294688  0.2002994
## 3 -1.2101485 -1.2476750 -0.7277747 -1.1534135
## 4 -0.7314218 -0.6647441  1.0032068 -0.3335272
## 
## Memberships:
##                        1           2          3          4
## Iowa         0.005939255 0.009155372 0.96585947 0.01904590
## Rhode Island 0.104616576 0.098854401 0.20500209 0.59152694
## Maryland     0.697459281 0.227720539 0.02731256 0.04750762
## Tennessee    0.078024194 0.872296030 0.02111342 0.02856636
## Utah         0.049301432 0.044484100 0.08442894 0.82178552
## Arizona      0.740498081 0.118781050 0.03988867 0.10083220
## Mississippi  0.179555100 0.624367937 0.10296383 0.09311313
## Wisconsin    0.024017906 0.033630983 0.83136508 0.11098604
## Virginia     0.155690387 0.395730684 0.19167059 0.25690834
## Maine        0.021165990 0.034336946 0.89152511 0.05297195
## Texas        0.545608753 0.240753676 0.05410235 0.15953522
## Louisiana    0.275003950 0.617629141 0.04197257 0.06539434
## Montana      0.062161310 0.135620851 0.66557661 0.13664123
## Michigan     0.848927329 0.096168273 0.01784963 0.03705477
## Arkansas     0.131803310 0.565593614 0.18039386 0.12220922
## New York     0.694179984 0.131927283 0.04157413 0.13231860
## Florida      0.711655719 0.173670792 0.03979837 0.07487512
## Alaska       0.369474028 0.381553979 0.11356564 0.13540635
## Hawaii       0.064103932 0.066647766 0.14874490 0.72050340
## New Jersey   0.082015921 0.059546923 0.05743425 0.80100291
## 
## Closest hard clustering:
##         Iowa Rhode Island     Maryland    Tennessee         Utah 
##            3            4            1            2            4 
##      Arizona  Mississippi    Wisconsin     Virginia        Maine 
##            1            2            3            2            3 
##        Texas    Louisiana      Montana     Michigan     Arkansas 
##            1            2            3            1            2 
##     New York      Florida       Alaska       Hawaii   New Jersey 
##            1            1            2            4            4 
## 
## Available components:
## [1] "centers"     "size"        "cluster"     "membership"  "iter"       
## [6] "withinerror" "call"

fviz_cluster(list(data = df, cluster=cm$cluster), frame.type = "norm",
             frame.level = 0.68)

4 Infos

This analysis has been performed using R software (ver. 3.2.4)

• J. C. Dunn (1973): A Fuzzy Relative of the ISODATA Process and Its Use in Detecting Compact Well-Separated Clusters. Journal of Cybernetics 3: 32-57
• J. C. Bezdek (1981): Pattern Recognition with Fuzzy Objective Function Algorithms. Plenum Press, New York Tariq Rashid: “Clustering”

1 Data

The built-in mtcars R data is used:

df <- as.matrix(scale(mtcars))

2 Draw a heat map using R base function

The built-in R heatmap function [in stats package] can be used.

A simplified format is:

heatmap(x, scale = "row")

• x: a numeric matrix
• scale: a character indicating if the values should be centered and scaled in either the row direction or the column direction, or none. Allowed values are in c(“row”, “column”, “none”). Default is “row”.

# Default plot
heatmap(df, scale = "none")

Heatmap - R data visualization

# Use custom colors
col<- colorRampPalette(c("red", "white", "blue"))(256)
heatmap(scale(as.matrix(mtcars)), scale = "none",
        col =  col)

Heatmap - R data visualization

The R code below will customize the heatmap as follow:

1. An RColorBrewer color palette name is used to change the appearance
2. The argument RowSideColors and ColSideColors are used to annotate rows and columns respectively. The expected values for these options are a vector containing color names specifying the classes for rows/columns.

# Use RColorBrewer color palette names
library("RColorBrewer")
col <- colorRampPalette(brewer.pal(10, "RdYlBu"))(256)
heatmap(df, scale = "none", col =  col, 
        RowSideColors = rep(c("blue", "pink"), each = 16),
        ColSideColors = c(rep("purple", 5), rep("orange", 6)))

Heatmap - R data visualization

3 Enhanced heat map

The function heatmap.2() [in gplots package] provides many extensions to the standard R heatmap() function presented in the previous section.

# install.packages("gplots")
library("gplots")
heatmap.2(df, scale = "none", col = bluered(100), 
          trace = "none", density.info = "none")

Heatmap - R data visualization

Other arguments can be used including:

• labRow, labCol
• margins
• hclustfun: hclustfun=function(x) hclust(x, method=“ward”)
• keysize

In the R code above, bluered() function [in gplots package] is used to generate a smoothly varying set of colors. You can also use the following color generator functions:

• colorpanel(n, low, mid, high)
  • n: Desired number of color elements to be generated
  • low, mid, high: Colors to use for the Lowest, middle, and highest values. mid may be omitted.
• redgreen(n)
• greenred(n)
• bluered(n)
• redblue(n)

4 Interactive heatmap

The package d3heatmap can be used to produce an interactive heatmap:

It can be installed as follow:

if (!require("devtools")) install.packages("devtools")
devtools::install_github("rstudio/d3heatmap")

The function d3heatmap() is used to create the interactive heatmap:

The possibilities below are provided:

• Put the mouse on a heatmap cell of interest to view the row and the column names as well as the corresponding value.
• select an area for zooming. After zooming, click on the heatmap again to go back to the previous display

library("d3heatmap")
d3heatmap(scale(mtcars), colors = "RdBu",
          k_row = 4, k_col = 2)

• colors: Either an RColorBrewer color palette name (e.g. “YlOrRd” or “Blues”), or a vector of colors to interpolate in hexadecimal “#RRGGBB” format, or a color interpolation function like colorRamp. Read this: available colors in R
• k_row, k_col: an integer specifying the desired number of groups by which to color the dendrogram’s branches in row and column, respectively.

For further customizing the heatmap read ?d3heatmap. Possible options include:

5 Enhancing heatmaps using dendextend

The package dendextend can be used to enhance functions from other packages. The mtcars data is used in the following sections. We’ll start by defining the order and the appearance for rows and columns using dendextend. These results are used in others functions from others packages.

The order and the appearance for rows and columns can be defined as follow:

library(dendextend)
# order for rows
Rowv  <- mtcars %>% scale %>% dist %>% hclust %>% as.dendrogram %>%
   set("branches_k_color", k = 3) %>% set("branches_lwd", 1.2) %>%
   ladderize

# Order for columns
# We must transpose the data
Colv  <- mtcars %>% scale %>% t %>% dist %>% hclust %>% as.dendrogram %>%
   set("branches_k_color", k = 2, value = c("orange", "blue")) %>%
   set("branches_lwd", 1.2) %>%
   ladderize

The arguments above can be used in the functions below:

1. The standard heatmap() function [in **stats* package]:

heatmap(scale(mtcars), Rowv = Rowv, Colv = Colv,
        scale = "none")

2. The enhanced heatmap.2() function [in gplots package]:

library(gplots)
heatmap.2(scale(mtcars), scale = "none", col = bluered(100), 
          Rowv = Rowv, Colv = Colv,
          trace = "none", density.info = "none")

Heatmap - R data visualization

3. The interactive heatmap generator d3heatmap() function [in d3heatmap package]:

library("d3heatmap")
d3heatmap(scale(mtcars), colors = "RdBu",
          Rowv = Rowv, Colv = Colv)

6 Complex heatmap

ComplexHeatmap is an R/bioconductor package, developed by Zuguang Gu, which provides a flexible solution to arrange and annotate multiple heatmaps. It allows also to visualize the association between different data from different sources.

if (!require("devtools")) install.packages("devtools")
devtools::install_github("jokergoo/ComplexHeatmap")

library("ComplexHeatmap")

Heatmap(matrix, col, name)

6.3 Single heatmap

A single heatmap can be used to visualize a data set containing continuous or discrete values.

In the example below we’ll visualize the built-in mtcars data set.

Recall that, the mtcars data comprises fuel consumption and 10 aspects of automobile design and performance for 32 automobiles (1973–74 models).

data(mtcars)
head(mtcars[, 1:6])

##                    mpg cyl disp  hp drat    wt
## Mazda RX4         21.0   6  160 110 3.90 2.620
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875
## Datsun 710        22.8   4  108  93 3.85 2.320
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215
## Hornet Sportabout 18.7   8  360 175 3.15 3.440
## Valiant           18.1   6  225 105 2.76 3.460

Before drawing the heatmap, the data is firstly scaled, using the R base scale() function.

df <- scale(mtcars)
Heatmap(df, name = "mtcars")

Heatmap - R data visualization

6.3.1 Colors

The argument col is used to specify colors. As our data matrix contains continuous values, the option col should be a color mapping function. In this case, the colorRamp2() function [in circlize] can be used as follow:

library("circlize")
Heatmap(df, name = "mtcars",
        col = colorRamp2(c(-2, 0, 2), c("green", "white", "red")))

Heatmap - R data visualization

The two arguments for colorRamp2() is a vector of breaks values and corresponding colors.

It’s also possible to use RColorBrewer color palettes:

library("RColorBrewer")
Heatmap(df, name = "mtcars",
        col = colorRamp2(c(-2, 0, 2), brewer.pal(n=3, name="RdBu")))

Heatmap - R data visualization

In the next sections, we’ll use the following custom color palette:

mycol <- colorRamp2(c(-2, 0, 2), c("blue", "white", "red"))

6.3.2 Titles

The heatmap name, column title and row title can be changed as follow:

Heatmap(df, name = "mtcars", col = mycol,
        column_title = "Column title",
        row_title = "Row title")

Heatmap - R data visualization

Note that, the default side for the row title is “left” and the default side for the column title is “top”. This can be changed using the options below:

• row_title_side: Allowed values are “left” or “right” (e.g.: row_title_side = “right” )
• column_title_side: Allowed values are “top” or “bottom” (e.g.: column_title_side = “bottom” )

It’s also possible to modify the font size and face of titles using the options:

• row_title_gp: graphic parameters for drawing row text
• column_title_gp: graphic parameters for drawing column text

For instance,

Heatmap(df, name = "mtcars", col = mycol,
        column_title = "Column title",
        column_title_gp = gpar(fontsize = 14, fontface = "bold"),
        row_title = "Row title",
        row_title_gp = gpar(fontsize = 14, fontface = "bold"))

Heatmap - R data visualization

In the R code above, the possible values for fontface can be an integer or string: 1 = plain, 2 = bold, 3 = italic and 4 = bold italic. If a string, then valid values are: “plain”, “bold”, “italic”, “oblique”, and “bold.italic”.

6.3.3 Row and column names

1. Show row/column names:
   • show_row_names: whether to show row names. Default value is TRUE
   • show_column_names: whether to show column names. Default value is TRUE

Heatmap(df, name = "mtcars", show_row_names = FALSE)

2. Change font size and face:
   • row_names_gp: graphical parameters for drawing row names
   • column_names_gp: graphical parameters for drawing column names

Heatmap(df, name = "mtcars", 
        row_names_gp = gpar(fontsize = 14, fontface = "bold",
                            col = c("blue", "red")))

6.3.4 Clustering

6.3.4.1 Change the appearance of clustering

By default, rows and columns are clustered. This can be inactivated using the argument:

• cluster_rows = FALSE. If TRUE, makes cluster on rows
• cluster_columns = FALSE. If TRUE, makes cluster on columns

Cluster on rows are inactivated using the R code below:

# Inactivate cluster on rows
Heatmap(df, name = "mtcars", col = mycol, cluster_rows = FALSE)

Heatmap - R data visualization

In some cases, we want to make the cluster on rows/columns, but we don’t want to show the dendogram on the final image. In this case, use the options:

• show_row_hclust: logical value; whether to show row clusters
• show_column_hclust: logical value; whether to show column clusters

It’s also possible to change the side of row and column clusters using the arguments:

• row_hclust_side: The allowed values are “left” or “right”
• column_hclust_side: The allowed values are “top” or “bottom”

If you want to change the height of column clusters or the width of row clusters, you can use the option column_dend_height and row_dend_width as follow:

Heatmap(df, name = "mtcars", col = mycol,
        column_dend_height = unit(2, "cm"),
        row_dend_width = unit(2, "cm") )

Heatmap - R data visualization

We can also customize the appearance of dendogram using the function color_branches() [in dendextend package]:

# install.packages("dendextend")
library(dendextend)
row_dend = hclust(dist(df)) # row clustering
col_dend = hclust(dist(t(df))) # column clustering
Heatmap(df, name = "mtcars", col = mycol,
        cluster_rows = color_branches(row_dend, k = 4),
        cluster_columns = color_branches(col_dend, k = 2))

Heatmap - R data visualization

6.3.4.2 Metric for clustering

The arguments clustering_distance_rows and clustering_distance_columns are used to specify the metric for row and column clustering, respectively. Default values are “euclidean”.

Allowed values are:

1. A Pre-defined character which is in (“euclidean”, “maximum”, “manhattan”, “canberra”, “binary”, “minkowski”, “pearson”, “spearman”, “kendall”):

Heatmap(df, name = "mtcars", clustering_distance_rows = "pearson",
        clustering_distance_columns = "pearson")

2. A Pre-defined function, such as dist(), to calculate distance from matrix (m):

Heatmap(df, name = "mtcars", 
        clustering_distance_rows = function(m) dist(m))

3. A Self-defined function which calculates distance from two vectors:

Heatmap(df, name = "mtcars", 
        clustering_distance_rows = function(x, y) 1 - cor(x, y))

Note that, in the R code above, the example is generally shown for the argument clustering_distance_rows which specify the metric for row clustering. I recommend to use the same metric for the argument clustering_distance_columns (metric for column clustering).

As an illustration, the R code below applies a self defined function for clustering which is robust to outliers based on the pair-wise distance:

# Clustering metric function
robust_dist = function(x, y) {
    qx = quantile(x, c(0.1, 0.9))
    qy = quantile(y, c(0.1, 0.9))
    l = x > qx[1] & x < qx[2] & y > qy[1] & y < qy[2]
    x = x[l]
    y = y[l]
    sqrt(sum((x - y)^2))
}
# Heatmap
Heatmap(df, name = "mtcars", 
    clustering_distance_rows = robust_dist,
    clustering_distance_columns = robust_dist,
    col = colorRamp2(c(-2, 0, 2), c("purple", "white", "orange")))

Heatmap - R data visualization

6.3.4.3 Clustering methods

The arguments clustering_method_rows and clustering_method_columns can be used to specify the method for making hierarchical clustering. Allowed values are those supported by hclust() function including “ward.D”, “ward.D2”, “single”, “complete”, “average”, … (see ?hclust).

As an example:

Heatmap(df, name = "mtcars", clustering_method_rows = "ward.D",
        clustering_method_columns = "ward.D")

6.3.5 Split heatmap by rows

There are many ways to split the heatmap. One solution is to apply k-means using the argument km.

It’s important to use the set.seed() function when performing k-means so that the results obtained can be reproduced precisely at a later time.

set.seed(2)
# split into 2 groups
Heatmap(df, name = "mtcars", col = mycol, k = 2)

Heatmap - R data visualization

It’s also possible to use split argument to specify row classes as a vector. In the following example we’ll use the levels of the factor variable cyl [in mtcars] to split the heatmap by rows. Recall that cyl corresponds to the number of cylinders.

# split by a vector specifying row classes
Heatmap(df, name = "mtcars", col = mycol, 
        split = mtcars$cyl )

Heatmap - R data visualization

Note that, split can be also a data frame in which different combinations of levels split the rows of the heatmap.

# Split by combining multiple variables
Heatmap(df, name ="mtcars", col = mycol,
        split = data.frame(cyl = mtcars$cyl, am = mtcars$am))

Heatmap - R data visualization

# Combine km and split
Heatmap(df, name ="mtcars", col = mycol,
        km = 2, split =  mtcars$cyl)

Heatmap - R data visualization

If you want to use other partitioning method, rather than k-means, you can easily do it by just assigning the partitioning vector to split. In the R code below, we’ll use pam() function [in cluster package]. pam() stands for Partitioning (clustering) of the data into k clusters “around medoids”, a more robust version of K-means.

# install.packages("cluster")
library("cluster")
set.seed(2)
pa = pam(df, k = 3)
Heatmap(df, name = "mtcars", col = mycol,
        split = paste0("pam", pa$clustering))

Heatmap - R data visualization

It’s also possible to combine user defined dendrograms and split. In this case, split can be specified as a single number:

library(dendextend)
row_dend = hclust(dist(df)) # row clustering
row_dend = color_branches(row_dend, k = 4)
Heatmap(df, name = "mtcars", col = mycol,
        cluster_rows = row_dend, split = 2)

Heatmap - R data visualization

6.4 Heatmap annotation

The HeatmapAnnotation class is used to define annotation on row or column. A simplified format is:

HeatmapAnnotation(df, name, col, show_legend)

• df: a data.frame with column names
• name: the name of the heatmap annotation
• col: a list of colors which contains color mapping to columns in df

For the example below, we’ll transpose our data to have the observations in columns and the variables in rows.

6.4.1 Prepare the data

# Transpose
df <- t(df)
# Heatmap of the transposed data
Heatmap(df, name ="mtcars", col = mycol)

Heatmap - R data visualization

6.4.2 Simple annotation

In simple annotation a vector, containing discrete or continuous values, is used to annotate rows or columns.

We’ll use the qualitative variables cyl (levels = “4”, “5” and “8”) and am (levels = “0” and “1”), and the continuous variable mpg to annotate columns.

For each of these 3 variables, custom colors are defined as follow:

# Annotation data frame
annot_df <- data.frame(cyl = mtcars$cyl, am = mtcars$am,  mpg = mtcars$mpg)

# Define colors for each levels of qualitative variables
# Define gradient color for continuous variable (mpg)
col = list(cyl = c("4" = "green", "6" = "gray", "8" = "darkred"),
            am = c("0" = "yellow", "1" = "orange"),
            mpg = colorRamp2(c(17, 25), c("lightblue", "purple"))
            )

# Create the heatmap annotation
ha <- HeatmapAnnotation(annot_df, col = col)

# Combine the heatmap and the annotation
Heatmap(df, name = "mtcars", col = mycol,
        top_annotation = ha)

Heatmap - R data visualization

It’s possible to hide the annotation legend using the argument show_legend = FALSE as follow:

ha <- HeatmapAnnotation(annot_df, col = col, show_legend = FALSE)
Heatmap(df, name = "mtcars", col = mycol, top_annotation = ha)

Annotation names can be added using the R code hereafter. The function qq() [in GetoptLong package], for simple variable interpolation in texts, is required.

library("GetoptLong")
# Combine Heatmap and annotation
ha <- HeatmapAnnotation(annot_df, col = col, show_legend = FALSE)
Heatmap(df, name = "mtcars", col = mycol, top_annotation = ha)
# Add annotation names on the right
for(an in colnames(annot_df)) {
    seekViewport(qq("annotation_@{an}"))
    grid.text(an, unit(1, "npc") + unit(2, "mm"), 0.5,
              default.units = "npc", just = "left")
}

Heatmap - R data visualization

To add annotation names on the left, use the code below:

# Annotation names on the left
for(an in colnames(annot_df)) {
    seekViewport(qq("annotation_@{an}"))
    grid.text(an, unit(1, "npc") - unit(2, "mm"), 0.5,
              default.units = "npc", just = "left")
}

6.4.3 Complex annotation

In this section we’ll see how to combine heatmap and some basic graphs to show the data distribution. For simple annotation graphics, the following functions can be used: anno_points(), anno_barplot(), anno_boxplot(), anno_density() and anno_histogram().

An example is shown below:

# Define some graphics to display the distribution of columns
.hist = anno_histogram(df, gp = gpar(fill = "lightblue"))
.density = anno_density(df, type = "line", gp = gpar(col = "blue"))
ha_mix_top = HeatmapAnnotation(hist = .hist, density = .density)

# Define some graphics to display the distribution of rows
.violin = anno_density(df, type = "violin", 
                       gp = gpar(fill = "lightblue"), which = "row")
.boxplot = anno_boxplot(df, which = "row")
ha_mix_right = HeatmapAnnotation(violin = .violin, bxplt = .boxplot,
                              which = "row", width = unit(4, "cm"))

# Combine annotation with heatmap
Heatmap(df, name = "mtcars", col = mycol,
        column_names_gp = gpar(fontsize = 8),
        top_annotation = ha_mix_top, 
        top_annotation_height = unit(4, "cm")) + ha_mix_right

Heatmap - R data visualization

Note that, it’s also possible to use the argument bottom_annotation.

6.5 Combine multiple heatmaps

Multiple heatmaps can be arranged as follow:

# Heatmap 1
ht1 = Heatmap(df, name = "ht1", col = mycol, km = 2,
              column_names_gp = gpar(fontsize = 9))
# Heatmap 2
ht2 = Heatmap(df, name = "ht2", 
        col = colorRamp2(c(-2, 0, 2), c("green", "white", "red")),
        column_names_gp = gpar(fontsize = 9))
# Combine the two heatmaps
ht1 + ht2

Heatmap - R data visualization

You can use the option width = unit(3, “cm”)) to control the size of the heatmaps.

Note that when combining multiple heatmaps, the first heatmap is considered as the main heatmap. Some settings of the remaining heatmaps are auto-adjusted according to the setting of the main heatmap. These include: removing row clusters and titles, and adding splitting

The draw() function can be used to customize the appearance of the final image:

draw(ht1 + ht2, 
    # Titles
    row_title = "Two heatmaps, row title", 
    row_title_gp = gpar(col = "red"),
    column_title = "Two heatmaps, column title", 
    column_title_side = "bottom",
    # Gap between heatmaps
    gap = unit(0.5, "cm"))

Legends can be removed using the arguments show_heatmap_legend = FALSE, show_annotation_legend = FALSE.

6.7 Gene expression matrix

In gene expression data, rows are genes and columns are samples. More information about genes can be attached after the expression heatmap such as gene length and type of genes.

expr = readRDS(paste0(system.file(package = "ComplexHeatmap"),
                      "/extdata/gene_expression.rds"))
mat = as.matrix(expr[, grep("cell", colnames(expr))])

type = gsub("s\\d+_", "", colnames(mat))
ha = HeatmapAnnotation(df = data.frame(type = type))

Heatmap(mat, name = "expression", km = 5, top_annotation = ha, 
    top_annotation_height = unit(4, "mm"), 
    show_row_names = FALSE, show_column_names = FALSE) +
Heatmap(expr$length, name = "length", width = unit(5, "mm"),
        col = colorRamp2(c(0, 100000), c("white", "orange"))) +
Heatmap(expr$type, name = "type", width = unit(5, "mm")) +
Heatmap(expr$chr, name = "chr", width = unit(5, "mm"),
        col = rand_color(length(unique(expr$chr))))

Heatmap - R data visualization

It’s also possible to visualize genomic alterations and to integrate different molecular levels (gene expression, DNA methylation, …). Read the vignette for further examples.

6.8 Visualize distribution of column in matrix

The function densityHeatmap() is used.

densityHeatmap(scale(mtcars))

Heatmap - R data visualization

The dashed lines on the heatmap correspond to the five quantile numbers. The text for the five quantile levels are added in the right of the heatmap.

7 Infos

This analysis has been performed using R software (ver. 3.2.3) and ComplexHeatmap (ver. )

1 Required packages

The following packages will be used:

• cluster for clustering analyses
• factoextra for visualizing clusters using ggplot2 plotting system

Install factoextra package as follow:

if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/factoextra")

The cluster package can be installed using the code below:

install.packages("cluster")

Load packages:

library(cluster)
library(factoextra)

2 Data preparation

We’ll use the built-in R data set USArrests, which can be loaded and prepared as follow:

# Load the data set
data(USArrests)

# Remove any missing value (i.e, NA values for not available)
# That might be present in the data
USArrests <- na.omit(USArrests)

# View the firt 6 rows of the data
head(USArrests, n = 6)

##            Murder Assault UrbanPop Rape
## Alabama      13.2     236       58 21.2
## Alaska       10.0     263       48 44.5
## Arizona       8.1     294       80 31.0
## Arkansas      8.8     190       50 19.5
## California    9.0     276       91 40.6
## Colorado      7.9     204       78 38.7

In this data set, columns are variables and rows are observations (i.e., samples).

To inspect the data before the K-means clustering we’ll compute some descriptive statistics such as the mean and the standard deviation of the variables.

The apply() function is used to apply a given function (e.g : min(), max(), mean(), …) on the data set. The second argument can take the value of:

• 1: for applying the function on the rows
• 2: for applying the function on the columns

desc_stats <- data.frame(
  Min = apply(USArrests, 2, min), # minimum
  Med = apply(USArrests, 2, median), # median
  Mean = apply(USArrests, 2, mean), # mean
  SD = apply(USArrests, 2, sd), # Standard deviation
  Max = apply(USArrests, 2, max) # Maximum
  )
desc_stats <- round(desc_stats, 1)
head(desc_stats)

##           Min   Med  Mean   SD   Max
## Murder    0.8   7.2   7.8  4.4  17.4
## Assault  45.0 159.0 170.8 83.3 337.0
## UrbanPop 32.0  66.0  65.5 14.5  91.0
## Rape      7.3  20.1  21.2  9.4  46.0

Note that the variables have a large different means and variances. They must be standardized to make them comparable.

Standardization consists of transforming the variables such that they have mean zero and standard deviation one. The scale() function can be used as follow:

df<- scale(USArrests)

3 Assessing the clusterability

The function get_clust_tendency() [in factoextra] can be used. It computes Hopkins statistic and provides a visual approach.

library("factoextra")
res <- get_clust_tendency(df, 40, graph = FALSE)
# Hopskin statistic
res$hopkins_stat

## [1] 0.3440875

# Visualize the dissimilarity matrix
res$plot

## NULL

The value of Hopkins statistic is significantly < 0.5, indicating that the data is highly clusterable. Additionally, It can be seen that the ordered dissimilarity image contains patterns (i.e., clusters).

4 Estimate the number of clusters in the data

As k-means clustering requires to specify the number of clusters to generate, we’ll use the function clusGap() [in cluster] to compute gap statistics for estimating the optimal number of clusters . The function fviz_gap_stat() [in factoextra] is used to visualize the gap statistic plot.

library("cluster")
set.seed(123)
# Compute the gap statistic
gap_stat <- clusGap(df, FUN = kmeans, nstart = 25, 
                    K.max = 10, B = 500) 
# Plot the result
library(factoextra)
fviz_gap_stat(gap_stat)

The gap statistic suggests a 4 cluster solutions.

It’s also possible to use the function NbClust() [in NbClust] package.

5 Compute k-means clustering

K-means clustering with k = 4:

# Compute k-means
set.seed(123)
km.res <- kmeans(df, 4, nstart = 25)
head(km.res$cluster, 20)

##     Alabama      Alaska     Arizona    Arkansas  California    Colorado 
##           4           3           3           4           3           3 
## Connecticut    Delaware     Florida     Georgia      Hawaii       Idaho 
##           2           2           3           4           2           1 
##    Illinois     Indiana        Iowa      Kansas    Kentucky   Louisiana 
##           3           2           1           2           1           4 
##       Maine    Maryland 
##           1           3

# Visualize clusters using factoextra
fviz_cluster(km.res, USArrests)

6 Cluster validation statistics: Inspect cluster silhouette plot

Recall that the silhouette measures (\(S_i\)) how similar an object \(i\) is to the the other objects in its own cluster versus those in the neighbor cluster. \(S_i\) values range from 1 to - 1:

• A value of \(S_i\) close to 1 indicates that the object is well clustered. In the other words, the object \(i\) is similar to the other objects in its group.
• A value of \(S_i\) close to -1 indicates that the object is poorly clustered, and that assignment to some other cluster would probably improve the overall results.

sil <- silhouette(km.res$cluster, dist(df))
rownames(sil) <- rownames(USArrests)
head(sil[, 1:3])

##            cluster neighbor  sil_width
## Alabama          4        3 0.48577530
## Alaska           3        4 0.05825209
## Arizona          3        2 0.41548326
## Arkansas         4        2 0.11870947
## California       3        2 0.43555885
## Colorado         3        2 0.32654235

fviz_silhouette(sil)

##   cluster size ave.sil.width
## 1       1   13          0.37
## 2       2   16          0.34
## 3       3   13          0.27
## 4       4    8          0.39

It can be seen that there are some samples which have negative silhouette values. Some natural questions are :

Which samples are these? To what cluster are they closer?

This can be determined from the output of the function silhouette() as follow:

neg_sil_index <- which(sil[, "sil_width"] < 0)
sil[neg_sil_index, , drop = FALSE]

##          cluster neighbor   sil_width
## Missouri       3        2 -0.07318144

7 eclust(): Enhanced clustering analysis

The function eclust() [in factoextra] provides several advantages compared to the standard packages used for clustering analysis:

• It simplifies the workflow of clustering analysis
• It can be used to compute hierarchical clustering and partitioning clustering in a single line function call
• The function eclust() computes automatically the gap statistic for estimating the right number of clusters.
• It automatically provides silhouette information
• It draws beautiful graphs using ggplot2

# Compute k-means
res.km <- eclust(df, "kmeans")

# Gap statistic plot
fviz_gap_stat(res.km$gap_stat)

# Silhouette plot
fviz_silhouette(res.km)

##   cluster size ave.sil.width
## 1       1   13          0.27
## 2       2   13          0.37
## 3       3    8          0.39
## 4       4   16          0.34

 # Enhanced hierarchical clustering
res.hc <- eclust(df, "hclust") # compute hclust
fviz_dend(res.hc, rect = TRUE) # dendrogam

fviz_silhouette(res.hc) # silhouette plot
fviz_cluster(res.hc) # scatter plot

8 Infos

This analysis has been performed using R software (ver. 3.2.3)

1 Concept

The traditional clustering methods such as hierarchical clustering and partitioning algorithms (k-means and others) are heuristic and are not based on formal models.

An alternative is to use model-based clustering, in which, the data are considered as coming from a distribution that is mixture of two or more components (i.e. clusters) (Chris Fraley and Adrian E. Raftery, 2002 and 2012).

Each component k (i.e. group or cluster) is modeled by the normal or Gaussian distribution which is characterized by the parameters:

• \(\mu_k\): mean vector,
• \(\sum_k\): covariance matrix,
• An associated probability in the mixture. Each point has a probability of belonging to each cluster.

2 Model parameters

The model parameters can be estimated using the EM (Expectation-Maximization) algorithm initialized by hierarchical model-based clustering. Each cluster k is centered at the means \(\mu_k\), with increased density for points near the mean.

Geometric features (shape, volume, orientation) of each cluster are determined by the covariance matrix \(\sum_k\).

Different possible parameterizations of \(\sum_k\) are available in the R package mclust (see ?mclustModelNames).

The available model options, in mclust package, are represented by identifiers including: EII, VII, EEI, VEI, EVI, VVI, EEE, EEV, VEV and VVV.

The first identifier refers to volume, the second to shape and the third to orientation. E stands for “equal”, V for “variable” and I for “coordinate axes”.

For example:

• EVI denotes a model in which the volumes of all clusters are equal (E), the shapes of the clusters may vary (V), and the orientation is the identity (I) or “coordinate axes.
• EEE means that the clusters have the same volume, shape and orientation in p-dimensional space.
• VEI means that the clusters have variable volume, the same shape and orientation equal to coordinate axes.

The mclust package uses maximum likelihood to fit all these models, with different covariance matrix parameterizations, for a range of k components. The “best model” is selected using the Bayesian Information Criterion or BIC. A large BIC score indicates strong evidence for the corresponding model.

3 Advantage of model-based clustering

The key advantage of model-based approach, compared to the standard clustering methods (k-means, hierarchical clustering, …), is the suggestion of the number of clusters and an appropriate model.

4 Example of data

We’ll use the bivariate faithful data set which contains the waiting time between eruptions and the duration of the eruption for the Old Faithful geyser in Yellowstone National Park (Wyoming, USA).

# Load the data
data("faithful")
head(faithful)

##   eruptions waiting
## 1     3.600      79
## 2     1.800      54
## 3     3.333      74
## 4     2.283      62
## 5     4.533      85
## 6     2.883      55

An illustration of the data can be drawn using ggplot2 package as follow:

library("ggplot2")
ggplot(faithful, aes(x=eruptions, y=waiting)) +
  geom_point() +  # Scatter plot
  geom_density2d() # Add 2d density estimation

5 Mclust(): R function for computing model-based clustering

The function Mclust() [in mclust package] can be used to compute model-based clustering.

Install and load the package as follow:

# Install
install.packages("mclust")

# Load
library("mclust")

The function Mclust() provides the optimal mixture model estimation according to BIC. A simplified format is:

Mclust(data, G = NULL)

The function Mclust() returns an object of class ‘Mclust’ containing the following elements:

• modelName: A character string denoting the model at which the optimal BIC occurs.
• G: The optimal number of mixture components (i.e: number of clusters)
• BIC: All BIV values
• bic Optimal BIC value
• loglik: The loglikelihood corresponding to the optimal BIC
• df: The number of estimated parameters
• Z: A matrix whose \([i,k]^{th}\) entry is the probability that observation \(i\) in the test data belongs to the \(k^{th}\) class. Column names are cluster numbers, and rows are observations
• classification: The cluster number of each observation, i.e. map(z)
• uncertainty: The uncertainty associated with the classification

6 Example of cluster analysis using Mclust()

library(mclust)
# Model-based-clustering
mc <- Mclust(faithful)
# Print a summary
summary(mc)

## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm 
## ----------------------------------------------------
## 
## Mclust EEE (ellipsoidal, equal volume, shape and orientation) model with 3 components:
## 
##  log.likelihood   n df       BIC       ICL
##       -1126.361 272 11 -2314.386 -2360.865
## 
## Clustering table:
##   1   2   3 
## 130  97  45

# Values returned by Mclust()
names(mc)

##  [1] "call"           "data"           "modelName"      "n"             
##  [5] "d"              "G"              "BIC"            "bic"           
##  [9] "loglik"         "df"             "hypvol"         "parameters"    
## [13] "z"              "classification" "uncertainty"

# Optimal selected model
mc$modelName

## [1] "EEE"

# Optimal number of cluster
mc$G

## [1] 3

# Probality for an observation to be in a given cluster
head(mc$z)

##           [,1]         [,2]         [,3]
## 1 2.181744e-02 1.130837e-08 9.781825e-01
## 2 2.475031e-21 1.000000e+00 3.320864e-13
## 3 2.521625e-03 2.051823e-05 9.974579e-01
## 4 6.553336e-14 9.999998e-01 1.664978e-07
## 5 9.838967e-01 7.642900e-20 1.610327e-02
## 6 2.104355e-07 9.975388e-01 2.461029e-03

# Cluster assignement of each observation
head(mc$classification, 10)

##  1  2  3  4  5  6  7  8  9 10 
##  3  2  3  2  1  2  1  3  2  1

# Uncertainty associated with the classification
head(mc$uncertainty)

##            1            2            3            4            5 
## 2.181745e-02 3.321787e-13 2.542143e-03 1.664978e-07 1.610327e-02 
##            6 
## 2.461239e-03

Model-based clustering results can be drawn using the function plot.Mclust():

plot(x, what = c("BIC", "classification", "uncertainty", "density"),
     xlab = NULL, ylab = NULL, addEllipses = TRUE, main = TRUE, ...)

# BIC values used for choosing the number of clusters
plot(mc, "BIC")

# Classification: plot showing the clustering
plot(mc, "classification")

# Classification uncertainty
plot(mc, "uncertainty")

# Estimated density. Contour plot
plot(mc, "density")

Clusters generated by Mclust() can be drawn using the function fviz_cluster() [in factoextra package]. Read more about [factoextra](https://www.sthda.com/english/wiki/factoextra-r-package-quick-multivariate-data-analysis-pca-ca-mca-and-visualization-r-software-and-data-mining.

library(factoextra)
fviz_cluster(mc, frame.type = "norm", geom = "point")

7 Infos

This analysis has been performed using R software (ver. 3.2.3)

• Chris Fraley, A. E. Raftery, T. B. Murphy and L. Scrucca (2012). mclust Version 4 for R: Normal Mixture Modeling for Model-Based Clustering, Classification, and Density Estimation. Technical Report No. 597, Department of Statistics, University of Washington. pdf
• Chris Fraley and A. E. Raftery (2002). Model-based clustering, discriminant analysis, and density estimation. Journal of the American Statistical Association 97:611:631.

1 plot.hclust(): R base function

As you already know, the standard R function plot.hclust() can be used to draw a dendrogram from the results of hierarchical clustering analyses (computed using hclust() function).

A simplified format is:

plot(x, labels = NULL, hang = 0.1, 
     main = "Cluster dendrogram", sub = NULL,
     xlab = NULL, ylab = "Height", ...)

# Default plot
plot(hc)

# Put the labels at the same height: hang = -1
plot(hc, hang = -1, cex = 0.6)

2 plot.dendrogram() function

In order to visualize the result of a hierarchical clustering analysis using the function plot.dendrogram(), we must firstly convert it as a dendrogram.

The format of the function plot.dendrogram() is:

plot(x, type = c("rectangle", "triangle"), horiz = FALSE)

# Convert hclust into a dendrogram and plot
hcd <- as.dendrogram(hc)
# Default plot
plot(hcd, type = "rectangle", ylab = "Height")

# Triangle plot
plot(hcd, type = "triangle", ylab = "Height")

# Zoom in to the first dendrogram
plot(hcd, xlim = c(1, 20), ylim = c(1,8))

The above dendrogram can be customized using the arguments:

• nodePar: a list of plotting parameters to use for the nodes (see ?points). Default value is NULL. The list may contain components named pch, cex, col, xpd, and/or bg each of which can have length two for specifying separate attributes for inner nodes and leaves.
• edgePar: a list of plotting parameters to use for the edge segments (see ?segments). The list may contain components named col, lty and lwd (for the segments). As with nodePar, each can have length two for differentiating leaves and inner nodes.
• leaflab: a string specifying how leaves are labeled. The default “perpendicular” write text vertically; “textlike” writes text horizontally (in a rectangle), and “none” suppresses leaf labels.

# Define nodePar
nodePar <- list(lab.cex = 0.6, pch = c(NA, 19), 
                cex = 0.7, col = "blue")
# Customized plot; remove labels
plot(hcd, ylab = "Height", nodePar = nodePar, leaflab = "none")

# Horizontal plot
plot(hcd,  xlab = "Height",
     nodePar = nodePar, horiz = TRUE)

# Change edge color
plot(hcd,  xlab = "Height", nodePar = nodePar, 
     edgePar = list(col = 2:3, lwd = 2:1))

3 Phylogenetic trees

The package ape (Analyses of Phylogenetics and Evolution) can be used to produce a more sophisticated dendrogram.

The function plot.phylo() can be used for plotting a dendrogram. A simplified format is:

plot(x, type = "phylogram", show.tip.label = TRUE,
     edge.color = "black", edge.width = 1, edge.lty = 1,
     tip.color = "black")

# install.packages("ape")
library("ape")
# Default plot
plot(as.phylo(hc), cex = 0.6, label.offset = 0.5)

# Cladogram
plot(as.phylo(hc), type = "cladogram", cex = 0.6, 
     label.offset = 0.5)

# Unrooted
plot(as.phylo(hc), type = "unrooted", cex = 0.6,
     no.margin = TRUE)

# Fan
plot(as.phylo(hc), type = "fan")

# Radial
plot(as.phylo(hc), type = "radial")

# Cut the dendrogram into 4 clusters
colors = c("red", "blue", "green", "black")
clus4 = cutree(hc, 4)
plot(as.phylo(hc), type = "fan", tip.color = colors[clus4],
     label.offset = 1, cex = 0.7)

# Change the appearance
# change edge and label (tip)
plot(as.phylo(hc), type = "cladogram", cex = 0.6,
     edge.color = "steelblue", edge.width = 2, edge.lty = 2,
     tip.color = "steelblue")

4 ggdendro package : ggplot2 and dendrogram

The R package ggdendro can be used to extract the plot data from dendrogram and for drawing a dendrogram using ggplot2.

install.packages("ggdendro")

library("ggplot2")
library("ggdendro")

# Visualization using the default theme named theme_dendro()
ggdendrogram(hc)

# Rotate the plot and remove default theme
ggdendrogram(hc, rotate = TRUE, theme_dendro = FALSE)

# Build dendrogram object from hclust results
dend <- as.dendrogram(hc)

# Extract the data (for rectangular lines)
# Type can be "rectangle" or "triangle"
dend_data <- dendro_data(dend, type = "rectangle")
# What contains dend_data
names(dend_data)

## [1] "segments"    "labels"      "leaf_labels" "class"

# Extract data for line segments
head(dend_data$segments)

##           x         y     xend      yend
## 1 19.771484 13.516242 8.867188 13.516242
## 2  8.867188 13.516242 8.867188  6.461866
## 3  8.867188  6.461866 4.125000  6.461866
## 4  4.125000  6.461866 4.125000  2.714554
## 5  4.125000  2.714554 2.500000  2.714554
## 6  2.500000  2.714554 2.500000  1.091092

# Extract data for labels
head(dend_data$labels)

##   x y          label
## 1 1 0        Alabama
## 2 2 0      Louisiana
## 3 3 0        Georgia
## 4 4 0      Tennessee
## 5 5 0 North Carolina
## 6 6 0    Mississippi

# Plot line segments and add labels
p <- ggplot(dend_data$segments) + 
  geom_segment(aes(x = x, y = y, xend = xend, yend = yend))+
  geom_text(data = dend_data$labels, aes(x, y, label = label),
            hjust = 1, angle = 90, size = 3)+
  ylim(-3, 15)
print(p)

5 dendextend package: Extending R’s dendrogram functionality

The package dendextend contains many functions for changing the appearance of a dendrogram and for comparing dendrograms.

In this section we’ll use the chaining operator (%>%) to simplify our code.

data <- scale(USArrests)
dist.res <- dist(data)
hc <- hclust(dist.res, method = "ward.D2")
dend <- as.dendrogram(hc)
plot(dend)

dend <- USArrests[1:5,] %>% # data
        scale %>% # Scale the data
        dist %>% # calculate a distance matrix, 
        hclust(method = "ward.D2") %>% # Hierarchical clustering 
        as.dendrogram # Turn the object into a dendrogram.
plot(dend)

install.packages('dendextend')

library(dendextend)

set(object, what, value)

# Create a dendrogram and plot it
dend <- USArrests[1:5,] %>%  scale %>% 
        dist %>% hclust %>% as.dendrogram

dend %>% plot

# Get the labels of the tree
labels(dend)

## [1] "Alaska"     "Arizona"    "California" "Alabama"    "Arkansas"

# Change the labels, and then plot:
dend %>% set("labels", c("a", "b", "c", "d", "e")) %>% plot

# Change color and size for labels
dend %>% set("labels_col", c("green", "blue")) %>% # change color
  set("labels_cex", 2) %>% # Change size
  plot(main = "Change the color \nand size") # plot

# Color labels by specifying the number of cluster (k)
dend %>% set("labels_col", value = c("green", "blue"), k=2) %>% 
          plot(main = "Color labels \nper cluster")
abline(h = 2, lty = 2)

# Change the type, the color and the size of node points
# +++++++++++++++++++++++++++++
dend %>% set("nodes_pch", 19) %>%  # node point type
  set("nodes_cex", 2) %>%  # node point size
  set("nodes_col", "blue") %>% # node point color
  plot(main = "Node points")

# Change the type, the color and the size of leave points
# +++++++++++++++++++++++++++++
dend %>% set("leaves_pch", 19) %>%  # node point type
  set("leaves_cex", 2) %>%  # node point size
  set("leaves_col", "blue") %>% # node point color
  plot(main = "Leaves points")

# Specify different point types and colors for each leave
dend %>% set("leaves_pch", c(17, 18, 19)) %>%  # node point type
  set("leaves_cex", 2) %>%  # node point size
  set("leaves_col", c("blue", "red", "green")) %>% #node point color
  plot(main = "Leaves points")

# Default colors
dend %>% set("branches_k_color", k = 2) %>% 
  plot(main = "Default colors")

# Customized colors
dend %>% set("branches_k_color", 
             value = c("red", "blue"), k = 2) %>% 
   plot(main = "Customized colors")

# Vertical plot
dend %>% set("branches_k_color", k = 3) %>% plot
dend %>% rect.dendrogram(k=3, border = 8, lty = 5, lwd = 2)

# Horizontal plot
dend %>% set("branches_k_color", k = 3) %>% plot(horiz = TRUE)
dend %>% rect.dendrogram(k = 3, horiz = TRUE, border = 8, lty = 5, lwd = 2)

grp <- c(1,1,1, 2,2)
k_3 <- cutree(dend,k = 3, order_clusters_as_data = FALSE) 
# The FALSE above makes sure we get the clusters in the order of the
# dendrogram, and not in that of the original data. It is like:
# cutree(dend, k = 3)[order.dendrogram(dend)]

the_bars <- cbind(grp, k_3)

dend %>% set("labels", "") %>% plot
colored_bars(colors = the_bars, dend = dend)

dend <- iris[1:30,-5] %>% scale %>% dist %>% 
   hclust %>% as.dendrogram %>%
   set("branches_k_color", k=3) %>% set("branches_lwd", 1.2) %>%
   set("labels_colors") %>% set("labels_cex", c(.9,1.2)) %>% 
   set("leaves_pch", 19) %>% set("leaves_col", c("blue", "red"))
# plot the dend in usual "base" plotting engine:
plot(dend)

library(ggplot2)
# Rectangle dendrogram using ggplot2
ggd1 <- as.ggdend(dend)
ggplot(ggd1) 

# Change the theme to the default ggplot2 theme
ggplot(ggd1, horiz = TRUE, theme = NULL) 

# Theme minimal
ggplot(ggd1, theme = theme_minimal()) 

# Create a radial plot and remove labels
ggplot(ggd1, labels = FALSE) + 
  scale_y_reverse(expand = c(0.2, 0)) +
  coord_polar(theta="x")

library(pvclust)
data(lung) # 916 genes for 73 subjects
set.seed(1234)
result <- pvclust(lung[1:100, 1:10], method.dist="cor", 
                  method.hclust="average", nboot=10)

## Bootstrap (r = 0.5)... Done.
## Bootstrap (r = 0.6)... Done.
## Bootstrap (r = 0.7)... Done.
## Bootstrap (r = 0.8)... Done.
## Bootstrap (r = 0.9)... Done.
## Bootstrap (r = 1.0)... Done.
## Bootstrap (r = 1.1)... Done.
## Bootstrap (r = 1.2)... Done.
## Bootstrap (r = 1.3)... Done.
## Bootstrap (r = 1.4)... Done.

# Default plot of the result
plot(result)
pvrect(result)

# pvclust and dendextend
result %>% as.dendrogram %>% 
  set("branches_k_color", k = 2, value = c("purple", "orange")) %>%
  plot
result %>% text
result %>% pvrect

6 Infos

This analysis has been performed using R software (ver. 3.2.1)

1 Required packages

The following package will be used:

• cluster for computing pam and for analyzing cluster silhouettes
• factoextra for visualizing clusters using ggplot2 plotting system
• NbClust for finding the optimal number of clusters

Install factoextra package as follow:

if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/factoextra")

The remaining packages can be installed using the code below:

pkgs <- c("cluster",  "NbClust")
install.packages(pkgs)

Load packages:

library(factoextra)
library(cluster)
library(NbClust)

2 Data preparation

The data set iris is used. We start by excluding the species column and scaling the data using the function scale():

# Load the data
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])

This iris data set gives the measurements in centimeters of the variables sepal length and width and petal length and width, respectively, for 50 flowers from each of 3 species of iris. The species are Iris setosa, versicolor, and virginica.

3 Example of partitioning method results

The functions kmeans() [in stats package] and pam() [in cluster package] are described in this section. We’ll split the data into 3 clusters as follow:

# K-means clustering
set.seed(123)
km.res <- kmeans(iris.scaled, 3, nstart = 25)
# k-means group number of each observation
km.res$cluster

##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 3 3 3 2 3 3 3 3 3 3 3 3 2 3 3 3 3
##  [71] 2 3 3 3 3 2 2 2 3 3 3 3 3 3 3 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 2 3 2 2 2
## [106] 2 3 2 2 2 2 2 2 3 3 2 2 2 2 3 2 3 2 3 2 2 3 2 2 2 2 2 2 3 3 2 2 2 3 2
## [141] 2 2 3 2 2 2 3 2 2 3

# Visualize k-means clusters
fviz_cluster(km.res, data = iris.scaled, geom = "point",
             stand = FALSE, frame.type = "norm")

# PAM clustering
library("cluster")
pam.res <- pam(iris.scaled, 3)
pam.res$cluster

##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 3 3 3 2 3 3 3 3 3 3 3 3 2 3 3 3 3
##  [71] 3 3 3 3 3 2 2 2 3 3 3 3 3 3 3 3 2 3 3 3 3 3 3 3 3 3 3 3 3 3 2 3 2 2 2
## [106] 2 3 2 2 2 2 2 2 3 2 2 2 2 2 3 2 3 2 3 2 2 3 3 2 2 2 2 2 3 3 2 2 2 3 2
## [141] 2 2 3 2 2 2 3 2 2 3

# Visualize pam clusters
fviz_cluster(pam.res, stand = FALSE, geom = "point",
             frame.type = "norm")

Read more about partitioning methods: Partitioning clustering

4 Example of hierarchical clustering results

The built-in R function hclust() is used:

# Compute pairewise distance matrices
dist.res <- dist(iris.scaled, method = "euclidean")
# Hierarchical clustering results
hc <- hclust(dist.res, method = "complete")
# Visualization of hclust
plot(hc, labels = FALSE, hang = -1)
# Add rectangle around 3 groups
rect.hclust(hc, k = 3, border = 2:4) 

# Cut into 3 groups
hc.cut <- cutree(hc, k = 3)
head(hc.cut, 20)

##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

Read more about hierarchical clustering: Hierarchical clustering

5 Three popular methods for determining the optimal number of clusters

In this section we describe the three most popular methods including: i) Elbow method, ii) silhouette method and iii) gap statistic.

set.seed(123)
# Compute and plot wss for k = 2 to k = 15
k.max <- 15 # Maximal number of clusters
data <- iris.scaled
wss <- sapply(1:k.max, 
        function(k){kmeans(data, k, nstart=10 )$tot.withinss})

plot(1:k.max, wss,
       type="b", pch = 19, frame = FALSE, 
       xlab="Number of clusters K",
       ylab="Total within-clusters sum of squares")
abline(v = 3, lty =2)

fviz_nbclust(x, FUNcluster, method = c("silhouette", "wss"))

fviz_nbclust(iris.scaled, kmeans, method = "wss") +
    geom_vline(xintercept = 3, linetype = 2)

fviz_nbclust(iris.scaled, pam, method = "wss") +
  geom_vline(xintercept = 3, linetype = 2)

fviz_nbclust(iris.scaled, hcut, method = "wss") +
  geom_vline(xintercept = 3, linetype = 2)

library(cluster)
k.max <- 15
data <- iris.scaled
sil <- rep(0, k.max)

# Compute the average silhouette width for 
# k = 2 to k = 15
for(i in 2:k.max){
  km.res <- kmeans(data, centers = i, nstart = 25)
  ss <- silhouette(km.res$cluster, dist(data))
  sil[i] <- mean(ss[, 3])
}

# Plot the  average silhouette width
plot(1:k.max, sil, type = "b", pch = 19, 
     frame = FALSE, xlab = "Number of clusters k")
abline(v = which.max(sil), lty = 2)

require(cluster)
fviz_nbclust(iris.scaled, kmeans, method = "silhouette")

require(cluster)
fviz_nbclust(iris.scaled, pam, method = "silhouette")

require(cluster)
fviz_nbclust(iris.scaled, hcut, method = "silhouette",
             hc_method = "complete")

clusGap(x, FUNcluster, K.max, B = 100, verbose = TRUE, ...)

# Compute gap statistic
library(cluster)
set.seed(123)
gap_stat <- clusGap(iris.scaled, FUN = kmeans, nstart = 25,
                    K.max = 10, B = 50)

# Print the result
print(gap_stat, method = "firstmax")

## Clustering Gap statistic ["clusGap"].
## B=50 simulated reference sets, k = 1..10
##  --> Number of clusters (method 'firstmax'): 3
##           logW   E.logW       gap     SE.sim
##  [1,] 4.534565 4.754595 0.2200304 0.02504585
##  [2,] 4.021316 4.489687 0.4683711 0.02742112
##  [3,] 3.806577 4.295715 0.4891381 0.02384746
##  [4,] 3.699263 4.143675 0.4444115 0.02093871
##  [5,] 3.589284 4.052262 0.4629781 0.02036366
##  [6,] 3.519726 3.972254 0.4525278 0.02049566
##  [7,] 3.448288 3.905945 0.4576568 0.02106987
##  [8,] 3.398210 3.850807 0.4525967 0.01969193
##  [9,] 3.334279 3.802315 0.4680368 0.01905974
## [10,] 3.250246 3.759661 0.5094149 0.01928183

# Base plot of gap statistic
plot(gap_stat, frame = FALSE, xlab = "Number of clusters k")
abline(v = 3, lty = 2)

# Use factoextra
fviz_gap_stat(gap_stat)

# Print
print(gap_stat, method = "Tibs2001SEmax")
# Plot
fviz_gap_stat(gap_stat, 
              maxSE = list(method = "Tibs2001SEmax"))
# Relaxed the gap test to be within two standard deviations
fviz_gap_stat(gap_stat, 
          maxSE = list(method = "Tibs2001SEmax", SE.factor = 2))

# Compute gap statistic
set.seed(123)
gap_stat <- clusGap(iris.scaled, FUN = pam, K.max = 10, B = 50)
# Plot gap statistic
fviz_gap_stat(gap_stat)

# Compute gap statistic
set.seed(123)
gap_stat <- clusGap(iris.scaled, FUN = hcut, K.max = 10, B = 50)
# Plot gap statistic
fviz_gap_stat(gap_stat)

6 NbClust: A Package providing 30 indices for determining the best number of clusters

NbClust(data = NULL, diss = NULL, distance = "euclidean",
        min.nc = 2, max.nc = 15, method = NULL, index = "all")

library("NbClust")
set.seed(123)
res.nb <- NbClust(iris.scaled, distance = "euclidean",
                  min.nc = 2, max.nc = 10, 
                  method = "complete", index ="gap") 
res.nb # print the results

## $All.index
##       2       3       4       5       6       7       8       9      10 
## -0.2899 -0.2303 -0.6915 -0.8606 -1.0506 -1.3223 -1.3303 -1.4759 -1.5551 
## 
## $All.CriticalValues
##       2       3       4       5       6       7       8       9      10 
## -0.0539  0.4694  0.1787  0.2009  0.2848  0.0230  0.1631  0.0988  0.1708 
## 
## $Best.nc
## Number_clusters     Value_Index 
##          3.0000         -0.2303 
## 
## $Best.partition
##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [36] 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 3 3 3 2 3 2 3 2 3 2 2 3 2 3 3 3 3 2 2 2
##  [71] 3 3 3 3 3 3 3 3 3 2 2 2 2 3 3 3 3 2 3 2 2 3 2 2 2 3 3 3 2 2 3 3 3 3 3
## [106] 3 2 3 3 3 3 3 3 3 3 3 3 3 3 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [141] 3 3 3 3 3 3 3 3 3 3

# All gap statistic values
res.nb$All.index

# Best number of clusters
res.nb$Best.nc

# Best partition
res.nb$Best.partition

nb <- NbClust(iris.scaled, distance = "euclidean", min.nc = 2,
        max.nc = 10, method = "complete", index ="all")

# Print the result
nb

fviz_nbclust(nb) + theme_minimal()

## Among all indices: 
## ===================
## * 2 proposed  0 as the best number of clusters
## * 1 proposed  1 as the best number of clusters
## * 2 proposed  2 as the best number of clusters
## * 18 proposed  3 as the best number of clusters
## * 3 proposed  10 as the best number of clusters
## 
## Conclusion
## =========================
## * Accoridng to the majority rule, the best number of clusters is  3 .

7 Infos

This analysis has been performed using R software (ver. 3.2.1)

• Charrad M., Ghazzali N., Boiteau V., Niknafs A. (2014). NbClust: An R Package for Determining the Relevant Number of Clusters in a Data Set. Journal of Statistical Software, 61(6), 1-36.
• Kaufman, L. and Rousseeuw, P.J. (1990). Finding Groups in Data: An Introduction to Cluster Analysis. Wiley, New York.
• Tibshirani, R., Walther, G. and Hastie, T. (2001). Estimating the number of data clusters via the Gap statistic. Journal of the Royal Statistical Society B, 63, 411–423. PDF