Clustering: hierarchical, k-means, and model-based (mclust)

Introduction

In this practical, we apply three clustering approaches. In Part 1, we use hierarchical and k-means clustering on synthetic bivariate data, to build intuition about how distance metrics, linkage, and the number of clusters affect the result, and about the instability of k-means across runs. In Part 2, we move to a real dataset of Swiss banknote measurements and use model-based clustering (mclust) to fit and compare probabilistic cluster models, including choosing the number of clusters and comparing model fit via BIC.

library(MASS)
library(tidyverse)
library(patchwork)
library(ggdendro)
library(mclust)

Make sure to load MASS before tidyverse, otherwise the function MASS::select() will overwrite dplyr::select().

Part 1: Hierarchical and k-means clustering

Before we start, set a seed for reproducibility, we use 123, and also use options(scipen = 999) to suppress scientific notations, making it easier to compare and interpret results later in the session.

set.seed(123)
options(scipen = 999)

Data processing

1. The data can be generated by running the code below. Try to understand what is happening as you run each line of the code below.

# randomly generate bivariate normal data
set.seed(123)
sigma      <- matrix(c(1, .5, .5, 1), 2, 2)
sim_matrix <- mvrnorm(n = 100, mu = c(5, 5), Sigma = sigma)
colnames(sim_matrix) <- c("x1", "x2")

# change to a data frame (tibble) and add a cluster label column
sim_df <- 
  sim_matrix |> 
  as_tibble() |>
  mutate(class = sample(c("A", "B", "C"), size = 100, replace = TRUE))

# Move the clusters to generate separation
sim_df_small <- 
  sim_df |>
  mutate(x2 = case_when(class == "A" ~ x2 + .5,
                        class == "B" ~ x2 - .5,
                        class == "C" ~ x2 + .5),
         x1 = case_when(class == "A" ~ x1 - .5,
                        class == "B" ~ x1 - 0,
                        class == "C" ~ x1 + .5))
sim_df_large <- 
  sim_df |>
  mutate(x2 = case_when(class == "A" ~ x2 + 2.5,
                        class == "B" ~ x2 - 2.5,
                        class == "C" ~ x2 + 2.5),
         x1 = case_when(class == "A" ~ x1 - 2.5,
                        class == "B" ~ x1 - 0,
                        class == "C" ~ x1 + 2.5))

2. Prepare two unsupervised datasets by removing the class feature.

3. For each of these datasets, create a scatterplot. Combine the two plots into a single frame (look up the “patchwork” package to see how to do this!) What is the difference between the two datasets?

Hierarchical clustering

4. Run a hierarchical clustering on these datasets and display the result as dendrograms. Use euclidian distances and the complete agglomeration method. Make sure the two plots have the same y-scale. What is the difference between the dendrograms?

Hint: functions you will need are hclust, ggdendrogram, and ylim.

5. For the dataset with small differences, also run a complete agglomeration hierarchical cluster with manhattan distance.

6. Use the cutree() function to obtain the cluster assignments for three clusters and compare the cluster assignments to the 3-cluster euclidian solution. Do this comparison by creating two scatter plots with cluster assignment mapped to the colour aesthetic. Which difference do you see?

K-means clustering

7. Create k-means clustering with 2, 3, 4, and 6 clusters on the large difference data. Again, create coloured scatter plots for these clustering results.

8. Do the same thing again a few times. Do you see the same results every time? where do you see differences?

9. OPTIONAL (if time permits): Find a way online to perform bootstrap stability assessment for the 3 and 6-cluster solutions.

Hierarchical and k-means clustering both partition observations into hard, non-overlapping groups based on distance. In Part 2, we take a different, probabilistic approach: model-based clustering, which fits a mixture of distributions to the data and lets us evaluate fit statistically (via BIC) rather than only visually.

Part 2: Model-based clustering using mclust

We apply model-based clustering on a data set of bank note measurements. The data is built into the mclust package and can be loaded as a tibble:

df <- as_tibble(banknote)

Data exploration

10. Read the help file of the banknote data set to understand what it’s all about.

11. Create a scatter plot of the left (x-axis) and right (y-axis) measurements on the data set. Map the Status column to colour. Jitter the points to avoid overplotting. Are the classes easy to distinguish based on these features?

12. From now on, we will assume that we don’t have the labels. Remove the Status column from the data set.

13. Create density plots for all columns in the data set. Which single feature is likely to be best for clustering?

Univariate model-based clustering

14. Use Mclust to perform model-based clustering with 2 clusters on the feature you chose. Assume equal variances. Name the model object fit_E_2. What are the means and variances of the clusters?

15. Use the formula from the slides and the model’s log-likelihood (fit_E_2$loglik) to compute the BIC for this model. Compare it to the BIC stored in the model object (fit_E_2$bic). Explain how many parameters (m) you used and which parameters these are.

16. Plot the model-implied density using the plot() function. Afterwards, add rug marks of the original data to the plot using the rug() function from the base graphics system.

17. Use Mclust to perform model-based clustering with 2 clusters on this feature again, but now assume unequal variances. Name the model object fit_V_2. What are the means and variances of the clusters? Plot the density again and note the differences.

18. How many parameters does this model have? Name them.

19. According to the deviance, which model fits better?

20. According to the BIC, which model is better?

Multivariate model-based clustering

We will now use all available information in the data set to cluster the observations.

21. Use Mclust with all 6 features to perform clustering. Allow all model types (shapes), and from 1 to 9 potential clusters. What is the optimal model based on the BIC?

22. How many mean parameters does this model have?

23. OPTIONAL (if time permits): Run a 2-component VVV model on this data. Create a matrix of bivariate contour (“density”) plots using the plot() function. Which features provide good component separation? Which do not?

24. OPTIONAL (if time permits): Create a scatter plot just like the first scatter plot in this tutorial, but map the estimated class assignments to the colour aesthetic.

Map the uncertainty (part of the fitted model list) to the size aesthetic, such that larger points indicate more uncertain class assignments. Jitter the points to avoid overplotting. What do you notice about the uncertainty?

(This uses the VVV2 model from the previous optional exercise. Skip both together if you’re short on time.)