STM1001: Computer Lab 7B Solutions

Data Science Stream

Topic 7B: Big Data I (Clustering)

Example R code solutions for the Data Science Computer Lab 7, which uses penguin data from the palmerpenguins R package¹ (Horst, Hill, and Gorman 2020), are presented below.

1 Cluster Analysis

Purpose

No answer required.

Clustering Techniques

No answer required.

2 Preparations

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

2.1 Penguin Data

No answer required.

2.2

plot(penguins)

2.3

For several of the continuous random variable scatter plots, there does appear to be visible clustering occurring - e.g. flipper_length_mm vs bill_depth_mm suggests that there are two distinct clusters. Therefore applying a clustering technique does seem to be warranted.

2.4 Preparing data for Cluster Analysis

No answer required.

2.4.1

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

2.4.2

No answer required.

2.4.3

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

3.1

set.seed(1) # included for reproducibility purposes
kmeans_fit3 <- kmeans(penguins_scaled, 3)

3.2

We can check the cluster membership as follows:

kmeans_fit3$cluster

##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [38] 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1
##  [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 3
## [149] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [186] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [223] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [260] 3 3 3 3 3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 1 2 2 2 2 2 2 2 1
## [297] 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2

The between_ss /total_ss value is 72.1%.

3.3 Visualising our results

plot(penguins[, c(1,3:6)], col = kmeans_fit3$cluster)

As we can see, the clusters are actually doing a pretty good job of separating penguins by species!

3.3.1

plot(as.numeric(penguins$species), col = kmeans_fit3$cluster, 
     xlab = "penguin ID", 
     ylab = "(Adelie = 1, Chinstrap = 2, Gentoo = 3", las = 1)

Based on these plots, it seems that the \(k\)-means clustering has done a great job of clustering the penguins into clusters of different species. If you look carefully, there is some error when distinguishing between Chinstrap and Adelie penguins. If we also include more data such as sex or island, we might be able to more clearly distinguish between these two species.

3.3.2 Two-way Table Summary

💻

two_way_table <- table(as.numeric(penguins$species), kmeans_fit3$cluster)
dimnames(two_way_table) <- list(levels(penguins$species),
                                c("Cluster 1", "Cluster 2", "Cluster 3")
                                )
two_way_table

##           Cluster 1 Cluster 2 Cluster 3
## Adelie           22         0       124
## Chinstrap        63         0         5
## Gentoo            0       119         0

The two-way table highlights the accuracy of the clustering - only a few Chinstrap penguins and some Adelie penguins have been misclassified.

3.3.3 Cluster plots

fviz_cluster(kmeans_fit3, data = penguins_scaled)

3.4

kmeans_fit2 <- kmeans(penguins_scaled, 2)
kmeans_fit4 <- kmeans(penguins_scaled, 4)

The between_ss /total_ss results for these new cluster fits are:

• 58.5% for \(k=2\).
• 77.1% for \(k=4\).

The between_ss /total_ss value actually increases as we increase the number of clusters, which isn’t necessarily accurate - we know e.g. that there are not four species of penguin in our data (what happens if we use \(k=10\)?!).

3.5 Diagnostics

No answer required.

3.5.1 The silhouette method

kmeans_sil <- fviz_nbclust(penguins_scaled, kmeans, method = "silhouette")
kmeans_sil

Based on this plot, it seems that the optimal number of clusters is two.

3.5.2 The gap statistic method

kmeans_gap <- fviz_nbclust(penguins_scaled, kmeans, method = "gap")
kmeans_gap

Based on this plot, it seems that the optimal number of clusters is four.

3.5.3

fviz_cluster(kmeans_fit2, data = penguins_scaled)

fviz_cluster(kmeans_fit4, data = penguins_scaled)

3.6

Answers may vary here. Three or four clusters seem reasonable.

4 Extension 1: PAM Clustering

4.1

pam_fit2 <- pam(penguins_scaled, 2)
pam_fit3 <- pam(penguins_scaled, 3)
pam_fit4 <- pam(penguins_scaled, 4)

4.2

pam_sil <- fviz_nbclust(penguins_scaled, pam, method = "silhouette")
pam_sil

These results suggest that two clusters are optimal.

4.3

The R code below shows results for a cluster choice of three.

fviz_cluster(pam_fit3, data = penguins_scaled)

plot(as.numeric(penguins$species), col = pam_fit3$clustering)

We can see from these plots that the PAM clustering does a very good job of clustering penguins into species (although there is still some trouble distinguishing Adelie and Chinstrap penguins). For the larger numbers of clusters, the algorithm may also be distinguishing between male and female penguins - try checking if this is true.

5 Extension 2: Hierarchical Clustering

5.1

hier_fit <- agnes(penguins_scaled)

5.2

plot(hier_fit, which = 2)

5.3

hier_fit <- agnes(penguins_scaled)
hier_fit_single <- agnes(penguins_scaled, method = "single")
hier_fit_ward <- agnes(penguins_scaled, method = "ward")

plot(hier_fit_single, which = 2)

plot(hier_fit_ward, which = 2)

The dendrograms all look quite different. The ward method one looks perhaps the cleanest and easiest to assess.

5.4

The code shown below produces a coloured dendrogram with three clusters, for the hierarchical clustering algorithm using the “ward” measurement method.

penguin_dendro <- hcut(penguins_scaled, 
                k = 3, 
                hc_func = "agnes",
                hc_method = "ward",
                hc_metric = "euclidean") 

fviz_dend(penguin_dendro)

6 Extension 3: Fuzzy Clustering

6.1

The R code below shows fits for the different cluster number options.

fuzzy_fit2 <- fanny(penguins_scaled, 2)
fuzzy_fit3 <- fanny(penguins_scaled, 3)
fuzzy_fit4 <- fanny(penguins_scaled, 4)

Please note that output for these fits is omitted due to the length of the results.

6.2

fuzzy_sil <- fviz_nbclust(penguins_scaled, fanny, method = "silhouette")
fuzzy_sil

We obtain the same results as for the \(k\)-means and PAM clustering here, i.e. we have a result of two clusters.

6.3

The R code below shows results for a cluster choice of two.

fviz_cluster(fuzzy_fit2, data = penguins_scaled)

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

From these plots we can see that the fuzzy clustering performs similarly to the previous methods.

That’s all the content covered.

References

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.

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 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. We also utilise the packages factoextra (Kassambara and Mundt 2020), ggfortify (Horikoshi et al. 2021) and cluster(R Core Team 2021).↩︎