---
title: "🗓️ Week 08: Introduction and Clustering"
subtitle: "Theme: Unsupervised Learning"
author: "Dr. Ghita Berrada"
format: 
  html:
    from: markdown+emoji
    page-layout: full
    toc: true
    toc-depth: 2
self-contained: true
---

Welcome to Week 08! This week, we continue our exploration of unsupervised learning techniques with clustering...

# :gear: Setup

Before you can load the libraries below, you might have to do the following package installs (feel free to skip if you already have done the installs before)

``` r
install.packages("fpc")
install.packages("factoextra")
install.packages("dbscan")
install.packages("rnaturalearth") #if unable to install, use the instruction install.packages("rnaturalearth",type="binary")
install.packages("rnaturalearthdata")
install.packages("sf")
install.packages("tidytext")
install.packages("ClusterR")
```

Then, let's load the libraries we'll need today.

```{r, echo=FALSE, message=FALSE, warning=FALSE}
# Tidyverse packages
library(tidyr)
library(dplyr)
library(ggplot2)
library(readr)
library(skimr)
library(tidytext)

# Tidymodels packages, spelt out individually for clarity
library(recipes)

#new libraries for clustering
library(fpc) #library for Calinski-Harabasz index
library(factoextra) # library for elbow and average silhouette plots
library(dbscan) # library for DBSCAN
library(ClusterR)


#libraries for nicer plots
library(rnaturalearth)
library(rnaturalearthdata)
library(sf)
 
```

## Part 1: Looking at earthquakes

Just as usual, we'll introduce theoretical concepts through an example dataset/case study.

The original data is a dataset that comes pre-loaded with R: `quakes` and the case study was inspired by [this tutorial](https://www.r-bloggers.com/2021/03/k-means-101-an-introductory-guide-to-k-means-clustering-in-r/).

The dataset describes earthquakes near Fiji since 1964.

```{r}
skim(quakes)
```

This dataset is formed of 1000 rows and 5 columns: latitude (degrees), longitude (degrees), depth of the earthquake (km), its magnitude (Richter scale) and the number of stations reporting the quake.

We'll only focus on depth and magnitude for now. Let's visualise this data.

```{r}

quakes %>% select(depth,mag) %>%
          ggplot(aes(x=depth,y=mag)) +
          geom_point(alpha=0.4, size=2, color='#9e303b') +
  labs(
    title = "Is there a way to segment earthquakes\n based on depth and magnitude?",
    x = "Earthquake depth (km)",
    y = "Earthquake magnitude (Richter scale)",
  ) 
       
```

# Part 2: Meet the k-means algorithm

## How does it work

Let's take a look at this [visualisation of K-Means](https://www.naftaliharris.com/blog/visualizing-k-means-clustering/).

[Artwork by [\@allison_horst](https://twitter.com/allison_horst)]{style="font-size:0.75em;color:gray"}

![](../../../2023/autumn-term/weeks/week08/figures/k_means1.jpg)

![](images/clipboard-3149433162.png)

![](../../../2023/autumn-term/weeks/week08/figures/k_means2.jpg)

![](../../../2023/autumn-term/weeks/week08/figures/k_means3.jpg)

![](../../../2023/autumn-term/weeks/week08/figures/k_means4.jpg)

![](../../../2023/autumn-term/weeks/week08/figures/k_means5.jpg)

![](../../../2023/autumn-term/weeks/week08/figures/k_means6.jpg)

![](../../../2023/autumn-term/weeks/week08/figures/k_means7.jpg)

![](../../../2023/autumn-term/weeks/week08/figures/k_means8.jpg)

![](../../../2023/autumn-term/weeks/week08/figures/k_means9.jpg)

![](../../../2023/autumn-term/weeks/week08/figures/k_means10.jpg)

![](../../../2023/autumn-term/weeks/week08/figures/k_means11.jpg)

![](../../../2023/autumn-term/weeks/week08/figures/k_means12.jpg)

## k-means in more mathematical terms

Κ-means clustering uses **iterative refinement** to produce a final result.

*Algorithm inputs*: number of clusters **k** and data set.

-   The data set is a collection of features for each data point.
-   The algorithms starts with initial estimates for the Κ centroids, which can either be randomly generated or randomly selected from the data set.
-   The algorithm then iterates between two steps:

1.  Data assignment step:

Each centroid defines one of the clusters. In this step, each data point is assigned to its nearest centroid, based on the **squared Euclidean distance**. More formally, if $c_i$ is the collection of centroids in set $C$, then each data point $x$ is assigned to a cluster based on

$\underset{c_i \in C}{\arg\min} \; dist(c_i,x)^2$

where dist( · ) is the standard (**L2**) Euclidean distance. Let the set of data point assignments for each i-th cluster centroid be $S_i$.

2.  Centroid update step:

In this step, the centroids are recomputed. This is done by taking the mean of all data points assigned to that centroid's cluster.

$c_i=\frac{1}{|S_i|}\sum_{x_i \in S_i} x_i$

The algorithm iterates between steps one and two until a stopping criteria is met (i.e., no data points change clusters, the sum of the distances is minimized, or some maximum number of iterations is reached).

This algorithm is guaranteed to converge to a result. The result may be a local optimum (i.e. not necessarily the best possible outcome), meaning that assessing more than one run of the algorithm with randomized starting centroids may give a better outcome.

:warning:

-   K-means centroids generally initialized randomly, which doesn't guarantee an optimal clustering:

    -   A variant initialization technique is called **K-means++**. The basic idea is to select the first centroid randomly then assign the remaining centroids based on the maximum squared distance: the goal is to push the centroids as far as possible from one another.

    Here is the detailed K-means++ initialization algorithm:

    1.  Assign $c_1$ randomly ($c_1$ being the first centroid)
    2.  Calculate the distance between all data points and the selected centroid to get $D_i=\max_{j\in[1,k]} ||x_i-c_j||^2$ (i.e the distance of a data point $x_i$ from the farthest centroid $c_j$).
    3.  Initialize the data point $x_i$ as the next new centroid
    4.  Repeat steps 2 and 3 till all the defined $K$ centroids are found

See (Arthur and Vassilvitskii 2007) for more details.

-   K-means is sensitive to outliers. One way to avoid that is to select medoids (i.e the actual points -that exist in the data- within each cluster, whose dissimilarities with all the other points in the cluster are minimum as opposed to centroids i.e the average of all points in each cluster in k-means). This is the k-medoid algorithm.

## Part 3: Using the k-means algorithm

Let's see what k-means comes up when applied to our depth+magnitude data.

:dart: **ACTION POINT**

(In pairs)

1.  Play with the number of clusters and see what the k-means function returns. Can you make sense of the output?

```{r}
# Let's filter our data frame to focus on depth and magnitude

n_clusters <- 4 # How many do you think we should have?

quakes_df <- 
  quakes %>% 
  select(depth,mag) 

### THIS IS HOW YOU RUN K-MEANS IN BASE R ###
fitted_kmeans <- kmeans(quakes_df, centers = n_clusters, iter.max = 10, nstart = 5)
fitted_kmeans
```

2.  Visualize the clusters

```{r}
# Visualize the clusters
quakes_df <- 
  quakes_df %>%
  mutate(cluster = as.factor(fitted_kmeans$cluster))

g <-
  ggplot(quakes_df, aes(x = depth, y = mag, color = cluster)) +
  geom_point(alpha=0.4, size=2) +
  labs(
    title = paste("Clusters found by k-means with", n_clusters, "clusters"),
    subtitle = "Do you like the clusters?\n Do they make any sense to you?",
    x = "Earthquake depth (km)",
    y = "Earthquake magnitude (Richter scale)",
    color = "Cluster"
  ) +
  theme_bw(base_size=15) 

g
```

It looks as if only the variable depth was used when doing the clustering. Why is that?

Let's change this!

```{r}
quakes_features<- recipe(~ ., data = quakes) %>%
  step_rm(stations) %>% 
  step_normalize(all_predictors()) %>% 
  prep() %>% 
  bake(new_data = NULL)

quakes_features %>% 
  slice_head(n = 5)
```

```{r}
# Let's filter our normalized data frame to focus on depth and magnitude

n_clusters <- 4 # How many do you think we should have?

normalized_quakes_df <- 
  quakes_features %>% 
  select(depth,mag) 

### THIS IS HOW YOU RUN K-MEANS IN BASE R ###
normalized_fitted_kmeans <- kmeans(normalized_quakes_df, centers = n_clusters, iter.max = 10, nstart = 5)
normalized_fitted_kmeans
```

```{r}
# Visualize the clusters
normalized_quakes_df <- 
  normalized_quakes_df %>%
  mutate(cluster = as.factor(normalized_fitted_kmeans$cluster))

g <-
  ggplot(normalized_quakes_df, aes(x = depth, y = mag, color = cluster)) +
  geom_point(alpha=0.4, size=2) +
  labs(
    title = paste("Clusters found by k-means with", n_clusters, "clusters"),
    subtitle = "Do you like the clusters?\n Do they make any sense to you?",
    x = "Earthquake depth (km)",
    y = "Earthquake magnitude (Richter scale)",
    color = "Cluster"
  ) +
  theme_bw(base_size=15) 

g
```

# Part 3: Finding optimal number of clusters

How do we choose k in k-means?

## The "elbow" method

The elbow method, also known as *total within sum of squares*, is a *heuristic* used to determine the *optimal number of clusters* for k-means clustering.

**Idea behind elbow method:**

-   run k-means clustering on dataset for a range of values of $k$ (number of clusters)
-   and for each $k$ calculate the sum of squared distances of each point from its closest centroid (SSE): the elbow point is the point on the plot of SSE against the number of clusters ($k$) where the change in SSE begins to level off, indicating that adding more clusters does not improve the model much.

The steps to perform the elbow method are:

-   Select a range of $k$ values, generally from 1 to 10 or the square root of the number of observations/data points in the dataset.
-   Run k-means for each $k$ value and calculate the SSE (sum of squared distances of each point from its closest centroid).
-   Plot the SSE for each $k$ value: the point on the plot where the SSE starts to decrease at a slower rate is the elbow point, and the corresponding number of clusters is the optimal value for $k$

![Source: [Dangeti, Pratap. 2017. Statistics for Machine Learning: Techniques for Exploring Supervised, Unsupervised, and Reinforcement Learning Models with Python and R. Birmingham, UK: Packt Publishing]{style="font-size:0.4em;"}](../../../2023/autumn-term/weeks/week08/figures/elbow_method.png){style="text-align:right;font-size:0.9em"}

## The silhouette score

<div>

-   The silhouette value/score measures how similar a point is to its own cluster (cohesion) compared to other clusters (separation)
-   Value between -1 and +1: a **high value** is desirable and indicates that the point is placed in the correct cluster. If many points have negative Silhouette values, it may indicate that we have created too many or too few clusters.
-   The Silhouette Value $s(i)$ for each data point $i$ is defined as follows: $s(i)=\frac{b(i)-a(i)}{\max(a(i), b(i))} \textrm{ if } |C_i|>1 \textrm{ and } s(i)=0 \textrm{ if } |C_i|=1$ (note that the silhouette value is 0 for clusters with one single element,)

where $a(i)$ is the measure of similarity of the point $i$ to its own cluster, measured as the average distance of $i$ from other points in the cluster i.e:

for each data point $i \in C_i$ (data point $i$ in cluster $C_i$), $a_i=\frac{1}{|C_i|-1}\sum_{j\in C_i,i\neq j} d(i,j)$

and $b(i)$ is the measure of dissimilarity of $i$ from points of other clusters i.e: <br/>for each data point $i \in C_i$, $b(i)=\min_{i\neq j}\frac{1}{|C_j|}\sum_{j\in C_j}d(i,j)$

$d(i,j)$ is the distance between points $i$ and $j$, with the Euclidean Distance generally used as the distance metric.

</div>

## Average silhouette method

<div>

-   Compute the average silhouette score for all data points for a range $k$ of clusters
-   Plot the average silhouette score for each $k$ value: the number of clusters with the highest silhouette value is the optimal number of clusters we are looking for

</div>

![Source: <https://uc-r.github.io/kmeans_clustering>](../../../2023/autumn-term/weeks/week08/figures/silhouette_method.png)

## Other approaches

-   Calinski — Harabasz Method:

    -   It relies on computing the Calinski—Harabasz index $CH$ i.e the ratio of between-cluster variance (i.e between cluster sum of squares (BCSS)) and the within-cluster variance (i.e within cluster sum of squares (WCSS)) for a range of values of $k$:

    $CH=\frac{BCSS/k-1}{WCSS/N-k}=\frac{\sum_{i=1}^k n_i||c_i-c||^2/(k-1)}{\sum_{i=1}^k\sum_{x\in C_i}||x-c_i||^2/N-k}$

    where $n_i$ and $c_i$ are the number of points and centroid of the $i$-th cluster respectively, $c$ is the global centroid, $N$ is the total number of points and ${x_1,...x_N}$ are the data points that compose the dataset.

    -   then we plot the index against $k$: A higher value of the CH index means the clusters are dense and well separated, though there is no "acceptable" cut-off value. We need to choose that solution which gives a peak or at least an abrupt elbow on the line plot of CH index

See [this link](https://medium.com/@ozturkfemre/unsupervised-learning-determination-of-cluster-number-be8842cdb11) or [this link](https://towardsdatascience.com/are-you-still-using-the-elbow-method-5d271b3063bd) for more approaches to determine $k$

## Time to practice

As before, we can use recipes to pre-process our data before applying k-means:

```{r}
quakes_raw <- quakes %>% select(-stations)
# Recipes can still be used for unsupervised models
recipe_kmeans <- recipe(~ depth + mag, data = quakes_raw) %>% 
                          step_normalize(depth,mag) %>%
                          prep() %>% 
                          bake(new_data = NULL)
```

```{r}
recipe_kmeans %>% summary()
```

### Finding the optimal number of cluster for the k-means

Let's plot the elbow plot

```{r}
fviz_nbclust(recipe_kmeans, kmeans, method = "wss")
```

:dart: **Action Point**:

1.  Can you now plot the the curve based on the silhouette score instead? What is the optimal number of clusters based on this new curve?

```{r}
fviz_nbclust(recipe_kmeans, kmeans, method = "silhouette")
```

Let's have a quick look at other indices:

```{r}
# here's how we compute the Calinski-Harabasz index for various values of k

k_values <- 2:10
ch_indices <- sapply(k_values, function(k) {
  kmeans_results <- kmeans(recipe_kmeans, centers = k)
  cluster.stats(dist(recipe_kmeans), kmeans_results$cluster)$ch
})
ch_indices

# And now we plot the index against k

ch_df <- data.frame(k = k_values, CH_Index = ch_indices)

# Create the ggplot
ggplot(ch_df, aes(x = k, y = CH_Index)) +
  geom_line() +  # Adds a line
  geom_point() +  # Adds points at each k
  labs(x = "Number of Clusters (k)", 
       y = "CH Index", 
       title = "CH Index vs. Number of Clusters") +
  theme_minimal()  # Clean theme
```

Instead of plotting curves, we can compute optimal indices directly using the `NbClust` library:

```{r}
#Optimal number of clusters according to silhouette method
nb_sil <- NbClust(recipe_kmeans,distance = "euclidean",min.nc=2, max.nc=9, method = "kmeans",index="silhouette")

#Optimal number of clusters according to Calinski-Harabasz method
nb_ch <- NbClust(recipe_kmeans,distance = "euclidean",min.nc=2, max.nc=9, method = "kmeans",index="ch")

#Optimal number of clusters according to majority of available indices (does not include elbow method!)
nb_all <- NbClust(recipe_kmeans,distance = "euclidean",min.nc=2, max.nc=9, method = "kmeans",index="all")

cat("The optimal number of clusters according the average silhouette method is: ",nb_sil$Best.nc[1],"\nAccording to the Calinski-Harabasz index method, it is:",nb_ch$Best.nc[1])
```

```{r}
nb_all$Best.nc[1,] #print optimal number of cluster per index
hist(nb_all$Best.nc[1,], 
     breaks = max(na.omit(nb_all$Best.nc[1,])), 
     col = "skyblue",            # Change the color of the bars
     border = "black",           # Color of the bar borders
     main = "Histogram of Optimal Number of Clusters per Index ",  # Title of the plot
     xlab = "Optimal number of cluster per index",  # Custom x-axis label
     ylab = "Frequency") 
```

#### What if I want to initialize my k-means with k-means++ instead of the usual random initialization?

Let's use 4 clusters as before.

```{r}
kmeans_result <- KMeans_rcpp(recipe_kmeans, clusters = 4, initializer = "kmeans++")
kmeansplus_res <- recipe_kmeans %>%
    mutate(cluster = as.factor(kmeans_result$clusters))
g <- ggplot(kmeansplus_res, aes(x = depth, y = mag, color = cluster)) +
  geom_point(alpha = 0.4, size = 2) +
  labs(
    title = paste("Clusters found by k-means++ with 4 clusters"),
    subtitle = "Do you like the clusters?\n Do they make any sense to you?",
    x = "Earthquake depth (km)",
    y = "Earthquake magnitude (Richter scale)",
    color = "Cluster"
  ) +
  theme_bw(base_size = 15)

g
```

#### What if I want to perform k-medoid on the quakes data instead?

You can still use the elbow method (with a slight tweak!):

```{r}
fviz_nbclust(recipe_kmeans, pam, method = "wss") #pam indicates that the method we're using to compute the within-sum-of-squares is k-medoid!
```

The optimal number of clusters according to this elbow plot is 4 (once again).

```{r}
#make this example reproducible
set.seed(1)

#perform k-medoids clustering with k = 4 clusters
kmed <- pam(recipe_kmeans, k = 4)
```

Let's visualise the k-medoid results

```{r}
fviz_cluster(kmed, data = recipe_kmeans, palette = "Set2", ggtheme = theme_minimal(), geom = "point")
```

## Part 4: What if I have way more than two dimensions?

For this part, we look at a dataset of 26 numeric indicators (some originating from [Worldwide Governance Indicators (WGI) project](https://info.worldbank.org/governance/wgi/) at the World Bank provides a single score for several different governance indicators, including adherence to the rule of law, how well designed regulations are and how effective government in each country is and some from [Freedom House](https://freedomhouse.org/) relating to political freedoms) across 170 countries. The data relates to either 2019 or (if this was unavailable) the latest year no earlier than 2016. These 170 countries are those with a complete dataset - any country missing any data is excluded from the analysis.

The dataset indicators has from six categories:

-   Economic fundamentals (e.g GDP, GDP per capita)

-   Demographics (e.g total population, population density, urban population (% of total population)

-   Freedom and Rights (e.g Civil liberties score, Political rights score)

-   Governance (e.g WGI-GovernanceEffectiveness, WGI - RegulatoryQuality)

-   Health and Education (e.g Literacy rate, adult total (% of people ages 15 and above), Current health expenditure (% of GDP))

-   Infrastructure (e.g Fixed broadband subscriptions (per 100 people))

The analysis that follows is inspired by [this post](this%20post) and the dataset was sourced from the GitHub repository linked to it.

Let's load the data

```{r}
# Consolidated indicators
development_indicators <- read_csv("data/development-indicators-latest.csv")

# Indicators grouped into categories
indicators <- read_csv("data/indicators-all.csv")
```

Here's a quick overview of the data:

```{r}
skim(development_indicators)
```

The chart below shows the distribution of each numeric variable we use in the analysis.

```{r, fig.height=15}
development_indicators %>% 
  pivot_longer(cols = where(is.numeric), names_to = 'name', values_to = 'value') %>% 
  ggplot(aes(value)) + 
  geom_density(fill = '#98A6D4') + 
  facet_wrap(~str_wrap(name, 30), scales = 'free', ncol = 3) + 
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(), 
        strip.text.x = element_text(size = 10))  + 
  labs(x = "", y = "", title = "Distribution of indicators")
```

### Dimensionality reduction with PCA

Many of the features in our data vary together (co-vary), and the variation in one variable is likely to be duplicated in another. Principal component analysis (PCA) is a way of identifying the ways in which numeric predictors co-vary.

Additionally, applying k-means on this many features is likely to be too computationally expensive (k-means is $O(n^2)$)!

We also saw that most of our variables had a distribution which was not a normal distribution. However, k-means performs best when the data's shape is spherical i.e when the underlying distribution is normal (and obviously when the scale of all features is the same)!

So, aside before reducing our data's dimensionality with PCA, we'll have to transform the data to ensure that the distributions of the variables in it are as close to a normal distribution as possible and then scale (i.e normalize) the data. The first goal is accomplished with a log transformation of the variables that have a skewed distribution (i.e `GDP (current US$)`, `Land area (sq. km)`, `Population, total` and `Population density (people per sq. km of land area)` ).

We then apply PCA.

```{r}
development_recipe <- recipe(~ ., data = development_indicators) %>% 
  update_role(country_code, country_name, region, sub_region, new_role = 'id') %>% 
  step_log(`GDP (current US$)`, 
           `Land area (sq. km)`,
           `Population, total`, 
           `Population density (people per sq. km of land area)`, base = 2) %>% 
  step_normalize(all_predictors()) %>% 
  step_pca(all_predictors())

development_prepped <- prep(development_recipe)

pca_fit <- development_prepped$steps[[3]]$res

# the eigenvalues determine the amount of variance explained
eigenvalues <-
  pca_fit %>%
  broom::tidy(matrix = "eigenvalues")

eigenvalues %>%
  filter(PC <= 30) %>%
  ggplot(aes(PC, cumulative)) +
  geom_col()+
  geom_text(aes(label = cumulative),angle = 90)
```

We can retain 77% of our data variance with 4 components. So let's keep 4 components.

Let's have a look at the most important variables to the first 4 principal components. The graph below shows which of the original variables contribute the most to the overall component scores.

```{r,fig.width=12, fig.height=14}
tidyied_pca <- tidy(development_prepped, 3)

tidyied_pca <- tidyied_pca %>% 
  inner_join(indicators, by = c("terms" = "indicator_name"))

tidyied_pca %>% 
  filter(str_extract(component, "(\\d)+") %>% as.numeric() <= 4) %>% 
  group_by(component) %>% 
  top_n(n = 8, wt = abs(value)) %>% 
  mutate(component = fct_inorder(component)) %>% 
  ungroup() %>% 
  mutate(terms = reorder_within(terms, abs(value), component)) %>% 
  ggplot(aes(abs(value), terms)) + 
  geom_col(aes(fill = indicator_category)) + 
  geom_vline(xintercept = 0) + 
  scale_y_reordered() + 
  guides(fill = guide_legend(title = "Category")) + 
  theme(axis.text = element_text(size = 9), legend.position = 'top') + 
  facet_wrap(~component, ncol = 2, scales = 'free') + 
  labs(x = "", y = "")
```

Governance indicators feature prominently in PC1, as well as infrastructure variables related to Internet connectivity. The number of people under the age of 40 also contributes heavily to this component. Given that PC1 explains almost half of the variance in the data, these are some of the most important predictors in the analysis to explaining variance between countries.

PC2 contains several variables related to absolute size or 'bigness', such as a country's population, economic output and land area. This component also features variables such as civil liberties and political rights scores of countries. This component account for (about) another 11% of the variance of our data.

This scatter plot shows the projection of each country into the first two principal components. We can see many European nations to the left hand side of the chart, these are countries which are likely to score well in the governance indicators, life expectancy, and Internet usage. The right hand side features many African nations, who will perform less well on these indicators.

As PC2 heavily relates to size, we can see towards the top of the visualisation many small nations (such as Grenada, Tonga and Kiribati) and large ones at the other end of the y axis (such as China, Russia and Brazil).

```{r}
development_bake <- bake(development_prepped, new_data = NULL)

development_bake %>% 
  ggplot(aes(PC1, PC2)) + 
  geom_point(aes(color = region)) + 
  geom_text(aes(label = country_name), check_overlap = TRUE)
```

### Clustering the countries with k-means

We can now the countries in the principal component space using k-means clustering.

Let's determine the optimal k

```{r}
set.seed(2021)

development_kmeans <- development_bake %>% 
  select(where(is.numeric))
fviz_nbclust(recipe_kmeans, kmeans, method = "wss")
```

```{r}
set.seed(42)
#Applying k-means
dvlpt_kmeans <- kmeans(development_kmeans, centers = n_clusters, iter.max = 10, nstart = 5)
dvlpt_kmeans
```

```{r}
development_bake <- 
  development_bake %>%
  mutate(cluster = as.factor(dvlpt_kmeans$cluster))
```

Let's redo the earlier plot per country, now colored with clusters:

```{r}
development_bake %>% 
  ggplot(aes(PC1, PC2, color = cluster)) + 
  geom_point() + 
  geom_text(aes(label = country_name), color ='black', check_overlap = TRUE) + 
  geom_rug()
```

Finally, let's give some names to the clusters and re-plot this time on a map!

```{r}
development_bake <- development_bake %>% 
  mutate(cluster_label = case_when(cluster == 1 ~ "Rich, well governed and free",
                                   cluster == 2 ~ "Poor, but young",
                                   cluster == 3 ~ "Small, (mostly) free and developing",
                                   cluster == 4 ~ "(mostly) big and developing"))

development_indicators <- development_indicators %>% 
  inner_join(development_bake %>% select(country_code, cluster_label), by = c("country_code"))

world <- ne_countries(scale = "small", returnclass = "sf") %>% 
  filter(geounit != 'Antarctica')

world_with_clusters <- world %>% 
  select("country_code" = gu_a3, admin) %>% 
  left_join(development_indicators)

world_with_clusters %>% 
  ggplot() + 
  geom_sf(aes(fill = factor(cluster_label)))
```

### Part 4 : K-Means is not always optimal!

Let's go back to [this page](https://www.naftaliharris.com/blog/visualizing-k-means-clustering/) and check the smiley face example. How does K-means perform?

Now let's look at the same example with [DBSCAN](https://www.naftaliharris.com/blog/visualizing-dbscan-clustering/)

#### Limitations of K-Means

-   Selecting k is not straightforward

-   It assumes the shape of the data is spherical

-   The initialization is random and can have an impact on the final result. It might be mitigated with K-means++

-   It is sensitive to outliers. This could be mitigated with the use of K-medoids instead of K-means.

    k-means needs the number of clusters `k` as a parameter. However, as we saw this week, selecting `k` is not exactly straightforward.

#### How does DBSCAN work?

![](./figures_W09/dbscan-principle.png)

Fig 1: Different types of points and parameters in DBSCAN

*A few definitions to start with*

-   Any point $x$ in the data set, with a neighbour count greater than or equal to `MinPts`, is a **core point**
-   $x$ is a **border point**, if the number of its neighbors is less than `MinPts`, but it belongs to the $\epsilon$-neighborhood of some core point $z$.
-   If a point is neither a core nor a border point, then it is called a **noise point** or an **outlier**.

E.g: The figure (Fig. 1) shows the different types of points (core, border and outlier points) using `MinPts = 6`:

-   $x$ is a core point because $neighbours_\epsilon(x)=6$

-   $y$ is a border point because $neighbours_\epsilon(y)<MinPts$, but it belongs to the $\epsilon$-neighborhood of the core point $x$.

-   Finally, $z$ is a noise point.

*A few more definitions:*

-   **Direct density reachable**: A point “A” is directly density reachable from another point “B” if: i) “A” is in the $\epsilon$-neighborhood of “B” and ii) “B” is a core point. E.g $y$ is directly density reachable
-   **Density reachable**: A point “A” is density reachable from “B” if there are a set of core points leading from “B” to “A.
-   **Density connected**: Two points “A” and “B” are density connected if there are a core point “C”, such that both “A” and “B” are density reachable from “C”.

In density-based clustering, a cluster is defined as a group of density connected points.

### The DBSCAN Algorithm

1.  For each point $x_i$, compute the distance between $x_i$ and the other points. Find all neighbor points within distance `eps` of the starting point ($x_i$). Each point, with a neighbor count greater than or equal to `MinPts`, is marked as core point or visited.

2.  For each core point, if it’s not already assigned to a cluster, create a new cluster. Find recursively all its density connected points and assign them to the same cluster as the core point.

3.  Iterate through the remaining unvisited points in the data set.

The points that do not belong to any cluster are treated as outliers or noise.

If this explanation confuses you, you can have a look at [this Youtube video](https://www.youtube.com/watch?app=desktop&v=RDZUdRSDOok)

### How would we apply DBSCAN to the quake example?

```{r}
df<-normalized_quakes_df %>% select(-cluster) %>% select(where(is.numeric))
db <- dbscan(df, eps = 0.26, MinPts = 4)

# Plot DBSCAN results
#plot(db, clustering_data, main = "DBSCAN results", frame = FALSE)
fviz_cluster(db, data = df, stand = FALSE,
             ellipse = FALSE, show.clust.cent = FALSE,
             geom = "point", ggtheme = theme_bw(base_size=20))
```

### How do we choose the `eps` parameter in DBSCAN

Now that we now how DBSCAN works, the million dollar question is how do we set its parameters? Well, once again, it's rather tricky...

The value of `minPts` should be at least one greater than the number of dimensions of the dataset, i.e., **`minPts>=Dimensions+1`**. It does not make sense to take `minPts` as 1 because the result will be each point being a separate cluster. Therefore, it must be at least 3. Generally, we go for at least twice the dimensions, e.g if we only only select two variables in our data (i.e `depth` and `mag` ), we'd select `minPts` as 4.

Now that `minPts` is set, we can select `eps` using the k-distance graph. In practice, this means that:

-   you compute the average distance (by default Euclidean distance) between each point and its k nearest neighbors, where k = the `minPts` value you selected.

-   the average k-distances are then plotted in ascending order on a k-distance graph. You’ll find the optimal value for `eps` at the point of maximum curvature (i.e. where the graph has the greatest slope).

    The code below gives an example of the k-distance graph for our ongoing quakes clustering example...

```{r}
dbscan::kNNdistplot(df, k = 4)+
abline(h = 0.26, lty = 2) +
scale_y_continuous(n.breaks = 40)     
```

As we can see, the point of maximum curvature is around 25, which sets our `eps` value...

### Quick note: Similarity metrics

Underlying the discussion about clustering, are similarity metrics. We mostly assumed until now that the similarity metric that was used was the Euclidean distance (and for k-means, this is a fair assumption to make). But there are cases where the Euclidean distance would not be an appropriate choice (e.g when comparing strings or sets or for computer vision or recommendation engines applications).

Here is a quick run-through of some of the most common similarity metrics:

*Real-valued data:*

-   Euclidian distance: $d(A,B)=\sqrt{\sum_{i=1}^n |A_i-B_i|^2}$

-   Cosine distance: $S_C(A,B)=\frac{\sum_{i=1}^n A_iB_i}{\sqrt{\sum_{i=1}^n A_i^2}\sqrt{\sum_{i=1}^n B_i^2}}$

where **A** and **B** are two n-dimensional vectors of attributes, and $A_i$ and $B_i$ the $i$th components of vectors **A** and **B** respectively

(metric frequently used in NLP and computer vision - see \[this page\](<https://www.datastax.com/guides/real-world-applications-of-cosine-similarity>) for more information)

*String similarity:*

-   Hamming distance: given two strings of equal length, it is the number of positions at which the corresponding symbols (i.e either letters, bits, or decimal digits) are different e.g:

    -   [0000]{style="color:#BA160C"} and [1111]{style="color:#BA160C"} is 4
    -   ca[t]{style="color:#BA160C"} and ca[r]{style="color:#BA160C"} is 1
    -   c[a]{style="color:#BA160C"}l[l]{style="color:#BA160C"} and c[o]{style="color:#BA160C"}l[d]{style="color:#BA160C"} is 2

-   Levenstein distance: given two strings, it's the minimal number of edits (substitutions, deletions, insertions) needed to turn one string into another e.g:

    The Levenshtein distance between *kitten* and *sitting* is 3, since the following 3 edits change one into the other, and there is no way to do it with fewer than 3 edits:

    1.  [k]{style="color:#BA160C"}itten → [s]{style="color:#BA160C"}itten (substitution of "s" for "k"),
    2.  sitt[e]{style="color:#BA160C"}n → sitt[i]{style="color:#BA160C"}n (substitution of "i" for "e"),
    3.  sittin → sittin[g]{style="color:#BA160C"} (insertion of "g" at the end).

*Set similarity:*

A set is a collection of items with no order or repetition. Sets are typically used to represent relationships or associations between objects or even people.

-   Jaccard index: it is often used in recommendation systems and social media analysis and measures the similarity between two sets based on the number of items present in both sets relative to the total number of items:

    $J(A,B)=\frac{|A\cap B|}{|A\cup B|}$

    where A and B are two sets, $|A\cap B|$ the number of items included in both sets and $|A\cup B|$ the total number of items

Another useful resources: \[this page\](<https://scikit-learn.org/stable/modules/clustering.html>) compares various clustering algorithms (among which k-means and DBSCAN) and shows their parameters, use cases,...

**Bibliography**

*About k-means++*

David Arthur and Sergei Vassilvitskii. 2007. *K-means++: the advantages of careful seeding*. In Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms (SODA '07). Society for Industrial and Applied Mathematics, USA, 1027–1035. <https://dl.acm.org/doi/10.5555/1283383.1283494>

*About k-medoid*

Leonard Kaufman and Peter J. Rousseeuw (1990). *Partitioning Around Medoids (Program PAM)*. In Finding Groups in Data (eds L. Kaufman and P.J. Rousseeuw). <https://doi.org/10.1002/9780470316801.ch2>

Erich Schubert, Peter J. Rousseeuw, *Fast and eager k-medoids clustering: O(k) runtime improvement of the PAM, CLARA, and CLARANS algorithms*, Information Systems, Volume 101, 2021, 101804, ISSN 0306-4379, <https://doi.org/10.1016/j.is.2021.101804.>

*About DBSCAN*

Martin Ester, Hans-Peter Kriegel, Jörg Sander, and Xiaowei Xu. 1996. *A density-based algorithm for discovering clusters in large spatial databases with noise*. In Proceedings of the Second International Conference on Knowledge Discovery and Data Mining (KDD'96). AAAI Press, 226–231. <https://dl.acm.org/doi/10.5555/3001460.3001507>

*About missing data mechanisms and missing data handling*

Enders, Craig K. *Applied missing data analysis*. Guilford Publications, 2022.

Scheffer, J. (2002), Dealing with missing data, Research Letters in the Information and Mathematical Sciences, 3, 153-160. <http://hdl.handle.net/10179/4355>

McKnight, Patrick E., et al. *Missing data: A gentle introduction*. Guilford Press, 2007.
