2 Clustering

2.1 Motivation

The image above is from a microarray experiment ¹. The intensity of each dot represents gene expression of a single gene (how “active” the gene is) for a particular individual/sample. The resultant data is therefore a big matrix of numbers, where each column represents a gene and each row an individual/sample.

There are two questions of interest:

1. Which genes are co-regulated, that is, behave in the same way?
2. Which individuals are similar to each other, that is, have a similar gene expression profile?

In both cases we want to discover some underlying structure in unlabelled data. Structure means patterns in the data that are sufficiently different from pure unstructured noise. Here we introduce clustering, a class of unsupervised learning methods that try to answer these questions.

2.2 What is clustering?

The goal of clustering is to find groups that share similar properties. The data in each group should be similar (minimise intracluster distance), but each cluster should be sufficiently different (maximise intercluster similarity).

2.3 What problems can clustering solve?

Clustering is particularly useful in applications where labelling the data is very time consuming/expensive.

• Gene expression: discovering co-regulated genes.

• Biological systematics: finding organisms sharing similar attributes.

• Computer vision: segmenting a digital image for object recognition.

• Epidemiology: identifying geographical clusters of diseases.

• Medical imaging: differentiating between tissues.

• Mathematical chemistry: grouping compounds by topological indices.

• Market basket analysis: determining which group of items tend to be bought together.

• Cybersecurity: detecting fraudulent activity.

• … and much more!

2.4 Types of clustering methods

• Partitional: the feature space is partitioned into \(k\) regions e.g \(k\)-means.

• Hierarchical: iteratively merging small clusters into larger ones (agglomerative) or breaking large clusters into smaller ones (divisive).

• Distribution-based: fit \(k\) multivariate statistical distributions e.g Gaussian mixture model (GMM).

2.5 Similarity measures

Most clustering methods rely on distance metrics that quantify how close two observations are. There are several ways to define this distance, which has a direct effect on the clustering result.

The Euclidean distance (think Pythagoras theorem) is depicted below, together with the Manhatttan distance (named after the journey a taxi has to follow in grid-like streets of cities like Manhattan).

The correlation coefficient is also another popular way to measure similarity.

There are various other distance metrics, please see dist in R or pdist in Python. In this introductory workshop we will focus on continuous features, but be aware that distance measures for categorical variables exists, such as, the Jaccard index, Gower distance and polychoric correlation.

2.6 The Iris dataset

To showcase some of the clustering methods, we will use the popular Iris dataset. The data set consists of 50 samples from three species of Iris flower (I. setosa, I. virginica and I. versicolor). Each flower is quantified by four measurements, length and width of sepal and petal.

Let us load this dataset:

# The iris dataset is preloaded in R
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

# The iris dataset is available from the sci-kit learn package
from sklearn import datasets
iris = datasets.load_iris()

# Print the first 6 rows
# Sepal Length, Sepal Width, Petal Length, Petal Width
iris.data[:6, ]

array([[5.1, 3.5, 1.4, 0.2],
       [4.9, 3. , 1.4, 0.2],
       [4.7, 3.2, 1.3, 0.2],
       [4.6, 3.1, 1.5, 0.2],
       [5. , 3.6, 1.4, 0.2],
       [5.4, 3.9, 1.7, 0.4]])

2.7 \(k\)-means

Arguably the most widely used partitioning clustering method. The feature space is divided into \(k\) regions as follows:

1. Select \(k\) centroids at random.

2. Compute the Euclidean distance between centroids and each data point.

3. Assign each data point to the closest centroid.

4. Compute new centroids; the average of all data points in that cluster.

5. Repeat steps 2 to 4 until data points remain in the same cluster or some maximum number of iterations reached.

Note: \(k\)-means clustering should only be used with continuous data!

For visualisation purposes let’s just use two features of the Iris dataset; sepal length and petal width.

# Fit k-means model
k <- 3
mdl <- kmeans(x=iris[, c(1, 4)], centers=k)

# Associate a colour with each cluster
library(RColorBrewer)
COL <- seq(k)
names(COL) <- brewer.pal(n=k, 'Set1')

# Plot results
plot(iris[, 4], iris[, 1], pch=19, col=names(COL[mdl$cluster]),
     xlab='Petal width (cm)', ylab='Sepal length (cm)')

# Fit k-means model
from sklearn.cluster import KMeans
k = 3
mdl = KMeans(n_clusters=k, n_jobs=-1) # -1 uses all cores
mdl.fit(X=iris.data[:, [0, 3]])

# Associate a colour with each cluster

from palettable.colorbrewer.qualitative import Set1_3
colDict = {0: Set1_3.hex_colors[0], 1: Set1_3.hex_colors[1], 2: Set1_3.hex_colors[2]}
myCol = [colDict[i] for i in mdl.labels_]

# Plot results
import matplotlib.pyplot as plt
plt.scatter(iris.data[:, 3], iris.data[:, 0], c=myCol)
plt.xlabel('Petal width (cm)')
plt.ylabel('Sepal length (cm)')

Pros	Cons
Simple and intuitive	\(k\) needs to be specified a priori
Computationally inexpensive/fast	Only applicable for continuous data where a mean is defined
	No guarantee of a global optimum solution

2.8 Agglomerative hiearchical clustering

In agglomerative hierarchical clustering small clusters are iteratively merged into larger ones. The clustering strategy is as follows:

1. Assign each datum as its own cluster.

2. Compute the distance between each cluster.

3. Merge the closest pair into a single cluster.

4. Repeat steps 2 to 3 until all clusters are merged together.

Step 3 is key, the distance metric and linkage function dictate the final result. The linkage function specifies how the inter-cluster distance is computed. There are various options:

• Centroid: mean of data points (same as in \(k\)-means).

• Single: distance between closest pair of points.

• Complete: distance between furthest pair of points.

• Average: mean pairwise distance between all points.

The distance can be computed using any similarity measure introduced previously.

# Compute distance matrix
d <- dist(x=iris[, c(1, 4)], method='euclidean') 

# Perform agglomerative hierarchical clustering
# Use 'average' link function
mdl <- hclust(d=d, method='average')

# Plot resultant dendrogram
plot(mdl, cex=0.6)

from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import linkage, dendrogram

# Compute distance matrix
d = pdist(X=iris.data[:, [0, 3]], metric="euclidean")

# Perform agglomerative hierarchical clustering
# Use 'average' link function
mdl = linkage(d, method='average')

# Plot resultant dendrogram
plt.figure(figsize=(9,6))
dendrogram(mdl)

plt.show()

The number at the end of each branch corresponds to the observation row number.

Pros	Cons
No need to specify \(k\)	Can be computationally expensive
Sub-groups within larger clusters can be easily identified	Interpretation is subjective. Where should we draw the line?
Dendrograms let us visualise results irrespective of number of features	Choice of distance method and linkage function can significantly change the result

2.9 Gaussian mixture model (GMM)

The GMM is a simple but powerful model that performs clustering via density estimation. The features’ histogram is modelled as the sum of multiple multivariate Gaussian distributions. Suppose we only had access to one feature, a GMM with \(k=2\) would look something like this:

The blue dashed lines represent the two individual univariate Gaussians, whilst the black line depicts the combined model. We can extend this to more features by using multivariate Gaussians. Mathematically this can be expressed as follows:

\[ p(x) = \sum_{i=1}^k \pi_i \mathcal{N}(x|\mu_i, \Sigma_i)\\ \sum_{i=1}^k \pi_i = 1 \] The Expectation-Maximisation (EM) algorithm is used to estimate the parameters \(\pi_i\) (known as mixing coefficients), \(\mu_i\) and \(\Sigma_i\).

library(mclust)

Package 'mclust' version 5.4.5
Type 'citation("mclust")' for citing this R package in publications.

# Fit Gaussian Mixture Model
k <- 3 # no. of clusters
mdl <- Mclust(data=iris[, c(4, 1)], G=3)

# Plot results
plot(mdl, what='classification',
     xlab='Petal width (cm)', 
     ylab='Sepal length (cm)')

import numpy as np
from matplotlib.colors import LogNorm
from sklearn.mixture import GaussianMixture as GMM

# Fit Gaussian Mixture Model
k = 3 # no. of clusters
mdl = GMM(n_components=3)
mdl.fit(X=iris.data[:, [3, 0]])

# Compute probability distribution function at each point on a gird

x = np.linspace(np.min(iris.data[:, 3]), np.max(iris.data[:, 3]), 100)
y = np.linspace(np.min(iris.data[:, 0]), np.max(iris.data[:, 0]), 100)
X, Y = np.meshgrid(x, y)
XX = np.array([X.ravel(), Y.ravel()]).T
Z = -mdl.score_samples(XX)
Z = Z.reshape(X.shape)

# Plot results
hPlot = plt.contour(X, Y, Z, norm=LogNorm(), 
                 levels=np.logspace(0, 3, 10))

plt.colorbar(hPlot, shrink=0.8, extend='both')

plt.scatter(iris.data[:, 3], iris.data[:, 0])
plt.xlabel('Petal width (cm)')
plt.ylabel('Sepal length (cm)')

Pros	Cons
Intuitive interpretation	\(k\) needs to be specified a priori
Computationally inexpensive	Strong assumption on the distribution of the feature space (multivariate Gaussian)
	No guarantee of a global optimum solution
	Fails when number of features is much greater than observations

GMMs offer a “soft” clustering approach, where every observation is part of every cluster but with varying levels of membership.

2.10 Determining the “correct” number of clusters

One of the biggest questions when it comes to clustering is “How many clusters do I have?”. The number of clusters \(k\) cannot be determined exactly, because the observations are unlabelled, so \(k\) is inherently ambiguous. Moreover, similarity is quite subjective and often we cannot define a clear cut-off.

For example, suppose that as part of a public health exercise we want to cluster a large group of individuals based on their health. Health is a multifaceted concept and cannot be observed directly; instead we measure various biomarkers, like body mass index (BMI), cholesterol levels, body composition, resting metabolic rate, etc. Although we would be able to differentiate between individuals at the two extremes (i.e athelete vs couch potato), most people will sit somewhere on a continuum. There isn’t a clear “line”, that once crossed an individual goes from being healthy to a bit unhealthy or moderately unhealthy etc. The number of clusters is therefore somewhat dictated by the problem at hand and the type of questions we’re trying to answer.

Nevertheless, there are various metrics that one can use to estimate the underlying number of clusters:

• Recall that the objective of clustering is to minimise the intracluster distance and maximise the intercluster similarity. Thus, we can plot the within and between clusters sum-of-squares distances as a function of \(k\). As we increase the number of clusters, there will be a point where the sum-of-squares distances will only change marginally, that is, adding more clusters does not improve these metrics significantly. The number of clusters is chosen to be the point at which the curve “plateaus” (5 in the synthetic example below). This is known as the “elbow criterion”. Please refer to the R or Python documentation on how to access these metrics.

• The silhouette width quantifies how similar an observation is to its own cluster compared to other clusters. This measure ranges from -1 (not compatible with that cluster) to 1 (extremely likely to be part of that cluster). The suggested configuration is the one that maximises the average silhouette width (3 in the synthetic example below). Please refer to the R or Python documentation on how to compute and plot these metrics.

• For distribution-based methods, choosing \(k\) can be framed as a model selection problem. We can plot the Akaike Information Criterion (AIC), Bayesian Information Criterior (BIC) or other information criterion measures. As we increase the number of clusters, there will be a point where the model fit will only improve marginally or start to decrease. The number of clusters is chosen to be the point at which the curve “plateaus” (the “elbow criterion”; 5 in the synthetic example below). Please refer to the R or Python documentation on how to compute these metrics.

There are myriad other metrics available in the literature, some related to specific clustering algorithms. You will also encounter methods that claim to automatically discover the optimal number of clusters for you. Although, this can be true in a mathematical sense, this estimate will still be based on various underlying assumptions and hyperparameters.

These cluster validity measures only give us a ballpark range for the “correct” number of clusters. Ultimately one needs to make use of prior knowledge to determine whether the number of clusters are practically relevant and if they make sense. For example, how many different phenotypes are you expecting in your population?

2.11 Tasks

2.11.1 Simulated data

Let’s start to get a feel for these clustering algorithms by simulating some data:

library(MASS) # mvrnorm (multivariate normal)

# Set simulation parameters
N <- 50 # no. of data points in each cluster
covMatrix <- matrix(data=c(1, 0, 0, 2), nrow=2)

# Simulate clusters (assume same covariances for now)
set.seed(1034) # to reproduce results
clustA <- mvrnorm(n=N, mu=c(6, 4), Sigma=covMatrix)
clustB <- mvrnorm(n=N, mu=c(3, 9), Sigma=covMatrix)
clustC <- mvrnorm(n=N, mu=c(9, 9), Sigma=covMatrix)

# Join all the data together and plot
xTrain <- rbind(clustA, clustB, clustC)
plot(xTrain[, 1], xTrain[, 2], pch=19,
xlab='Feature 1', ylab='Feature 2')

from numpy.random import multivariate_normal, seed

# Set simulation parameters
N = 50 # no. of data points in each cluster
covMatrix = np.array([[1, 0], [0, 2]], dtype='float')

# Simulate clusters (assume same covariances for now)
seed(1034) # to reproduce results
clustA = multivariate_normal(mean=np.array([6, 4]), cov=covMatrix, size=N)
clustB = multivariate_normal(mean=np.array([3, 9]), cov=covMatrix, size=N)
clustC = multivariate_normal(mean=np.array([9, 9]), cov=covMatrix, size=N)

# Join all the data together and plot
xTrain = np.vstack((clustA, clustB, clustC))
plt.figure(figsize=(6,6))
plt.scatter(xTrain[:, 0], xTrain[:, 1])
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')

Task

• Perform \(k\)-means clustering by allowing \(k\) to vary from 2 to 6 .

• Plot the intra and intercluster sum-of-squares as a function of \(k\) and deduce the “true” number of underlying clusters.

Solution

# Initialise some variables
kRange <- seq(from=2, to=6, by=1)
intra <- rep(NA, length(kRange)) # intracluster sum-of-squares
inter <- rep(NA, length(kRange)) # intercluster sum-of-squares

# Loop across desired range of ks
for (k in kRange) {
    mdl <- kmeans(x=xTrain , centers=k)
    intra[k-1] <- mdl$tot.withinss # it’s (k-1) because k starts from 2
    inter[k-1] <- mdl$betweenss 
}

# Plot inter/intercluster sum-of-squares as a function of $k$
yMin <- min(c(intra , inter)) 
yMax <- max(c(intra , inter))
plot(kRange, inter, type="o", pch=1, col="blue", lwd=3, lty=1, 
ylim=c(yMin, yMax), xlab="k", ylab='Sum-of-squares')
points(kRange, intra, type="o", pch=1, col="red", lwd=3)
legend("right", c("inter-cluster", "intra-cluster"), 
bty="n", col=c("blue", "red"), pch=1, lwd=3, lty=1)

The “correct” number of clusters is approximately 3. This matches our expectation.

from sklearn.preprocessing import scale

# Initialise some variables
kRange = range(2, 7)
intra = np.empty(len(kRange)) # intracluster sum-of-squares
inter = np.empty(len(kRange)) # intercluster sum-of-squares

# Loop across desired range of ks
for k in kRange:
    mdl = KMeans(n_clusters=k)
    mdl.fit(X=xTrain)
    intra[k-2] = mdl.inertia_  # it’s (k-2) because k starts from 2
    inter[k-2] = np.sum(scale(xTrain, with_std=False)**2) - mdl.inertia_

# Plot inter/intercluster sum-of-squares as a function of $k$

plt.figure(figsize=(5,5))
plt.plot(kRange, inter, 'bo-', label='inter-cluster')
plt.plot(kRange, intra, 'ro-', label='intra-cluster')
plt.legend(loc='center right')
plt.xlabel('k')
plt.ylabel('Sum-of-squares')

The “correct” number of clusters is approximately 3. This matches our expectation.

Task

• Fit a Gaussian Mixture Model (GMM) by allowing \(k\) to vary from 2 to 6 .

• Plot the AIC and BIC as a function of \(k\) and deduce the “true” number of underlying clusters.

Solution

# Initialise some variables
kRange <- seq(from=2, to=6, by=1)
AIC <- rep(NA, length(kRange)) # Akaike information criterion
BIC <- rep(NA, length(kRange)) # Bayesian information criterion


# Loop across desired range of ks
for (k in kRange) {
    mdl <- Mclust(data=xTrain, G=k)
    AIC[k-1] <- 2*mdl$loglik - 2*mdl$df
    BIC[k-1] <- mdl$bic
}

# Plot inter/intercluster sum-of-squares as a function of $k$
yMin <- min(c(AIC , BIC)) 
yMax <- max(c(AIC , BIC))
plot(kRange, AIC, type="o", pch=1, col="blue", lwd=3, lty=1, 
ylim=c(yMin, yMax), xlab="k", ylab='Information criterion')
points(kRange, BIC, type="o", pch=1, col="red", lwd=3)
legend("right", c("AIC", "BIC"), 
bty="n", col=c("blue", "red"), pch=1, lwd=3, lty=1)

The “correct” number of clusters is approximately 3. This matches our expectation.

# Initialise some variables
kRange = range(2, 7)
AIC = np.empty(len(kRange)) # Akaike information criterion
BIC = np.empty(len(kRange)) # Bayesian information criterion

# Loop across desired range of ks
for k in kRange:
    mdl = GMM(n_components=k)
    mdl.fit(X=xTrain)
    AIC[k-2] = -mdl.aic(xTrain)
    BIC[k-2] = -mdl.bic(xTrain)

# Plot inter/intercluster sum-of-squares as a function of $k$

plt.figure(figsize=(6,6))
plt.plot(kRange, AIC, 'bo-', label='AIC')
plt.plot(kRange, BIC, 'ro-', label='BIC')
plt.legend(loc='center right')
plt.xlabel('k')
plt.ylabel('Information criterion')

The “correct” number of clusters is approximately 3. This matches our expectation.

Task

Perform the above analysis on several sets of simulated data by changing the mean and covariance matrix of each simulated cluster. Try bringing the clusters closer together and then push them further apart. What happens to intra/intercluster distance and AIC/BIC plots?

Solution

• As clusters move closer together, the “elbow” becomes less pronounced due to the large overlap between clusters. The “correct” number of clusters will inevitably become more ambiguous.

• As clusters move further apart, the “elbow” becomes more pronounced, making it easier to determine the “correct” number of clusters. Note, however, that we could still have sub-groups within a particular cluster; it’s just that the within-cluster dissimilarity would be subtle compared to the across-clusters dissimilarity. For example, if one were to cluster a large set of images containing dogs, cats and horses, you would expect to find three distinct clusters. However, we know that within the dogs cluster you’d also have sub-clusters related to different species of dogs.

2.11.2 Gene expression

The file gene_expression.csv (all workshop datasets are available here), contains the acute lymphoblastic leukaemia (ALL) dataset which was published in the following study (see here and here). The dataset contains normalised gene expression values (measured using microarray) for 128 patients and 12,625 genes. The patients were diagnosed with either a B- or T-cell acute lymphocytic leukaemia.

Do not worry too much about the details (i.e what the genes are etc.), treat this dataset as a \(G \times N\) matrix where \(G\) is the total number of genes and \(N\) is the number of patients. We have access to the labels, type and stage of the disease (e.g B2). Thus, we can easily assess how well the clustering algorithm is doing, as we expect the B’s and T’s to cluster together.

xTrain <- read.csv('_data/gene_expression.csv', row.names=1)
print(dim(xTrain))

[1] 12625   128

import pandas as pd

xTrain = pd.read_csv('_data/gene_expression.csv', header=0, index_col=0)
print(xTrain.shape)

(12625, 128)

Task

Perform agglomerative hierarchical clustering on the 128 patients. Keep the distance method as euclidean, but change the linkage method (e.g single, average) and observe how the dendrogram changes.

Solution

linkMethods <- c('single', 'complete', 'average')
distMethod <- 'euclidean'
distance <- dist(as.matrix(t(xTrain), method=distMethod))
for (linkMethod in linkMethods) {
    mdl <- hclust(distance , method=linkMethod)
    plot(mdl, cex=0.5, xlab='distance', 
    main=paste0(distMethod, ' distance and ', linkMethod, ' link function'))
}

linkMethods = ['single', 'complete', 'average']
distMethod = 'euclidean'
distance = pdist(X=xTrain.T, metric=distMethod )
for linkMethod in linkMethods:
    mdl = linkage(distance, method=linkMethod)
    plt.figure(num=linkMethod, figsize=(9,6))
    dendrogram(mdl)
    plt.title('{0} distance and {1} link function'.format(distMethod, linkMethod))
    plt.show()

Task

Same as before, but now keep the linkage method as average and change the distance method (e.g euclidean, manhattan) and observe how the dendrogram changes.

Solution

linkMethod <- 'average'
distMethods <- c('euclidean', 'manhattan', 'canberra')
for (distMethod in distMethods) {
    distance <- dist(as.matrix(t(xTrain), method=distMethod))
    mdl <- hclust(distance , method=linkMethod)
    plot(mdl, cex=0.5, xlab='distance', 
    main=paste0(distMethod, ' distance and ', linkMethod, ' link function'))
}

linkMethod = 'average'
distMethods = ['euclidean', 'cityblock', 'canberra']
for distMethod in distMethods:
    distance = pdist(X=xTrain.T, metric=distMethod )
    mdl = linkage(distance, method=linkMethod)
    plt.figure(num=distMethod, figsize=(9,6))
    dendrogram(mdl)
    plt.title('{0} distance and {1} link function'.format(distMethod, linkMethod))
    plt.show()

2.11.3 Wine

The file wine.csv contains chemical analysis data of wines grown in the same region in Italy but from three different cultivars (see here for details).

xTrain <- read.csv('_data/wine.csv')
print(head(xTrain))

  WineType Alcohol MalicAcid  Ash AlcalinityAsh Magnesium TotalPhenols
1        A   14.23      1.71 2.43          15.6       127         2.80
2        A   13.20      1.78 2.14          11.2       100         2.65
3        A   13.16      2.36 2.67          18.6       101         2.80
4        A   14.37      1.95 2.50          16.8       113         3.85
5        A   13.24      2.59 2.87          21.0       118         2.80
6        A   14.20      1.76 2.45          15.2       112         3.27
  Flavanoids NonflavanoidPhenols Proanthocyanins ColorIntensity  Hue
1       3.06                0.28            2.29           5.64 1.04
2       2.76                0.26            1.28           4.38 1.05
3       3.24                0.30            2.81           5.68 1.03
4       3.49                0.24            2.18           7.80 0.86
5       2.69                0.39            1.82           4.32 1.04
6       3.39                0.34            1.97           6.75 1.05
  OD280_OD315 Proline
1        3.92    1065
2        3.40    1050
3        3.17    1185
4        3.45    1480
5        2.93     735
6        2.85    1450

xTrain = pd.read_csv('_data/wine.csv', header=0)
print(xTrain.head())

  WineType  Alcohol  MalicAcid  ...   Hue  OD280_OD315  Proline
0        A    14.23       1.71  ...  1.04         3.92     1065
1        A    13.20       1.78  ...  1.05         3.40     1050
2        A    13.16       2.36  ...  1.03         3.17     1185
3        A    14.37       1.95  ...  0.86         3.45     1480
4        A    13.24       2.59  ...  1.04         2.93      735

[5 rows x 14 columns]

There are thirteen variables (Alcohol, MalicAcid, etc.), together with WineType, which specifies the type of wine. Here, we are going to pretend we do not know that there are three types of wine, instead we’ll use clustering methods to uncover this information.

One thing to notice with this data, is that the units vary greatly across variables. For example, Proline ranges from 278 to 1680, whilst MalicAcid ranges from 0.74 to 5.8. So first we need to normalise the data, so that they’re on a common scale.

wineType <- xTrain$WineType # save for comparison
xTrain <- scale(xTrain[-1], center=TRUE, scale=TRUE)
print(head(xTrain))

       Alcohol   MalicAcid        Ash AlcalinityAsh  Magnesium TotalPhenols
[1,] 1.5143408 -0.56066822  0.2313998    -1.1663032 1.90852151    0.8067217
[2,] 0.2455968 -0.49800856 -0.8256672    -2.4838405 0.01809398    0.5670481
[3,] 0.1963252  0.02117152  1.1062139    -0.2679823 0.08810981    0.8067217
[4,] 1.6867914 -0.34583508  0.4865539    -0.8069748 0.92829983    2.4844372
[5,] 0.2948684  0.22705328  1.8352256     0.4506745 1.27837900    0.8067217
[6,] 1.4773871 -0.51591132  0.3043010    -1.2860793 0.85828399    1.5576991
     Flavanoids NonflavanoidPhenols Proanthocyanins ColorIntensity        Hue
[1,]  1.0319081          -0.6577078       1.2214385      0.2510088  0.3611585
[2,]  0.7315653          -0.8184106      -0.5431887     -0.2924962  0.4049085
[3,]  1.2121137          -0.4970050       2.1299594      0.2682629  0.3174085
[4,]  1.4623994          -0.9791134       1.0292513      1.1827317 -0.4263410
[5,]  0.6614853           0.2261576       0.4002753     -0.3183774  0.3611585
[6,]  1.3622851          -0.1755994       0.6623487      0.7298108  0.4049085
     OD280_OD315     Proline
[1,]   1.8427215  1.01015939
[2,]   1.1103172  0.96252635
[3,]   0.7863692  1.39122370
[4,]   1.1807407  2.32800680
[5,]   0.4483365 -0.03776747
[6,]   0.3356589  2.23274072

from sklearn.preprocessing import scale
wineType = xTrain['WineType'] # save for comparison
xTrain = scale(xTrain.drop(labels='WineType', axis=1), 
               with_mean=True, with_std=True)
print(xTrain)

[[ 1.51861254 -0.5622498   0.23205254 ...  0.36217728  1.84791957
   1.01300893]
 [ 0.24628963 -0.49941338 -0.82799632 ...  0.40605066  1.1134493
   0.96524152]
 [ 0.19687903  0.02123125  1.10933436 ...  0.31830389  0.78858745
   1.39514818]
 ...
 [ 0.33275817  1.74474449 -0.38935541 ... -1.61212515 -1.48544548
   0.28057537]
 [ 0.20923168  0.22769377  0.01273209 ... -1.56825176 -1.40069891
   0.29649784]
 [ 1.39508604  1.58316512  1.36520822 ... -1.52437837 -1.42894777
  -0.59516041]]

Task

• Perform \(k\)-means clustering by allowing \(k\) to vary from 1 to 10 .

• Plot the intra and intercluster sum-of-squares as a function of \(k\) and deduce the “true” number of underlying clusters.

Solution

# Initialise some variables
kRange <- seq(from=1, to=10, by=1)
intra <- rep(NA, length(kRange)) # intracluster sum-of-squares
inter <- rep(NA, length(kRange)) # intercluster sum-of-squares

# Loop across desired range of ks
for (k in kRange) {
    mdl <- kmeans(x=xTrain , centers=k)
    intra[k] <- mdl$tot.withinss
    inter[k] <- mdl$betweenss 
}

# Plot inter/intercluster sum-of-squares as a function of $k$
yMin <- min(c(intra , inter)) 
yMax <- max(c(intra , inter))
plot(kRange, inter, type="o", pch=1, col="blue", lwd=3, lty=1, 
     ylim=c(yMin, yMax), xlab="k", ylab='Sum-of-squares')
points(kRange, intra, type="o", pch=1, col="red", lwd=3)
legend("topright", c("inter-cluster", "intra-cluster"), 
       bty="n", col=c("blue", "red"), pch=1, lwd=3, lty=1)

The “correct” number of clusters is approximately 3. This matches our expectation.

from sklearn.preprocessing import scale

# Initialise some variables
kRange = range(1, 11)
intra = np.empty(len(kRange)) # intracluster sum-of-squares
inter = np.empty(len(kRange)) # intercluster sum-of-squares

# Loop across desired range of ks
for k in kRange:
    mdl = KMeans(n_clusters=k)
    mdl.fit(X=xTrain)
    intra[k-1] = mdl.inertia_ 
    inter[k-1] = np.sum(scale(xTrain, with_std=False)**2) - mdl.inertia_

# Plot inter/intercluster sum-of-squares as a function of $k$

plt.figure(figsize=(5,5))
plt.plot(kRange, inter, 'bo-', label='inter-cluster')
plt.plot(kRange, intra, 'ro-', label='intra-cluster')
plt.legend(loc='upper right')
plt.xlabel('k')
plt.ylabel('Sum-of-squares')

The “correct” number of clusters is approximately 3. This matches our expectation.

Task

• Perform silhouette analysis on the clusters obtained with \(k\)-means for \(k\) = 2 to 5.

• What’s the suggested number of clusters?

Solution

library(cluster) # for silhouette function

# Initialise some variables
kRange <- seq(from=2, to=5, by=1)
avgWidth <- rep(NA, length(kRange))

# Loop across desired range of ks
for (k in kRange) {
    mdl <- kmeans(x=xTrain , centers=k)
    silh <- silhouette(mdl$cluster, dist(xTrain))
    plot(silh, main=paste0('k=', k))
    avgWidth[k-1] <- mean(silh[, 3])
}

# Plot average width as a function of k 
plot(kRange, avgWidth, type="o", pch=1, col="blue", lwd=3, 
     lty=1, xlab="k", ylab="Average silhouette width")

The “correct” number of clusters is approximately 3. This matches our expectation. However, note that the average silhouette width is quite low, suggesting that some wines are ambiguously assigned to a cluster, and that in general the clusters are not very homogenous.

# Unfortunately in scikit-learn there is not a readily available function
# to plot silhouettes, so we will have to write one ourselves
# Adapted from scikit-learn:
from sklearn.metrics import silhouette_samples, silhouette_score
import numpy as np

#===================================================================#
def plot_silhouette(mdl, xTrain):
    """
    Plots a silhouette plot for a k-means object mdl
    It also requires the training data xTrain
    """
    # Compute the silhouette scores for all samples
    clustLabels = mdl.labels_
    silhScores = silhouette_samples(xTrain, clustLabels)
    
    # Loop across all clusters
    hFig = plt.figure(figsize=(6,6))
    yMin = 10
    for i in range(mdl.n_clusters):
        # Aggregate scores for that cluster and sort them
        thisSilhScores = silhScores[clustLabels == i]
        thisSilhScores.sort()
        
        # Set plot limits and plot
        yMax = yMin + len(thisSilhScores)
        plt.fill_betweenx(np.arange(yMin, yMax), 0, thisSilhScores, 
                          facecolor='grey', edgecolor='black', alpha=0.7)
        plt.xlabel('Silhouette width $s_i$')
        plt.title('k={}, n={}\n Average silhouette width: {:.2f}'\
                  .format(mdl.n_clusters, xTrain.shape[0], np.mean(silhScores)))

        # Leave space before plotting next cluster
        yMin = yMax + 10
    
    # Show plot
    plt.show()
    
    return hFig
    
#===================================================================#
    
# Initialise some variables
kRange = range(2, 6)
avgWidth = np.empty(len(kRange))

# Loop across desired range of ks
for k in kRange:
    mdl = KMeans(n_clusters=k).fit(X=xTrain)
    avgWidth[k-2] = silhouette_score(xTrain, mdl.labels_)
    plot_silhouette(mdl, xTrain)

# Plot average width as a function of k

plt.figure(figsize=(5,5))
plt.plot(kRange, avgWidth, 'bo-')
plt.xlabel('k')
plt.ylabel('Average silhouette width')

The “correct” number of clusters is approximately 3. This matches our expectation. However, note that the average silhouette width is quite low, suggesting that some wines are ambiguously assigned to a cluster, and that in general the clusters are not very homogenous.

---

1. Nowadays microarrays have been largely replaced by sequencing technologies. However, the problem remains exactly the same↩