STM1001: Computer Lab 7B

Data Science Stream

Topic 7B: Big Data I (Clustering)

Welcome to the seventh computer lab for the Data Science stream of STM1001.

In this computer lab we will carry out a variety of clustering analyses using penguin data from the palmerpenguins R package (Horst, Hill, and Gorman 2020).

By the end of this lab, you should feel comfortable applying various clustering techniques in R. Let’s get started!

1 Cluster Analysis

A cluster is a subset of a data set that consists of observations which, using a chosen metric, are similar to each other, and also dissimilar to other observations¹.

Cluster Analysis is the method of grouping data into clusters, using a clustering technique. There is a large variety of clustering techniques, and each of these techniques uses a different metric for determining the similarity between observations. Don’t worry, we won’t go into all the complicated mathematics involved in these techniques - our focus is on learning how to conduct cluster analyses in R.

Purpose

Clustering is often used as part of an exploratory data analysis procedure, to uncover hidden patterns in a new data set. The main purposes of clustering² are to:

1. Analyse the data structure
2. Relate the different elements of the data to each other, and then, most importantly
3. Aid in classifying the data into certain classes

It is worth noting that in general, we may not be aware of what these classes are, prior to the clustering analysis³.

Clustering Techniques

In this computer lab, we will apply the following popular clustering techniques:

• \(k\)-means Clustering
• PAM Clustering (optional)
• Hierarchical Clustering (optional)
• Fuzzy Clustering (optional)

We will provide brief explanations of each of these techniques as we work through this computer lab.

2 Preparations

In R, we can use the inbuilt cluster package to carry out a range of cluster analyses. However, we will also require a few additional packages⁴. Please run the R code below to install and load these packages.

# Specify required packages
cluster_packages <- c("cluster", "factoextra", "ggfortify", "palmerpenguins")
# Install missing packages
install.packages(setdiff(cluster_packages, rownames(installed.packages())))
# Load all packages
lapply(cluster_packages, library, character.only = TRUE)

Note: If you are in the Data Science stream, you should already have the palmerpenguins package installed.

2.1 Penguin Data

The penguins data set in the palmerpenguins R package (Horst, Hill, and Gorman 2020) contains recent data on 3 species of penguin (Adelie, Chinstrap and Gentoo) living on islands in the Palmer archipelago, off the coast of Antarctica. We will use this data for our clustering analyses.

This penguin data contains information on numerous variables. The table below provides an example of some recorded values for these different variables (the column names), for 3 different penguins:

species	island	bill_length_mm	bill_depth_mm	flipper_length_mm	body_mass_g	sex	year
Adelie	Torgersen	39.1	18.7	181	3750	male	2007
Gentoo	Biscoe	49.2	15.2	221	6300	male	2007
Chinstrap	Dream	45.5	17.0	196	3500	female	2008

2.2

Run the R code below to produce a matrix of scatter plots for the different variables in the penguins data set:

plot(penguins)

Note here that each scatter plot shows two variables from the penguins data set, plotted against each other.

The variable names are shown on the diagonal boxes - just cross reference a plot to the variable names in that row and column to determine the variables plotted on the \(y\)-axis and \(x\)-axis respectively.

Hint: For example, the plot shown in the sixth row, fifth column should be body_mass_g plotted on the y-axis against flipper_length_mm on the x-axis.

2.3

Looking at the scatter plots for the continuous random variables, does it appear that there are any natural clusters that will be easy to identify using a clustering technique?

2.4 Preparing data for Cluster Analysis

For the purposes of this computer lab, we will assume that we don’t actually have data on the species of each penguin in the penguins data set. Instead, we will conduct cluster analyses using the other information contained in this data set.

If the cluster technique works well, we may end up with clusters for each penguin species!

However, before we begin our clustering analyses, we will need to ensure that our data is in the appropriate format. It is important to note that clustering algorithms only work on continuous variables. As we have seen, the penguins data set contains both continuous and categorical (discrete) variables.

Therefore, for the time being, we will only use a subset of the variables in the penguins data set, namely bill_length_mm, bill_depth_mm, flipper_length_mm and body_mass_g.

2.4.1

Run the R code below to remove missing values from our data, and create a subset of the penguins data that only contains the continuous variables:

penguins <- na.omit(penguins)
penguins_subset <- penguins[, 3:6]

2.4.2

Before conducting a cluster analysis, it is also generally a good idea to normalise our data (this process is like the standardisation process covered in section 4.2 of Topic 3). This will remove the effect of different variables being measured on different scales (e.g. flipper length in mm and body mass in grams), and can prevent any single variable overpowering others when being assessed by the cluster algorithm.

Normalising the data consists of two steps:

1. First, we compute the sample mean and sample standard deviation for each variable.
2. Second, we subtract the relevant sample mean from each observation, and then divide by the relevant sample standard deviation.

For example, across the whole data set the sample mean for bill_depth_mm is roughly 17.151, and the sample standard deviation is roughly 1.975. The first penguin in the data set has a recorded bill_depth_mm value of 18.7. The normalised version of this value would therefore be \[\displaystyle \frac{18.7 - 17.151}{1.975} \approx 0.78.\]

Fortunately, we can conduct this normalisation process with one line of code in R, for all our variables, as outlined below.

2.4.3

Run the R code below to assign the normalised penguin data to the object penguins_scaled. Before moving on, inspect your normalised data using the head function.

penguins_scaled <- scale(penguins_subset)

3 \(k\)-means Clustering

Now that our data is suitably prepared, we can begin our cluster analyses.

The first clustering technique we will consider is \(k\)-means clustering, which is a form of centroid-based clustering. A simplified description of \(k\)-means clustering is presented below.

• We begin by arbitrarily choosing a number of clusters \(k\) to use.

• Each data point from our data set is assigned to one of these \(k\) initial clusters, based on the distance between the point and the mean of all points.

• The centroid, i.e. the mean of each cluster, is then calculated, and an iterative process begins:

  • Points are moved between clusters, one point at a time, depending on how close they are to each centroid.
  • This process continues, until no point can be moved between clusters without increasing the average distance between the points and the centroids⁵. At this point the clustering stops.

3.1

We can use the kmeans function from the inbuilt cluster R package to conduct a \(k\)-means clustering analysis. Run the following R code to conduct this analysis on our penguins_scaled data, with 3 clusters specified.

kmeans_fit3 <- kmeans(penguins_scaled, 3)

3.2

Let’s discuss the results. If we run the line kmeans_fit3, a list of output appears in R.

We are primarily interested in the cluster membership assigned by the algorithm. We can extract this cluster membership using kmeans_fit3$cluster.

It is also important to check the details in the Within cluster sum of squares by cluster section. Here:

• between_SS denotes the sum of squares between clusters (i.e., how well the clusters are separated from each other),
• total_ss denotes the total variability in the data.
• The calculation between_ss /total_ss tells us how much of the variability in the data is accounted for by the clusters found by the algorithm. This value can range from 0% to 100%, with larger values indicating a better result. Here, 72.1% is decent, although we might have expected a slightly higher result.

3.3 Visualising our results

We can visualise our results in several ways.

Firstly, run the following code to produce a matrix of scatter plots of the continuous variables, coloured by cluster.

windows() # or quartz() if you are on a Mac
plot(penguins[, c(1,3:6)], col = kmeans_fit3$cluster)

Don’t worry if the plots produced look a bit confusing. To make the plots easier to read, consider maximising the graphics window.

Recall from 2.2 that each scatter plot here shows two of the variables in the penguins data set plotted against each other. Now however, they are coloured by cluster (as determined by the \(k\)-means method).

If, by inspecting these plots, you can clearly distinguish between clusters, then the clustering method has performed well.

Hint: If you are having trouble interpreting the graph, this note may help. For example, the plot shown in the fifth row, fourth column should be body_mass_g plotted on the y-axis against flipper_length_mm on the x-axis.

3.3.1

We could also use the R code below to plot the species variable as a numeric value (Adelie = 1, etc), with the points coloured by cluster:

plot(as.numeric(penguins$species), col = kmeans_fit3$cluster)

Based on this plot, and the plot produced in 3.3, do you think the \(k\)-means clustering has performed well?

3.4

At this point, it might seem that our job is done. We have results for \(k=3\), and there are three species of penguins.

Remember though, normally we would not actually know what the number of clusters should be. Therefore, we would usually try a range of different \(k\) values.

Conduct \(k\)-means clustering on our penguins_scaled data, for \(k\) values of 2 and 4.

Note down the between_ss /total_ss results. What do you observe?

3.5 Diagnostics

You might have noticed that the between_ss /total_ss value increases as we increase the number of clusters. This can happen even if adding another cluster isn’t actually that helpful. Therefore, we should also consider some other assessment methods.

If we did not know the number of clusters, we could use several methods to determine the appropriate value of \(k\).

The factoextra R package contains a helpful function fviz_nbclust, which we can use to produce several informative plots.

Take a look at the comments in the Code chunk below (the ... signify that these are comments, and should be replaced by actual arguments):

example <- fviz_nbclust(...specify data set here... , 
                        ...specify clustering method here... , 
                        method = "...specify assessment method here...")

As you can see, the fviz_nbclust function takes three main arguments:

1. The data set,
2. The clustering method (i.e. kmeans), and
3. The assessment method (one of “wss” , “silhouette” , or “gap_stat”)

For the purposes of this lab, we will focus on just the silhouette (see Rousseeuw 1987) and gap methods.

3.5.1 The silhouette method

The silhouette assessment method measures the similarity of a point to other points in its cluster.

• The silhouette value for each observation can be between 0 and 1.
• Values close to 0 suggest clusters are very similar, while values close to 1 indicate clusters are very different. Thus, the higher the average silhouette width value (across all points in a cluster), the better.

Using the fviz_nbclust function, with the kmeans clustering method and the "silhouette" assessment method specified, determine the optimal number of clusters for our penguins_scaled data.

3.5.2 The gap statistic method

The second \(k\)-means clustering assessment method we will consider is the gap_stat assessment method. This computes a statistic known as the ‘gap statistic’. We won’t go into the mathematical details of this statistic; suffice to say that higher values are considered preferable.

Using the fviz_nbclust function, with the kmeans clustering method and the "gap_stat" assessment method specified, determine the optimal number of clusters for our penguins_scaled data.

3.5.3 Cluster plots

We can also use the fviz_cluster function from the factoextra R package to visualise the clusters.

Run the code below to visualise the three clusters found in 3.1⁶, and then, using this code as a base, visualise the sets of clusters found in 3.4.

fviz_cluster(kmeans_fit3, data = penguins_scaled)

3.6

In conclusion, based on all your \(k\)-means clustering analyses, which value of \(k\) would you recommend using, and why?

At this point, we’ve covered the main material for this computer lab - well done! If you would like, you can extend your knowledge on clustering by going through the Extension sections below. Otherwise, if you have time left, you might like to work on your assessments.

4 Extension 1: PAM Clustering

An alternative to \(k\)-means clustering is PAM (Partition around Medoids) clustering.

PAM clustering operates in a similar manner to \(k\)-means clustering, but uses medoids (median-like points) rather than centroids, when deciding cluster membership. This makes the PAM clustering technique more robust to outliers.

4.1

We can easily conduct PAM clustering using the pam function from the inbuilt cluster R package. The process is similar to the \(k\)-means clustering process conducted in 3.1.

Use the pam function to conduct PAM clustering for \(k\) values of 2, 3 and 4.

4.2

Next, use the fviz_nbclust function to determine the best number of PAM clusters to use, via the silhouette method.

Hint: Remember to use pam instead of kmeans for the clustering method specification.

4.3

Finally, use the fviz_cluster function and the code below to visualise the clusters for the PAM clustering result(s) you chose above in 4.2. What do you conclude?

# This code assumes you are assessing a result stored in the object pam_fit3
windows()
plot(as.numeric(penguins$species), col = pam_fit3$clustering)

5 Extension 2: Hierarchical Clustering

Unlike the previous clustering techniques discussed, hierarchical clustering begins by assigning each point to its own cluster.

This means that, for our penguins_scaled data set, we would start with 333 clusters - 1 for each penguin. Then, the two most similar clusters are merged (continuing our example, this would result in us now having 332 clusters). We repeat this process, until we have just one large cluster (which contains all the points).

As a result, we now have a hierarchy of clusters, which looks a bit like an upside-down tree (with observations from the same cluster on the same ‘branch’).

5.1

We can conduct hierarchical clustering using the agnes function from the cluster R package. Since we are not specifying the number of clusters, we only need to provide the one argument, i.e. the data set to analyse.

Use the agnes function to conduct hierarchical clustering on our penguins_scaled data set.

5.2

We can visualise the results of the hierarchical clustering using a dendrogram.

Complete the code below to create a dendrogram for your results (just replace the ... with your results)

plot(..., which = 2)

5.3

There are several linkage method options we can choose from when conducting hierarchical clustering. These relate to the way in which distance is measured between two points, to determine their level of similarity.

The default method in the agnes function is average. Try re-running your hierarchical clustering algorithm using some other options, by adding the argument method = "..." to your code in 5.1, and replacing the "..." with "single" and then "ward".

When you visualise the results, do the dendrograms look very different?

5.4

So far, the dendrograms we have produced looked rather unappealing - we can do better.

It is worth noting that the selection of an appropriate horizontal cut-off point on the dendrogram (thus selecting the number of clusters to use) is arbitrary. Based on your findings from the previous analyses however, you should now have decided on an ‘optimal’ number of \(k\) clusters to select.

We can visualise the dendrogram with the tree cut into these \(k\) clusters using the hcut and fviz_dend functions from the factoextra R package.

Take a look at the code below, and fill in the ... missing parts with the relevant details. Once you are happy with the code, run it. The resultant dendrogram should look much nicer.

penguin_dendro <- hcut(..., 
                       k = ..., 
                       hc_func = "agnes",
                       hc_method = "...",
                       hc_metric = "euclidean") 

fviz_dend(penguin_dendro)

Note: You will need to specify the measurement method, e.g. "ward" in the hc_method = argument.

6 Extension 3: Fuzzy Clustering

For the previous clustering methods, each data point belonged to only one cluster. This is sometimes referred to as hard clustering.

In contrast, in fuzzy clustering each point is allowed to belong to multiple clusters. The more similar the point is to other points in a cluster, the higher that point’s percentage of membership to that cluster⁷. This could be helpful when a point naturally should belong to multiple clusters (for instance, if penguins were clustered into 6 clusters, for males and females of each species, then it would be reasonable for each penguin to belong to two clusters).

6.1

We can easily conduct fuzzy clustering using the fanny function from the inbuilt cluster R package. The process is similar to the \(k\)-means clustering process (conducted in 3.1) and the PAM clustering process (conducted in 4.1).

Use the fanny function to conduct fuzzy clustering for \(k\) values of 2, 3 and 4.

Note: The output for these analyses will include a membership coefficients list of the membership percentages for each point. While it is not easy to visualise these, the output also includes the closest ‘hard’ cluster for each point, which can be called via $cluster, as per the previous clustering techniques.

6.2

Next, use the fviz_nbclust function to determine the best number of fuzzy clusters to use, via the silhouette method.

Note: Don’t worry if any warning messages like the one shown in the code chunk below appear - you should still be able to produce the plots:

## Warning in FUNcluster(x, i, ...): FANNY algorithm has not converged in 'maxit' =
## 500 iterations

6.3

Use the fviz_cluster function and the code below to visualise the clusters for the fuzzy clustering result you chose above in 6.2. What do you conclude?

# This code assumes you are assessing a result stored in the object fuzzy_fit3
windows()
plot(as.numeric(penguins$species), col = fuzzy_fit3$cluster)

Great work, that’s everything for today!

That concludes our work on clustering techniques. Did you have a preference for one of the four techniques introduced in this computer lab?

References

Gan, G., C. Ma, and J. Wu. 2007. Data Clustering: Theory, Algorithms, and Applications. ASA-SIAM Series on Statistics and Applied Probability ; 20. Philadelphia, Pa.: Society for Industrial; Applied Mathematics (SIAM, 3600 Market Street, Floor 6, Philadelphia, PA 19104). https://doi.org/https://doi.org/10.1137/1.9780898718348.

Horikoshi, M., Y. Tang, A. Dickey, M. Grenie, R. Thompson, L. Selzer, D. Strbenac, K. Voronin, and D. Pulatov. 2021. ggfortify: Data Visualization Tools for Statistical Analysis Results. https://github.com/sinhrks/ggfortify.

Horst, Allison Marie, Alison Presmanes Hill, and Kristen B Gorman. 2020. Palmerpenguins: Palmer Archipelago (Antarctica) Penguin Data. https://doi.org/10.5281/zenodo.3960218.

Kassambara, A., and F. Mundt. 2020. factoextra: Extract and Visualize the Results of Multivariate Data Analyses. http://www.sthda.com/english/rpkgs/factoextra.

Mirkin, B. 1996. Mathematical Classification and Clustering. 1st ed. 1996.. Nonconvex Optimization and Its Applications, 11.

R Core Team. 2021. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Rousseeuw, P. J. 1987. “Silhouette: A Graphical Aid to the Interpretation and Validation of Cluster Analysis.” Computational and Applied Mathematics 20: 53–65.

Thulin, M. 2021. Modern Statistics with R: From Wrangling and Exploring Data to Inference and Predictive Modelling.

These notes have been prepared by Rupert Kuveke. Please note that some of the content in these notes has been developed from content in Thulin (2021). The copyright for the material in these notes resides with the authors named above, with the Department of Mathematical and Physical Sciences and with La Trobe University. Copyright in this work is vested in La Trobe University including all La Trobe University branding and naming. Unless otherwise stated, material within this work is licensed under a Creative Commons Attribution-Non Commercial-Non Derivatives License BY-NC-ND.

---

1. Mirkin (1996), p.25↩︎

2. see e.g. Mirkin (1996), p.25↩︎

3. see e.g. Gan, Ma, and Wu (2007)↩︎

4. the factoextra package was created by Kassambara and Mundt (2020) and the ggfortify package was created by Horikoshi et al. (2021). The cluster package is part of the base R suite of packages (R Core Team 2021).↩︎

5. See e.g. section 4.10.3 of Thulin (2021)↩︎

6. Note that these might look a little different to those visualised in 3.3, since the axes variables are different↩︎

7. See e.g. Thulin (2021)↩︎