Activity 4.2 - Kmeans, PAM, and DBSCAN clustering

SUBMISSION INSTRUCTIONS

Render to html
Publish your html to RPubs
Submit a link to your published solutions

Loading required packages:

library(cluster)
library(dbscan)
library(factoextra)
library(tidyverse)
library(patchwork)
library(ggrepel)

Question 1

Reconsider the three data sets below. We will now compare kmeans, PAM, and DBSCAN to cluster these data sets.

three_spheres <- read.csv('Data/cluster_data1.csv')
ring_moon_sphere <- read.csv('Data/cluster_data2.csv')
two_spirals_sphere <- read.csv('Data/cluster_data3.csv')

A)

With kmeans and PAM, we can specify that we want 3 clusters. But recall with DBSCAN we select minPts and eps, and the number of clusters is determined accordingly. Use k-nearest-neighbor distance plots to determine candidate epsilon values for each data set if minPts = 4. Add horizontal line(s) to each plot indicating your selected value(s) of \(\epsilon.\)

library(dbscan)
kNNdistplot(multishapes[,1:2], minPts = 5)
abline(h = 0.35)

library(dbscan)
kNNdistplot(multishapes[,1:2], minPts = 5)
abline(h = 0.17)

library(dbscan)
kNNdistplot(multishapes[,1:2], minPts = 5)
abline(h = 0.30)

B)

Write a function called plot_dbscan_results(df, eps, minPts). This function takes a data frame, epsilon value, and minPts as arguments and does the following:

Runs DBSCAN on the inputted data frame df, given the eps and minPts values;
Creates a scatterplot of the data frame with points color-coded by assigned cluster membership. Make sure the title of the plot includes the value of eps and minPts used to create the clusters!!

Using this function, and your candidate eps values from A) as a starting point, implement DBSCAN to correctly identify the 3 cluster shapes in each of the three data sets. You will likely need to revise the eps values until you settle on a “correct” solution.

library(dbscan)

plot_dbscan_results <- function(df, eps, minPts) {

  X <- scale(df)

  db <- dbscan(X, eps = eps, minPts = minPts)

  plot(X,
       col = db$cluster + 1, pch = 19,
       main = paste("DBSCAN (eps =", eps, ", minPts =", minPts, ")"))

  db
}

minPts <- 4

plot_dbscan_results(three_spheres, eps = 0.3, minPts = minPts)

DBSCAN clustering for 300 objects.
Parameters: eps = 0.3, minPts = 4
Using euclidean distances and borderpoints = TRUE
The clustering contains 3 cluster(s) and 3 noise points.

  0   1   2   3 
  3  99 100  98 

Available fields: cluster, eps, minPts, metric, borderPoints

plot_dbscan_results(ring_moon_sphere, eps = 0.265, minPts = minPts)

DBSCAN clustering for 750 objects.
Parameters: eps = 0.265, minPts = 4
Using euclidean distances and borderpoints = TRUE
The clustering contains 3 cluster(s) and 1 noise points.

  0   1   2   3 
  1 250 249 250 

Available fields: cluster, eps, minPts, metric, borderPoints

plot_dbscan_results(two_spirals_sphere, eps = 0.30, minPts = minPts)

DBSCAN clustering for 815 objects.
Parameters: eps = 0.3, minPts = 4
Using euclidean distances and borderpoints = TRUE
The clustering contains 3 cluster(s) and 0 noise points.

  1   2   3 
257 258 300 

Available fields: cluster, eps, minPts, metric, borderPoints

C)

Compare your DBSCAN solutions to the 3-cluster solutions from k-means and PAM. Use the patchwork package and your function from B) to produce a 3x3 grid of plots: one plot per method/data set combo. Comment on your findings.

library(cluster)
library(dbscan)
library(factoextra)
library(tidyverse)
library(patchwork)
library(ggrepel)

three_spheres      <- read.csv("Data/cluster_data1.csv")
ring_moon_sphere   <- read.csv("Data/cluster_data2.csv")
two_spirals_sphere <- read.csv("Data/cluster_data3.csv")

plot_kmeans_results <- function(df, k = 3) {
  X <- scale(df %>% dplyr::select(where(is.numeric)))
  X <- as.data.frame(X)
  names(X)[1:2] <- c("X1","X2")

  km <- kmeans(X, centers = k, nstart = 20)

  ggplot(X, aes(x = X1, y = X2, col = factor(km$cluster))) +
    geom_point() +
    labs(title = "k-means", col = "cluster") +
    theme_minimal() +
    theme(legend.position = "none")
}

plot_pam_results <- function(df, k = 3) {
  X <- scale(df %>% dplyr::select(where(is.numeric)))
  X <- as.data.frame(X)
  names(X)[1:2] <- c("X1","X2")

  pm <- pam(X, k = k)

  ggplot(X, aes(x = X1, y = X2, col = factor(pm$clustering))) +
    geom_point() +
    labs(title = "PAM", col = "cluster") +
    theme_minimal() +
    theme(legend.position = "none")
}

plot_dbscan_results <- function(df, eps, minPts) {
  X <- scale(df %>% dplyr::select(where(is.numeric)))
  X <- as.data.frame(X)
  names(X)[1:2] <- c("X1","X2")

  db <- dbscan(X, eps = eps, minPts = minPts)

  X$cluster <- db$cluster
  X <- X %>% filter(cluster != 0) 

  ggplot(X, aes(x = X1, y = X2, col = factor(cluster))) +
    geom_point() +
    labs(title = paste("DBSCAN (eps =", eps, ")"), col = "cluster") +
    theme_minimal() +
    theme(legend.position = "none")
}

minPts <- 4

eps1 <- 0.3
eps2 <- 0.265
eps3 <- 0.30

p1 <- plot_kmeans_results(three_spheres, k = 3)
p2 <- plot_pam_results(three_spheres, k = 3)
p3 <- plot_dbscan_results(three_spheres, eps = eps1, minPts = minPts)

p4 <- plot_kmeans_results(ring_moon_sphere, k = 3)
p5 <- plot_pam_results(ring_moon_sphere, k = 3)
p6 <- plot_dbscan_results(ring_moon_sphere, eps = eps2, minPts = minPts)

p7 <- plot_kmeans_results(two_spirals_sphere, k = 3)
p8 <- plot_pam_results(two_spirals_sphere, k = 3)
p9 <- plot_dbscan_results(two_spirals_sphere, eps = eps3, minPts = minPts)


(p1 | p2 | p3) /
(p4 | p5 | p6) /
(p7 | p8 | p9)

All three methods perform well on the first dataset because the clusters are compact and convex. However, for the ring and spiral datasets, k-means and PAM struggle since they tend to split non-convex shapes into pieces. With this eps, DBSCAN does a much better job separating these non-convex clusters.

Question 2

In this question we will apply cluster analysis to analyze economic development indicators (WDIs) from the World Bank. The data are all 2020 indicators and include:

life_expectancy: average life expectancy at birth
gdp: GDP per capita, in 2015 USD
co2: CO2 emissions, in metric tons per capita
fert_rate: annual births per 1000 women
health: percentage of GDP spent on health care
imports and exports: imports and exports as a percentage of GDP
internet and electricity: percentage of population with access to internet and electricity, respectively
infant_mort: infant mortality rate, infant deaths per 1000 live births
inflation: consumer price inflation, as annual percentage
income: annual per-capita income, in 2020 USD

wdi <- read.csv('Data/wdi_extract_clean.csv') 
head(wdi)

      country life_expectancy        gdp      co2 fert_rate    health internet
1 Afghanistan        61.45400   527.8346 0.180555     5.145 15.533614  17.0485
2     Albania        77.82400  4437.6535 1.607133     1.371  7.503894  72.2377
3     Algeria        73.25700  4363.6853 3.902928     2.940  5.638317  63.4727
4      Angola        63.11600  2433.3764 0.619139     5.371  3.274885  36.6347
5   Argentina        75.87800 11393.0506 3.764393     1.601 10.450306  85.5144
6     Armenia        73.37561  4032.0904 2.334560     1.700 12.240562  76.5077
  infant_mort electricity  imports inflation  exports    income
1        55.3        97.7 36.28908  5.601888 10.42082  475.7181
2         8.1       100.0 36.97995  1.620887 22.54076 4322.5497
3        20.4        99.7 24.85456  2.415131 15.53520 2689.8725
4        42.3        47.0 27.62749 22.271539 38.31454 1100.2175
5         8.7       100.0 13.59828 42.015095 16.60541 7241.0303
6        10.2       100.0 39.72382  1.211436 29.76499 3617.0320

Focus on using kmeans for this problem.

A)

My claim: 3-5 clusters appear optimal for this data set. Support or refute my claim using appropriate visualizations.

library(tidyverse)
library(factoextra)
library(cluster)


wdi_num <- wdi%>% select(-country)
wdi_scaled <- scale(wdi_num)
set.seed(415)

fviz_nbclust(wdi_scaled, kmeans, method = "wss", k.max = 10, nstart = 50) +
  labs(title = "Elbow Plot")

fviz_nbclust(wdi_scaled, kmeans, method = "silhouette", k.max = 10, nstart = 50) +
    labs(title = "Silhouette scores")

In the elbow plot, the curve drops fast up to k = 3-4, then it starts flattening out so adding more cluster doesn’t buy you much. In the silhouette plot, the best score is around k = 4, and the 3-5 are all pretty similar. So, after 5, the score gets worse. So, your claim definetly holds up.

B)

Use k-means to identify 4 clusters. Characterize the 4 clusters using a dimension reduction technique. Provide examples of countries that are representative of each cluster. Be thorough.

library(tidyverse)
library(factoextra)

wdi <- read.csv("Data/wdi_extract_clean.csv") %>% drop_na()
Xsc <- scale(wdi %>% select(-country))

set.seed(415)
km4 <- kmeans(Xsc, centers = 4, nstart = 50)

pca <- prcomp(Xsc)
fviz_pca_biplot(pca,
                geom.ind = "point",
                label = "var",
                habillage = factor(km4$cluster),
                addEllipses = TRUE,
                repel = TRUE,
                alpha.ind = 0.35)

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
ℹ The deprecated feature was likely used in the ggpubr package.
  Please report the issue at <https://github.com/kassambara/ggpubr/issues>.

Too few points to calculate an ellipse

d <- sapply(1:nrow(Xsc), \(i){
  cl <- km4$cluster[i]
  sqrt(sum((Xsc[i,] - km4$centers[cl,])^2))
})

wdi %>%
  mutate(cluster = factor(km4$cluster), dist_to_center = d) %>%
  group_by(cluster) %>%
  arrange(dist_to_center) %>%
  slice_head(n = 5) %>%
  select(cluster, country)

# A tibble: 16 × 2
# Groups:   cluster [4]
   cluster country    
   <fct>   <chr>      
 1 1       Zimbabwe   
 2 2       Austria    
 3 2       Finland    
 4 2       Netherlands
 5 2       Belgium    
 6 2       Korea, Rep.
 7 3       Mexico     
 8 3       Tunisia    
 9 3       Mauritius  
10 3       Paraguay   
11 3       Romania    
12 4       Togo       
13 4       Zambia     
14 4       Gambia, The
15 4       Congo, Rep.
16 4       Madagascar

as.data.frame(Xsc) %>%
  mutate(cluster = factor(km4$cluster)) %>%
  group_by(cluster) %>%
  summarise(across(everything(), mean), .groups = "drop")

# A tibble: 4 × 13
  cluster life_expectancy    gdp    co2 fert_rate  health internet infant_mort
  <fct>             <dbl>  <dbl>  <dbl>     <dbl>   <dbl>    <dbl>       <dbl>
1 1               -1.53   -0.668 -0.703     0.976 -1.41     -1.37        1.50 
2 2                1.17    1.45   1.13     -0.715  0.689     1.02       -0.845
3 3                0.0605 -0.378 -0.200    -0.355 -0.0691    0.171      -0.280
4 4               -1.34   -0.660 -0.719     1.54  -0.531    -1.43        1.49 
# ℹ 5 more variables: electricity <dbl>, imports <dbl>, inflation <dbl>,
#   exports <dbl>, income <dbl>

C)

Remove Ireland, Singapore, and Luxembourg from the data set. Use k-means to find 4 clusters again, with these three countries removed. How do the cluster definitions change?

library(tidyverse)
library(factoextra)
library(ggrepel)

wdi_c <- wdi %>%
  filter(!country %in% c("Ireland", "Singapore", "Luxembourg")) %>%
  drop_na()

Xsc_c <- scale(wdi_c %>% select(-country))

set.seed(415)
km4_c <- kmeans(Xsc_c, centers = 4, nstart = 50)

pca_c <- prcomp(Xsc_c)

pcs_c <- as.data.frame(pca_c$x[, 1:2])
pcs_c$Country <- wdi_c$country
pcs_c$Cluster <- factor(km4_c$cluster)
pcs_c$outlier <- abs(scale(pcs_c$PC1)) > 1.4 | abs(scale(pcs_c$PC2)) > 1.2

fviz_pca_biplot(pca_c,
                geom.ind = "point",
                label = "var",
                habillage = factor(km4_c$cluster),
                addEllipses = TRUE,
                repel = TRUE,
                alpha.ind = 0.35)

After removing these countries, cluster 4 is now clear low development group as it sits on the left with high fertility and infant mortality. The cluster 3 separates out a trade-heavy/open economy group as it drops downward in the direction of imports and exports. The cluster 2 is the clear high-development group as it’s on the right, in the direction of high income/GDP and high access. The cluster 1 becomes the middle/mixed group near the center i.e. countries that aren’t extreme on development or trade fall on this one.