Description of pvclust() function

The function pvclust() can be used as follow:

pvclust(data, method.hclust = "average",
        method.dist = "correlation", nboot = 1000)

Note that, the computation time can be strongly decreased using parallel computation version called parPvclust(). (Read ?parPvclust() for more information.)

parPvclust(cl=NULL, data, method.hclust = "average",
           method.dist = "correlation", nboot = 1000,
           iseed = NULL)

The function pvclust() returns an object of class pvclust containing many elements including hclust which contains hierarchical clustering result for the original data generated by the function hclust().

Usage of pvclust() function

pvclust() performs clustering on the columns of the data set, which correspond to samples in our case. If you want to perform the clustering on the variables (here, genes) you have to transpose the data set using the function t().

The R code below computes pvclust() using 10 as the number of bootstrap replications (for speed):

library(pvclust)
set.seed(123)
res.pv <- pvclust(df, method.dist="cor", 
                  method.hclust="average", nboot = 10)

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

To extract the objects from the significant clusters, use the function pvpick():

clusters <- pvpick(res.pv)
clusters

Parallel computation can be applied as follow:

# Create a parallel socket cluster
library(parallel)
cl <- makeCluster(2, type = "PSOCK")
# parallel version of pvclust
res.pv <- parPvclust(cl, df, nboot=1000)
stopCluster(cl)

Silhouette coefficient

The silhouette analysis measures how well an observation is clustered and it estimates the average distance between clusters. The silhouette plot displays a measure of how close each point in one cluster is to points in the neighboring clusters.

For each observation \(i\), the silhouette width \(s_i\) is calculated as follows:

1. For each observation \(i\), calculate the average dissimilarity \(a_i\) between \(i\) and all other points of the cluster to which i belongs.

2. For all other clusters \(C\), to which i does not belong, calculate the average dissimilarity \(d(i, C)\) of \(i\) to all observations of C. The smallest of these \(d(i,C)\) is defined as \(b_i= \min_C d(i,C)\). The value of \(b_i\) can be seen as the dissimilarity between \(i\) and its “neighbor” cluster, i.e., the nearest one to which it does not belong.

3. Finally the silhouette width of the observation \(i\) is defined by the formula: \(S_i = (b_i - a_i)/max(a_i, b_i)\).

Silhouette width can be interpreted as follow:

Dunn index

The Dunn index is another internal clustering validation measure which can be computed as follow:

1. For each cluster, compute the distance between each of the objects in the cluster and the objects in the other clusters

2. Use the minimum of this pairwise distance as the inter-cluster separation (min.separation)

3. For each cluster, compute the distance between the objects in the same cluster.

4. Use the maximal intra-cluster distance (i.e maximum diameter) as the intra-cluster compactness

5. Calculate the Dunn index (D) as follow:

\[ D = \frac{min.separation}{max.diameter} \]

Required R packages

The following R packages are required in this chapter:

• factoextra for data visualization
• fpc for computing clustering validation statistics

• NbClust for determining the optimal number of clusters in the data set.

• Install the packages:

install.packages(c("factoextra", "fpc", "NbClust"))

• Load the packages:

library(factoextra)
library(fpc)
library(NbClust)

Data preparation

We’ll use the built-in R data set iris:

# Excluding the column "Species" at position 5
df <- iris[, -5]
# Standardize
df <- scale(df)

Clustering analysis

We’ll use the function eclust() [enhanced clustering, in factoextra] which provides several advantages:

• It simplifies the workflow of clustering analysis
• It can be used to compute hierarchical clustering and partitioning clustering in a single line function call
• Compared to the standard partitioning functions (kmeans, pam, clara and fanny) which requires the user to specify the optimal number of clusters, the function eclust() computes automatically the gap statistic for estimating the right number of clusters.
• It provides silhouette information for all partitioning methods and hierarchical clustering
• It draws beautiful graphs using ggplot2

The simplified format the eclust() function is as follow:

eclust(x, FUNcluster = "kmeans", hc_metric = "euclidean", ...)

The function eclust() returns an object of class eclust containing the result of the standard function used (e.g., kmeans, pam, hclust, agnes, diana, etc.).

It includes also:

• cluster: the cluster assignment of observations after cutting the tree
• nbclust: the number of clusters
• silinfo: the silhouette information of observations
• size: the size of clusters
• data: a matrix containing the original or the standardized data (if stand = TRUE)
• gap_stat: containing gap statistics

To compute a partitioning clustering, such as k-means clustering with k = 3, type this:

# K-means clustering
km.res <- eclust(df, "kmeans", k = 3, nstart = 25, graph = FALSE)
# Visualize k-means clusters
fviz_cluster(km.res, geom = "point", ellipse.type = "norm",
             palette = "jco", ggtheme = theme_minimal())

To compute a hierarchical clustering, use this:

# Hierarchical clustering
hc.res <- eclust(df, "hclust", k = 3, hc_metric = "euclidean", 
                 hc_method = "ward.D2", graph = FALSE)
# Visualize dendrograms
fviz_dend(hc.res, show_labels = FALSE,
         palette = "jco", as.ggplot = TRUE)

Cluster validation

fviz_silhouette(km.res, palette = "jco", 
                ggtheme = theme_classic())

##   cluster size ave.sil.width
## 1       1   50          0.64
## 2       2   47          0.35
## 3       3   53          0.39

# Silhouette information
silinfo <- km.res$silinfo
names(silinfo)
# Silhouette widths of each observation
head(silinfo$widths[, 1:3], 10)
# Average silhouette width of each cluster
silinfo$clus.avg.widths
# The total average (mean of all individual silhouette widths)
silinfo$avg.width
# The size of each clusters
km.res$size

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

##     cluster neighbor sil_width
## 112       2        3   -0.0106
## 128       2        3   -0.0249

cluster.stats(d = NULL, clustering, al.clustering = NULL)

library(fpc)
# Statistics for k-means clustering
km_stats <- cluster.stats(dist(df),  km.res$cluster)
# Dun index
km_stats$dunn

## [1] 0.0265

km_stats

External clustering validation

Among the values returned by the function cluster.stats(), there are two indexes to assess the similarity of two clustering, namely the corrected Rand index and Meila’s VI.

We know that the iris data contains exactly 3 groups of species.

Does the K-means clustering matches with the true structure of the data?

We can use the function cluster.stats() to answer to this question.

Let start by computing a cross-tabulation between k-means clusters and the reference Species labels:

table(iris$Species, km.res$cluster)

##             
##               1  2  3
##   setosa     50  0  0
##   versicolor  0 11 39
##   virginica   0 36 14

It can be seen that:

• All setosa species (n = 50) has been classified in cluster 1
• A large number of versicor species (n = 39 ) has been classified in cluster 3. Some of them ( n = 11) have been classified in cluster 2.
• A large number of virginica species (n = 36 ) has been classified in cluster 2. Some of them (n = 14) have been classified in cluster 3.

It’s possible to quantify the agreement between Species and k-means clusters using either the corrected Rand index and Meila’s VI provided as follow:

library("fpc")
# Compute cluster stats
species <- as.numeric(iris$Species)
clust_stats <- cluster.stats(d = dist(df), 
                             species, km.res$cluster)
# Corrected Rand index
clust_stats$corrected.rand

## [1] 0.62

# VI
clust_stats$vi

## [1] 0.748

The same analysis can be computed for both PAM and hierarchical clustering:

# Agreement between species and pam clusters
pam.res <- eclust(df, "pam", k = 3, graph = FALSE)
table(iris$Species, pam.res$cluster)
cluster.stats(d = dist(iris.scaled), 
              species, pam.res$cluster)$vi
# Agreement between species and HC clusters
res.hc <- eclust(df, "hclust", k = 3, graph = FALSE)
table(iris$Species, res.hc$cluster)
cluster.stats(d = dist(iris.scaled), 
              species, res.hc$cluster)$vi

Required R packages

We’ll use the following R packages:

• factoextra to determine the optimal number clusters for a given clustering methods and for data visualization.
• NbClust for computing about 30 methods at once, in order to find the optimal number of clusters.

To install the packages, type this:

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

Load the packages as follow:

library(factoextra)
library(NbClust)

Data preparation

We’ll use the USArrests data as a demo data set. We start by standardizing the data to make variables comparable.

# Standardize the data
df <- scale(USArrests)
head(df)

##            Murder Assault UrbanPop     Rape
## Alabama    1.2426   0.783   -0.521 -0.00342
## Alaska     0.5079   1.107   -1.212  2.48420
## Arizona    0.0716   1.479    0.999  1.04288
## Arkansas   0.2323   0.231   -1.074 -0.18492
## California 0.2783   1.263    1.759  2.06782
## Colorado   0.0257   0.399    0.861  1.86497

fviz_nbclust() function: Elbow, Silhouhette and Gap statistic methods

The simplified format is as follow:

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

The R code below determine the optimal number of clusters for k-means clustering:

# Elbow method
fviz_nbclust(df, kmeans, method = "wss") +
    geom_vline(xintercept = 4, linetype = 2)+
  labs(subtitle = "Elbow method")
# Silhouette method
fviz_nbclust(df, kmeans, method = "silhouette")+
  labs(subtitle = "Silhouette method")
# Gap statistic
# nboot = 50 to keep the function speedy. 
# recommended value: nboot= 500 for your analysis.
# Use verbose = FALSE to hide computing progression.
set.seed(123)
fviz_nbclust(df, kmeans, nstart = 25,  method = "gap_stat", nboot = 50)+
  labs(subtitle = "Gap statistic method")

## Clustering k = 1,2,..., K.max (= 10): .. done
## Bootstrapping, b = 1,2,..., B (= 50)  [one "." per sample]:
## .................................................. 50

According to these observations, it’s possible to define k = 4 as the optimal number of clusters in the data.

NbClust() function: 30 indices for choosing the best number of clusters

The simplified format of the function NbClust() is:

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

• To compute NbClust() for kmeans, use method = “kmeans”.
• To compute NbClust() for hierarchical clustering, method should be one of c(“ward.D”, “ward.D2”, “single”, “complete”, “average”).

The R code below computes NbClust() for k-means:

library("NbClust")
nb <- NbClust(df, distance = "euclidean", min.nc = 2,
        max.nc = 10, method = "kmeans")

The result of NbClust using the function fviz_nbclust() [in factoextra], as follow:

library("factoextra")
fviz_nbclust(nb)

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

Statistical methods

The Hopkins statistic (Lawson and Jurs 1990) is used to assess the clustering tendency of a data set by measuring the probability that a given data set is generated by a uniform data distribution. In other words, it tests the spatial randomness of the data.

For example, let D be a real data set. The Hopkins statistic can be calculated as follow:

1. Sample uniformly \(n\) points (\(p_1\),…, \(p_n\)) from D.

2. Compute the distance, \(x_i\), from each real point to each nearest neighbor: For each point \(p_i \in D\), find it’s nearest neighbor \(p_j\); then compute the distance between \(p_i\) and \(p_j\) and denote it as \(x_i = dist(p_i, p_j)\)

3. Generate a simulated data set (\(random_D\)) drawn from a random uniform distribution with \(n\) points (\(q_1\),…, \(q_n\)) and the same variation as the original real data set D.

4. Compute the distance, \(y_i\) from each artificial point to the nearest real data point: For each point \(q_i \in random_D\), find it’s nearest neighbor \(q_j\) in D; then compute the distance between \(q_i\) and \(q_j\) and denote it \(y_i = dist(q_i, q_j)\)

5. Calculate the Hopkins statistic (H) as the mean nearest neighbor distance in the random data set divided by the sum of the mean nearest neighbor distances in the real and across the simulated data set.

The formula is defined as follow:

\[H = \frac{\sum\limits_{i=1}^ny_i}{\sum\limits_{i=1}^nx_i + \sum\limits_{i=1}^ny_i}\]

How to interpret the Hopkins statistics? If \(D\) were uniformly distributed, then \(\sum\limits_{i=1}^ny_i\) and \(\sum\limits_{i=1}^nx_i\) would be close to each other, and thus \(H\) would be about 0.5. However, if clusters are present in D, then the distances for artificial points (\(\sum\limits_{i=1}^ny_i\)) would be substantially larger than for the real ones (\(\sum\limits_{i=1}^nx_i\)) in expectation, and thus the value of \(H\) will increase (Han, Kamber, and Pei 2012).

A value for \(H\) higher than 0.75 indicates a clustering tendency at the 90% confidence level.

The null and the alternative hypotheses are defined as follow:

• Null hypothesis: the data set D is uniformly distributed (i.e., no meaningful clusters)
• Alternative hypothesis: the data set D is not uniformly distributed (i.e., contains meaningful clusters)

Here, we present two R functions / packages to statistically evaluate clustering tendency by computing the Hopkins statistics:

1. get_clust_tendency() function [in factoextra package]. It returns the Hopkins statistics as defined in the formula above. The result is a list containing two elements:
   • hopkins_stat
   • and plot
2. hopkins() function [in clustertend package]. It implements 1- the definition of H provided here.

In the R code below, we’ll use the factoextra R package. Make sure that you have the latest version (or install: devtools::install_github("kassambara/factoextra")).

library(factoextra)
# Compute Hopkins statistic for iris dataset
res <- get_clust_tendency(df, n = nrow(df)-1, graph = FALSE)
res$hopkins_stat

## [1] 0.818

# Compute Hopkins statistic for a random dataset
res <- get_clust_tendency(random_df, n = nrow(random_df)-1,
                          graph = FALSE)
res$hopkins_stat

## [1] 0.466

Visual methods

The algorithm of the visual assessment of cluster tendency (VAT) approach (Bezdek and Hathaway, 2002) is as follow:

The algorithm of VAT is as follow:

For the visual assessment of clustering tendency, we start by computing the dissimilarity matrix between observations using the function dist(). Next the function fviz_dist() [factoextra package] is used to display the dissimilarity matrix.

fviz_dist(dist(df), show_labels = FALSE)+
  labs(title = "Iris data")

fviz_dist(dist(random_df), show_labels = FALSE)+
  labs(title = "Random data")

• Red: high similarity (ie: low dissimilarity) | Blue: low similarity

The color level is proportional to the value of the dissimilarity between observations: pure red if \(dist(x_i, x_j) = 0\) and pure blue if \(dist(x_i, x_j) = 1\). Objects belonging to the same cluster are displayed in consecutive order.

The VAT detects the clustering tendency in a visual form by counting the number of square shaped dark blocks along the diagonal in a VAT image.