STM1001: Computer Lab 6B

Science/Health Science/Data Science Module

Topic 6B: Big Data I (Clustering)

Welcome to the sixth computer lab for the Science, Health Science and Data Science modules.

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 which are also dissimilar to other observations, via the same metric¹.

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 purpose of clustering² is to:

1. Analyse the data structure
2. Relate the different elements of the data to each other, and then
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
• Fuzzy Clustering
• Hierarchical Clustering

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.

install.packages(c("factoextra", "ggfortify", "palmerpenguins"))
library(cluster)
library(factoextra)
library(ggfortify)
library(palmerpenguins)

Note: If you are in the Data Science module, 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 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 variables such as:

• flipper length
• body mass
• bill length
• sex
• species
• island on which the penguin lives

2.2

Run the R code below to take a look at the penguins data set, to familiarise yourself with the different variables.

head(penguins)
plot(penguins)

Note here that plot(penguins) will produce a matrix of scatter plots.

Each plot shows two of the variables in 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.

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.3 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, to try to group the penguins into clusters.

If the cluster technique works well, we could 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.3.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.3.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:

1. Computing the sample mean and sample standard deviation for each variable
2. Subtracting the relevant sample mean from each observation, and dividing by the relevant sample standard deviation

For example, 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.

2.3.3

Use the scale function to normalise our penguins_subset data, and assign the normalised data to the object penguins_scaled. Then, take a look at your normalised data using the head function.

Hint: Your output should look like the output in the Code chunk below:

penguins_scaled <- scale(penguins_subset)
head(penguins_scaled)

##      bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## [1,]     -0.8946955     0.7795590        -1.4246077  -0.5676206
## [2,]     -0.8215515     0.1194043        -1.0678666  -0.5055254
## [3,]     -0.6752636     0.4240910        -0.4257325  -1.1885721
## [4,]     -1.3335592     1.0842457        -0.5684290  -0.9401915
## [5,]     -0.8581235     1.7444004        -0.7824736  -0.6918109
## [6,]     -0.9312674     0.3225288        -1.4246077  -0.7228585

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.

In k-means clustering, we specify an arbitrary number of clusters \(k\), and each data point is assigned to one of these \(k\) clusters, based on the distance between the point and the mean of all points.

Once these initial clusters are created, the mean of each cluster is calculated, and an iterative process begins:

• Points are moved between clusters, one point at a time, depending on how close they are to each cluster mean.
• This process continues, until no point can be moved between clusters without increasing the average distance between the points and the cluster means⁵.

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 k=3 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), while 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.4

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.5

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, 4 and 5.

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

3.6

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”)

3.6.1

The wss assessment method produces a scree plot, with the number of clusters plotted on the horizontal axis, and the Total Within Sum of Squares (wss) variance plotted on the vertical axis. A lower variance is preferable, but after a certain point, adding more clusters will not lower this variability significantly. When assessing this plot, we select, as our optimal number of clusters, the number of clusters for which adding an extra cluster doesn’t greatly reduce this variability.

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

Hint: If you are not sure how to proceed, check the Code chunk below:

kmeans_wss <- fviz_nbclust(penguins_scaled, kmeans, method = "wss")
windows()
kmeans_wss

3.6.2

The silhouette assessment method measures the similarity of a point to other points in its cluster. The higher this ‘average silhouette width’ , 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.6.3

The final 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.7

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.5.

fviz_cluster(kmeans_fit3, data = penguins_scaled)

3.8

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

4 PAM Clustering

k-means clustering is a form of centroid-based clustering, whereby cluster means (i.e. centroids) are used as a metric in the clustering algorithm. An alternative to k-means clustering is PAM (Partion around Medoids) clustering. This 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, 4 and 5.

4.2

Next, use the fviz_nbclust function to determine the best number of PAM clusters to use (use all three methods).

Hint: 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)

At this point, we’ve covered the main material for this computer lab. 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.

5 Extension 1 - Hierarchical Clustering

Unlike the previous clustering techniques discussed, hierarchical clustering begins by assigning each point to its own cluster. For example, since we have 333 penguins in our filtered data set, we would start with 333 clusters. Then, the two most similar clusters are merged (continuing our eample, 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 scaled penguin data.

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 options we can use 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", "complete" 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 2 - 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 the 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 each penguin would 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, 4 and 5.

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. Just use the wss and silhouette methods (the gap_stat method can take a long time to compute for fuzzy clustering).

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/.

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 Mathematics and Statistics 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)↩︎