---
title: "Clustering"
subtitle: "Hierarchical clustering"
date: "June 15th, 2022"
---

```{r setup, echo = FALSE}
knitr::opts_chunk$set(fig.retina = 3, warning = FALSE, message = FALSE)
```

```{r xaringan-themer, include=FALSE, warning=FALSE}
library(xaringanthemer)
style_mono_accent(base_color = "#cc002b")
```


## Prep NBA player dataset

Created dataset of NBA player statistics per 100 possessions using [`ballr`](https://cran.r-project.org/web/packages/ballr/vignettes/use-ballr.html)

```{r load-data, warning = FALSE, message = FALSE}
library(tidyverse)
nba_pos_stats <- 
  read_csv("http://www.stat.cmu.edu/cmsac/sure/2022/materials/data/sports/clustering/nba_2022_player_per_pos_stats.csv")
# Find rows for players indicating a full season worth of stats
tot_players <- nba_pos_stats %>% filter(tm == "TOT")
# Stack this dataset with players that played on just one team
nba_player_stats <- nba_pos_stats %>% filter(!(player %in% tot_players$player)) %>% 
  bind_rows(tot_players)
# Filter to only players with at least 125 minutes played
nba_filtered_stats <- nba_player_stats %>% filter(mp >= 125)
head(nba_filtered_stats)
```


---

## Let's work from the bottom-up...

.pull-left[

- __Review__: We have $p$ variables for $n$ observations $x_1,\dots,x_n$,

- Compute the __distance / dissimilarity__ between observations

- e.g. __Euclidean distance__ between observations $i$ and $j$

$$d(x_i, x_j) = \sqrt{(x_{i1}-x_{j1})^2 + \cdots + (x_{ip}-x_{jp})^2}$$

__What are the distances between these NBA players using `x3pa` and `trb`?__ 

```{r nba-start-plot, eval = FALSE, echo = FALSE}
nba_filtered_stats %>%
  ggplot(aes(x = x3pa, y = trb)) +
  geom_point(alpha = 0.5) + 
  theme_bw() +
  coord_fixed()
```


]
.pull-right[

```{r ref.label='nba-start-plot', echo = FALSE}

```

]

---

## Remember to standardize!


.pull-left[

```{r nba-std-plot, eval = FALSE}
nba_filtered_stats <- nba_filtered_stats %>%
  mutate(std_x3pa = as.numeric(scale(x3pa)),
         std_trb = as.numeric(scale(trb)))
nba_filtered_stats %>%
  ggplot(aes(x = std_x3pa, y = std_trb)) +
  geom_point(alpha = 0.5) + 
  theme_bw() +
  coord_fixed()
```


]
.pull-right[

```{r ref.label='nba-std-plot', echo = FALSE}

```

]

---

## Compute the distance matrix using `dist()`

.pull-left[

- Compute pairwise Euclidean distance

```{r compute-dist}
player_dist <- dist(dplyr::select(nba_filtered_stats, #<<
                                  std_x3pa, std_trb))
```

- Returns an object of `dist` class - i.e., not a matrix

- Can convert to a matrix, then set the row and column names:

```{r dist-matrix}
player_dist_matrix <- as.matrix(player_dist) #<<
rownames(player_dist_matrix) <- nba_filtered_stats$player
colnames(player_dist_matrix) <- nba_filtered_stats$player
head(player_dist_matrix[1:3, 1:3])
```


]
.pull-right[

Can convert to a long table for plotting with `ggplot`:

```{r tidy-dist-plot, eval = FALSE}
long_dist_matrix <- 
  as_tibble(player_dist_matrix) %>%
  mutate(player1 = rownames(player_dist_matrix)) %>%
  pivot_longer(cols = -player1, #<<
               names_to = "player2", #<<
               values_to = "distance") #<<
long_dist_matrix %>%
  ggplot(aes(x = player1, y = player2, 
             fill = distance)) +
  geom_tile() +
  theme_bw() +
  theme(axis.text = element_blank(), 
        axis.ticks = element_blank(),
        legend.position = "bottom") +
  scale_fill_gradient(low = "darkorange", 
                      high = "darkblue")
```


]

---

## This is useless...

```{r ref.label='tidy-dist-plot', echo = FALSE, fig.align='center'}

```

---

### Code interlude: arrange your heatmap with [`seriation`](https://github.com/mhahsler/seriation)

.pull-left[
```{r seriate-plot, eval = FALSE}
library(seriation)
player_dist_seriate <- seriate(player_dist) #<<
player_order <- get_order(player_dist_seriate) #<<
player_names_order <- 
  nba_filtered_stats$player[player_order]

long_dist_matrix %>%
  mutate(player1 = 
           fct_relevel(player1, 
                       player_names_order),
         player2 = 
           fct_relevel(player2, 
                       player_names_order)) %>%
  ggplot(aes(x = player1, y = player2, 
             fill = distance)) +
  geom_tile() + theme_bw() +
  theme(axis.text = element_blank(), 
        axis.ticks = element_blank(),
        legend.position = "bottom") +
  scale_fill_gradient(low = "darkorange",
                      high = "darkblue")
```

]
.pull-right[
```{r ref.label='seriate-plot', echo = FALSE, fig.align='center'}

```
]

---

## (Agglomerative) [Hierarchical clustering](https://en.wikipedia.org/wiki/Hierarchical_clustering)

Let's pretend all $n$ observations are in their own cluster

--

- Step 1: Compute the pairwise dissimilarities between each cluster

  - e.g., distance matrix on previous slides
  
--
  
- Step 2: Identify the pair of clusters that are __least dissimilar__

--

- Step 3: Fuse these two clusters into a new cluster!

--

- __Repeat Steps 1 to 3 until all observations are in the same cluster__

--

__"Bottom-up"__, agglomerative clustering that forms a __tree / hierarchy__ of merging

No mention of any randomness!

No mention of the number of clusters $K$!

---

## (Agglomerative) [Hierarchical clustering](https://en.wikipedia.org/wiki/Hierarchical_clustering)

.pull-left[
Start with all observations in their own cluster

- Step 1: Compute the pairwise dissimilarities between each cluster

- Step 2: Identify the pair of clusters that are __least dissimilar__

- Step 3: Fuse these two clusters into a new cluster!

- __Repeat Steps 1 to 3 until all observations are in the same cluster__
]
.pull-right[
```{r out.width='70%', echo = FALSE, fig.align='center'}
knitr::include_graphics("https://upload.wikimedia.org/wikipedia/commons/thumb/b/b5/Clusters.svg/250px-Clusters.svg.png")
```
]

---

## (Agglomerative) [Hierarchical clustering](https://en.wikipedia.org/wiki/Hierarchical_clustering)

.pull-left[
Start with all observations in their own cluster

- Step 1: Compute the pairwise dissimilarities between each cluster

- Step 2: Identify the pair of clusters that are __least dissimilar__

- Step 3: Fuse these two clusters into a new cluster!

- __Repeat Steps 1 to 3 until all observations are in the same cluster__
]
.pull-right[
```{r out.width='85%', echo = FALSE, fig.align='center'}
knitr::include_graphics("https://upload.wikimedia.org/wikipedia/commons/thumb/a/ad/Hierarchical_clustering_simple_diagram.svg/418px-Hierarchical_clustering_simple_diagram.svg.png")
```

Forms a __dendrogram__ (typically displayed from bottom-up)
]


---

## How do we define dissimilarity between clusters?

We know how to compute distance / dissimilarity between two observations

__But how do we handle clusters?__

  - Dissimilarity between a cluster and an observation, or between two clusters
  
--

We need to choose a __linkage function__! Clusters are built up by __linking them together__

--

Compute all pairwise dissimilarities between observations in cluster 1 with observations in cluster 2

i.e. Compute the distance matrix between observations, $d(x_i, x_j)$ for $i \in C_1$ and $j \in C_2$

--

  - __Complete linkage__: Use the __maximum__ value of these dissimilarities: $\underset{i \in C_1, j \in C_2}{\text{max}} d(x_i, x_j)$
  
--

  - __Single linkage__: Use the __minimum__ value: $\underset{i \in C_1, j \in C_2}{\text{min}} d(x_i, x_j)$

--

  - __Average linkage__: Use the __average__ value: $\frac{1}{|C_1| \cdot |C_2|} \sum_{i \in C_1} \sum_{j \in C_2} d(x_i, x_j)$
  
--

Define dissimilarity between two clusters __based on our initial dissimilarity matrix between observations__

---

## Complete linkage example 

.pull-left[

- Use the `hclust` function with a `dist()` object

- Uses `complete` linkage by default

```{r nba-complete}
nba_complete_hclust <- 
  hclust(player_dist, method = "complete") #<<
```

- Need to use `cutree()` to return cluster labels:

```{r nba-complete-plot, eval = FALSE}
nba_filtered_stats %>%
  mutate(player_clusters = 
           as.factor(cutree(nba_complete_hclust, #<<
                            k = 4))) %>% #<<
  ggplot(aes(x = std_x3pa, y = std_trb,
             color = player_clusters)) +
  geom_point(alpha = 0.5) + 
  ggthemes::scale_color_colorblind() +
  theme_bw() +
  theme(legend.position = "bottom")
```


]
.pull-right[

Returns _compact_ clusters, similar to $K$-means

```{r ref.label = 'nba-complete-plot', echo = FALSE, fig.height=6, fig.align='center'}

```
]

---

## What are we cutting? Dendrograms

.pull-left[
Use the [`ggdendro`](https://cran.r-project.org/web/packages/ggdendro/index.html) package (instead of `plot()`)

```{r complete-dendro, eval = FALSE}
library(ggdendro)
ggdendrogram(nba_complete_hclust, theme_dendro = FALSE, #<<
             labels = FALSE, leaf_labels = FALSE) + #<< 
  labs(y = "Dissimilarity between clusters") +
  theme_bw() +
  theme(axis.text.x = element_blank(), 
        axis.title.x = element_blank(),
        axis.ticks.x = element_blank(),
        panel.grid = element_blank())
```

- Each __leaf__ is one observation

- __Height of branch indicates dissimilarity between clusters__

  - (After first step) Horizontal position along x-axis means nothing

]
.pull-right[
```{r ref.label = 'complete-dendro', echo = FALSE, fig.height=6, fig.align='center'}

```
]

---

## [Textbook example](https://bradleyboehmke.github.io/HOML/hierarchical.html)

```{r out.width='100%', echo = FALSE, fig.align='center'}
knitr::include_graphics("https://bradleyboehmke.github.io/HOML/19-hierarchical_files/figure-html/comparing-dendrogram-to-distances-1.png")
```

---

## Cut dendrograms to obtain cluster labels

.pull-left[
Specify the height to cut with `h` instead of `k`

```{r complete-dendro-cut, echo = FALSE, out.width='80%'}
ggdendrogram(nba_complete_hclust, theme_dendro = FALSE, #<<
             labels = FALSE, leaf_labels = FALSE) + #<< 
  labs(y = "Dissimilarity between clusters") +
  theme_bw() +
  theme(axis.text.x = element_blank(), 
        axis.title.x = element_blank(),
        axis.ticks.x = element_blank(),
        panel.grid = element_blank()) +
  geom_hline(yintercept = 6, linetype = "dashed",
             color = "darkred")
```

]
.pull-right[
```{r nba-h-cut, eval = FALSE}
cutree(nba_complete_hclust, h = 6) #<<
```

```{r nba-complete-cut-plot, echo = FALSE, out.width='80%'}
nba_filtered_stats %>%
  mutate(player_clusters = 
           as.factor(cutree(nba_complete_hclust, #<<
                            h = 6))) %>% #<<
  ggplot(aes(x = std_x3pa, y = std_trb,
             color = player_clusters)) +
  geom_point(alpha = 0.5) + 
  ggthemes::scale_color_colorblind() +
  theme_bw() +
  theme(legend.position = "bottom")
```
]

---

## Single linkage example

.pull-left[
Change the `method` argument to `single`

```{r single-dendro-cut, echo = FALSE, out.width='80%'}
nba_single_hclust <- 
  hclust(player_dist, method = "single")
ggdendrogram(nba_single_hclust, theme_dendro = FALSE, #<<
             labels = FALSE, leaf_labels = FALSE) + #<< 
  labs(y = "Dissimilarity between clusters") +
  theme_bw() +
  theme(axis.text.x = element_blank(), 
        axis.title.x = element_blank(),
        axis.ticks.x = element_blank(),
        panel.grid = element_blank())
```

]
.pull-right[
Results in a __chaining__ effect

```{r nba-single-plot, echo = FALSE, out.width='80%'}
nba_filtered_stats %>%
  mutate(player_clusters = 
           as.factor(cutree(nba_single_hclust, #<<
                            k = 4))) %>% #<<
  ggplot(aes(x = std_x3pa, y = std_trb,
             color = player_clusters)) +
  geom_point(alpha = 0.5) + 
  ggthemes::scale_color_colorblind() +
  theme_bw() +
  theme(legend.position = "bottom")
```
]

---

## Average linkage example

.pull-left[
Change the `method` argument to `average`

```{r average-dendro-cut, echo = FALSE, out.width='80%'}
nba_average_hclust <- 
  hclust(player_dist, method = "average")
ggdendrogram(nba_average_hclust, theme_dendro = FALSE, #<<
             labels = FALSE, leaf_labels = FALSE) + #<< 
  labs(y = "Dissimilarity between clusters") +
  theme_bw() +
  theme(axis.text.x = element_blank(), 
        axis.title.x = element_blank(),
        axis.ticks.x = element_blank(),
        panel.grid = element_blank())
```

]
.pull-right[
Closer to `complete` but varies in compactness

```{r nba-average-plot, echo = FALSE, out.width='80%'}
nba_filtered_stats %>%
  mutate(player_clusters = 
           as.factor(cutree(nba_average_hclust, #<<
                            k = 4))) %>% #<<
  ggplot(aes(x = std_x3pa, y = std_trb,
             color = player_clusters)) +
  geom_point(alpha = 0.5) + 
  ggthemes::scale_color_colorblind() +
  theme_bw() +
  theme(legend.position = "bottom")
```
]

---

## More linkage functions


- __Centroid linkage__: Computes the dissimilarity between the centroid for cluster 1 and the centroid for cluster 2

  - i.e. distance between the averages of the two clusters
  
  - use `method = centroid`
  
--

- __Ward’s linkage__: Merges a pair of clusters to minimize the within-cluster variance

  - i.e. aim is to minimize the objection function from $K$-means
  
  - can use `ward.D` or `ward.D2` (different algorithms)
  
--

```{r out.width='30%', echo = FALSE, fig.align='center'}
knitr::include_graphics("https://media1.tenor.com/images/bfb8e3e881ac4413ae12b61c4b982d60/tenor.gif?itemid=5140031")
```

---

## [Minimax linkage](http://statweb.stanford.edu/~tibs/sta306bfiles/minimax-clustering.pdf)

- Each cluster is defined __by a prototype__ observation (most representative)

- __Identify the point whose farthest point is closest__ (hence the minimax)

```{r out.width='60%', echo = FALSE, fig.align='center'}
knitr::include_graphics("https://europepmc.org/articles/PMC4527350/bin/nihms637357f2.jpg")
```

- Use this minimum-maximum distance as the measure of cluster dissimilarity

- Dendogram interpretation: each point point is $\leq h$ in dissimilarity to the __prototype__ of cluster

- __Cluster centers are chosen among the observations themselves - hence prototype__

---

## Minimax linkage example

.pull-left[

- Easily done in `R` via the [`protoclust`](https://github.com/jacobbien/protoclust) package

- Use the `protoclust()` function to apply the clustering to the `dist()` object

```{r nba-minimax, eval = FALSE}
library(protoclust)
nba_minimax <- protoclust(player_dist) #<<
ggdendrogram(nba_minimax, 
             theme_dendro = FALSE, 
             labels = FALSE, 
             leaf_labels = FALSE) + 
  labs(y = "Maximum dissimilarity from prototype") +
  theme_bw() +
  theme(axis.text.x = element_blank(), 
        axis.title.x = element_blank(),
        axis.ticks.x = element_blank(),
        panel.grid = element_blank())
```


]
.pull-right[

```{r ref.label = 'nba-minimax', echo = FALSE, fig.height=6, fig.align='center'}

```
]

---

## Minimax linkage example

.pull-left[

- Use the `protocut()` function to make the cut 

- But then access the cluster labels `cl`

```{r nba-minimax-cut, eval = FALSE}
minimax_player_clusters <- 
  protocut(nba_minimax, k = 4) #<<
nba_filtered_stats %>%
  mutate(player_clusters = 
           as.factor(minimax_player_clusters$cl)) %>% #<<
  ggplot(aes(x = std_x3pa, y = std_trb,
             color = player_clusters)) +
  geom_point(alpha = 0.5) + 
  ggthemes::scale_color_colorblind() +
  theme_bw() +
  theme(legend.position = "bottom")
```


]
.pull-right[

```{r ref.label = 'nba-minimax-cut', echo = FALSE, fig.height=6, fig.align='center'}

```
]

---

## Minimax linkage example

- Want to check out the prototypes for the three clusters

- `protocut` returns the indices of the prototypes (in order of the cluster labels)

```{r prototypes-indices}
minimax_player_clusters$protos
```

- View these player rows using `slice`:

```{r proto-players}
nba_filtered_stats %>%
  dplyr::select(player, pos, age, std_x3pa, std_trb) %>%
  slice(minimax_player_clusters$protos)
```

---

## Wrapping up...

- For context, how does player position (`pos`) relate to our clustering results?

```{r minimax-comparison}
table("Clusters" = minimax_player_clusters$cl, "Positions" = nba_filtered_stats$pos)
```

--

- Can see positions tend to fall within particular clusters...

- _What's the way to visually compare the two labels?_

--

- __We can easily include more variables__ - just changes our distance matrix

--

- But we might want to explore __soft__ assignments instead...



