```{r, include = FALSE}
source("../setup.R")
```

## Overview

We will cover:

* Defining distance measure
* $k$-means algorithm
* Hierarchical algorithms
* Making and using dendrograms

## Cluster analysis

<br>

- The aim of cluster analysis is to [group cases (objects) according to their similarity]{.monash-blue2} on the variables. It is also often called unsupervised classification, meaning that classification is the ultimate goal, but the classes (groups) are not known ahead of time. 
- Hence the first task in cluster analysis is to construct the class information. To determine closeness we start with [measuring the interpoint distances]{.monash-blue2}.

## Cluster these!

:::: {.columns}
::: {.column}

```{r}
#| message: false
#| echo: false
#| out-width: 100%
#| fig-width: 6
#| fig-height: 6

stdd <- function(x) (x-mean(x))/sd(x)
set.seed(20230513)
d_sep_cl <- matrix(rnorm(99*2), ncol=2)
d_sep_cl[1:33,1] <-
  d_sep_cl[1:33,1]+8
d_sep_cl[34:66,2] <-
  d_sep_cl[34:66,2]+8
d_sep_cl[34:66,1] <-
  d_sep_cl[34:66,1]+4
d_sep_cl <- data.frame(x1=stdd(d_sep_cl[,1]), 
                       x2=stdd(d_sep_cl[,2]))

x <- (runif(20)-0.5)*4
y <- x
d_nuis_pts <- data.frame(x1 = stdd(c(rnorm(50, -3), 
                            rnorm(50, 3), x)),
                 x2 = stdd(c(rnorm(50, -3), 
                            rnorm(50, 3), y)))

d_nuis_vars <- matrix(rnorm(99*2), ncol=2)
d_nuis_vars[1:49,1] <- d_nuis_vars[1:49,1]+8
d_nuis_vars <-
  data.frame(x1=stdd(d_nuis_vars[,1]),
             x2=stdd(d_nuis_vars[,2]))

d_odd_shapes <- matrix(rnorm(99*2),ncol=2)
d_odd_shapes[1:66,2] <- (d_odd_shapes[1:66,1])^2-5 + rnorm(66)*0.6
d_odd_shapes[1:66,1] <- d_odd_shapes[1:66,1]*3
d_odd_shapes <-
  data.frame(x1=stdd(d_odd_shapes[,1]),
             x2=stdd(d_odd_shapes[,2]))

p1 <- ggplot(d_sep_cl, aes(x=x1, y=x2)) + 
  geom_point(colour="#3B99B1", alpha=0.7) + 
  annotate("text", -2.5, 2.5, label="a") +
  xlim(-2.8, 2.8) + ylim(-2.8, 2.8) +
    theme(
        axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank())
p2 <- ggplot(d_nuis_pts, aes(x=x1, y=x2)) + 
  geom_point(colour="#3B99B1", alpha=0.7) + 
  annotate("text", -2.5, 2.5, label="b") +
  xlim(-2.8, 2.8) + ylim(-2.8, 2.8) +
    theme(
        axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank())
p3 <- ggplot(d_nuis_vars, aes(x=x1, y=x2)) + 
  geom_point(colour="#3B99B1", alpha=0.7) + 
  annotate("text", -2.5, 2.5, label="c") +
  xlim(-2.8, 2.8) + ylim(-2.8, 2.8) +
    theme(
        axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank())
p4 <- ggplot(d_odd_shapes, aes(x=x1, y=x2)) + 
  geom_point(colour="#3B99B1", alpha=0.7) + 
  annotate("text", -2.5, 2.5, label="d") +
  xlim(-2.8, 2.8) + ylim(-2.8, 2.8) +
    theme(
        axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank())
print(p1 + p2 + p3 + p4 + plot_layout(ncol=2))

```
:::
::: {.column style="font-size: 80%;"}
::: {.fragment}
```{r}
#| message: false
#| echo: false
#| out-width: 70%
#| fig-width: 6
#| fig-height: 6
set.seed(515)
km1 <- kmeans(d_sep_cl, 3)
d_sep_cl <- d_sep_cl |> 
  mutate(cl = factor(km1$cluster))

km2 <- kmeans(d_nuis_pts, 2)
d_nuis_pts <- d_nuis_pts |> 
  mutate(cl = factor(km2$cluster))

km3 <- kmeans(d_nuis_vars, 2)
d_nuis_vars <- d_nuis_vars |> 
  mutate(cl = factor(km3$cluster))

km4 <- kmeans(d_odd_shapes, 2)
d_odd_shapes <- d_odd_shapes |> 
  mutate(cl = factor(km4$cluster))

kp1 <- ggplot(d_sep_cl, aes(x=x1, y=x2, colour=cl)) + 
  geom_point(alpha=0.7) + 
  scale_color_discrete_divergingx(palette="Zissou 1") +
  annotate("text", -2.5, 2.5, label="a") +
  xlim(-2.8, 2.8) + ylim(-2.8, 2.8) +
  theme(legend.position = "none",
        axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank())
kp2 <- ggplot(d_nuis_pts, aes(x=x1, y=x2, colour=cl)) + 
  geom_point(alpha=0.7) + 
  scale_color_discrete_divergingx(palette="Zissou 1") +
  annotate("text", -2.5, 2.5, label="b") +
  xlim(-2.8, 2.8) + ylim(-2.8, 2.8) +
  theme(legend.position = "none",
        axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank())
kp3 <- ggplot(d_nuis_vars, aes(x=x1, y=x2, colour=cl)) + 
  geom_point(alpha=0.7) + 
  scale_color_discrete_divergingx(palette="Zissou 1") +
  annotate("text", -2.5, 2.5, label="c") +
  xlim(-2.8, 2.8) + ylim(-2.8, 2.8) +
  theme(legend.position = "none",
        axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank())
kp4 <- ggplot(d_odd_shapes, aes(x=x1, y=x2, colour=cl)) + 
  geom_point(alpha=0.7) + 
  scale_color_discrete_divergingx(palette="Zissou 1") +
  annotate("text", -2.5, 2.5, label="d") +
  xlim(-2.8, 2.8) + ylim(-2.8, 2.8) +
  theme(legend.position = "none",
        axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank())
print(kp1 + kp2 + kp3 + kp4 + plot_layout(ncol=2))
```

It's *easy* if we can [see the clusters]{.monash-blue2}, but what an [algorithm sees]{.monash-orange2} might be quite different. 

:::
:::
::::

## Seeing the clusters using spin-and-brush

```{r}
#| echo: true
#| eval: false
library(detourr)
set.seed(645)
detour(p_std[,2:5], 
       tour_aes(projection = bl:bm)) |>
       tour_path(grand_tour(2), fps = 60, 
                 max_bases=40) |>
       show_scatter(alpha = 0.7, 
                    axes = FALSE)
```


## How do you measure "close"? {.transition-slide .center}

## Common interpoint distance measures 

:::: {.columns}
::: {.column style="font-size: 70%;"}

Let $A=(x_{a1}, x_{a2}, ..., x_{ap})$ and $B=(x_{b1}, x_{b2}, ..., x_{bp})$.

[Euclidean]{.monash-blue2}

\begin{align*}
d_{EUC}(A, B) &= \sqrt{\sum_{j=1}^p (x_{aj}-x_{bj})^2} &\\
&= \sqrt{((X_A-X_B)^\top (X_A-X_B))}&
\end{align*}

[Other distance metrics]{.monash-blue2}

- Mahalanobis (or statistical) distance: $\sqrt{((X_A-X_B)^\top S^{-1} (X_A-X_B))}$

- Manhattan: $\sum_{j=1}^p|(X_{aj}-X_{bj})|$

- Minkowski: $(\sum_{j=1}^p|(X_{aj}-X_{bj})|^m)^{1/m}$

::: 

::: {.column style="font-size: 60%;"}

[Count data]{.monash-blue2}

- Canberra: $\frac{1}{n_z}\sum_{j=1}^p\frac{X_{aj}-X_{bj}}{X_{aj}+X_{bj}}$

- Bray-Curtis: $\frac{\sum_{j=1}^p|X_{aj}-X_{bj}|}{\sum_{j=1}^p(X_{aj}+X_{bj})}$

[Categorical variables]{.monash-blue2}

- 1- simple matching coefficient: $1-(\text{#matches})/p$
- Convert to dummy variables, and use Euclidean distance

[Mixed variable types]{.monash-blue2}

- Gower's distance

[Other]{.monash-blue2}

- Hamming: all binary variables, number of variables at which values are different.

- Cosine: $\frac{\sum_{j=1}^p X_{aj}X_{bj}}{||X_a|| ||X_b||}$ (How does this compare to a calculation of correlation??)

:::
::::

## Distance calculations

:::: {.columns}
::: {.column}

```{r}
#| message: FALSE
#| warning: false
#| echo: false
#| fig-width: 4
#| fig-height: 4
#| out-width: 70%
x <- data.frame(V1 = c(7.3, 7.4, 4.1),
                    V2 = c(7.6, 7.2, 4.6),
                    V3 = c(7.7, 7.3, 4.6),
                    V4 = c(8.0, 7.2, 4.8),
                    point = factor(c("a1", "a2", "a3")))
x
pscat <- ggpairs(x, columns=1:4,
                 upper=list(continuous="points"),
                 diag=list(continuous="blankDiag"),
                 axisLabels="internal",
                 ggplot2::aes(colour=point)) +
    scale_colour_discrete_divergingx(
      palette = "Zissou 1", nmax=4) +
    xlim(3.7, 8.5) + ylim(3.7, 8.5) + 
    theme_minimal() +
    theme(aspect.ratio=1)
pscat
```
:::
::: {.column}

Compute Euclidean distance between $a_1$ and $a_2$.

<br><br><br>

::: {.fragment}
🔑 [Standardise your variables!!!!]{.monash-blue2}

:::
::: {.fragment}
Could you compute a correlation distance? $d_{cor} = 1-|r|$
Is $a_1$ close to $a_3$ than $a_2$?
:::
:::
::::

## Basic rules of a distance metric

:::: {.columns}
::: {.column}

Anything can be a distance if it follows these rules:

1. $d(A, B) \geq 0$
2. $d(A, A) = 0$
3. $d(A, B) = d(B, A)$
4. Metric dissimilarity satisfies 
$d(A, B) \leq d(A, C) + d(C, B)$
:::
::: {.column style="font-size: 70%;"}
::: {.fragment}

<center>
![](../images/river-dist.png){width=450}
</center>

- If both points on left bank, or both on right bank, use Euclidean distance.
- If on opposite sides, Euclidean distances to bridge, plus length of bridge crossing.
- Does this satisfy the rules?
:::
:::
::::

## Dissimilarity vs similarity

- [Distance is a **dissimilarity**]{.monash-blue2} measure because small means close and large means far.

- [Correlation is a **similarity**]{.monash-blue2} measure because the larger the value the closer the objects. It can be converted to a dissimilarity with a transformation.

## $k$-means clustering {.transition-slide .center}

## k-means clustering - algorithm [(1/8)]{.smallest}

This is an [iterative]{.monash-blue2} procedure. To use it the number of clusters, $k$, must be decided first.  The stages of the iteration are:

1. [Initialize]{.monash-blue2} by either (a) partitioning the data into k groups, and compute the $k$ group means or (b) an [initial set of $k$ points]{.monash-blue2} as the first estimate of the [cluster means]{.monash-blue2} (seed points).
2. [Loop]{.monash-blue2} over all observations reassigning them to the group with the [closest mean]{.monash-blue2}.
3. [Recompute]{.monash-blue2} group [means]{.monash-blue2}.
4. Iterate steps 2 and 3 until [convergence]{.monash-blue2}.

[Thean C. Lim's blog post](https://theanlim.rbind.io/post/clustering-k-means-k-means-and-gganimate/)

## k-means clustering - algorithm [(2/8)]{.smallest}

:::: {.columns}
::: {.column}

```{r}
#| echo: false
set.seed(8) #6
df <- tibble(lbl=letters[1:12], 
             x1=sample(1:10, 12, replace=TRUE),
             x2=sample(1:10, 12, replace=TRUE))
df[1:4,2] <- df[1:4,2] + 12
df[5:8,3] <- df[5:8,3] + 12
kable(df) %>%
  kable_styling("striped", position = "center", 
                row_label_position = "c", 
                font_size=24) %>%
  row_spec(0, color = "white", background = "#3F9F7A") %>%
  column_spec(1:3, border_right=TRUE, width="2cm") 
```

:::
::: {.column}

```{r out.width="100%", fig.width=4, fig.height=4}
#| echo: false
ggplot(data=df, aes(x1, x2)) + geom_text(aes(label=lbl)) + 
  xlab("") + ylab("") + 
  xlim(c(1,22)) + ylim(c(1,22)) 
```
:::
::::

## k-means clustering - algorithm [(3/8)]{.smallest}

:::: {.columns}
::: {.column}

```{r}
#| echo: false
# Initial means
xb <- data.frame(cl = factor(c(1, 2)), x1 = c(10,11), x2 = c(11, 9))
```

<br><br><br>

Select $k=2$, and set initial seed means <br>
$\bar{x}_1=$ (`r xb[1,-1]`) ,
$\bar{x}_2=$ (`r xb[2,-1]`) <br>

:::
::: {.column}

```{r out.width="100%", fig.width=4, fig.height=4}
#| echo: false
mn <- data.frame(cl=xb$cl, lbl=c("m1", "m2"), 
                 x1=xb$x1, x2=xb$x2)
ggplot(data=df, aes(x1, x2)) + geom_text(aes(label=lbl)) + 
  xlab("") + ylab("") + theme_bw() + 
  xlim(c(1,22)) + ylim(c(1,22)) +
  geom_point(data=mn, aes(x=x1, y=x2, color=cl), 
             shape=3, size=3) + 
  scale_colour_brewer("", palette="Dark2") +
  theme(aspect.ratio=1, legend.position="none") 
```
:::
::::

## k-means clustering - algorithm [(4/8)]{.smallest}

:::: {.columns}
::: {.column}



<br><br><br>


Compute distances $(d_1, d_2)$ between each observation and each mean, 

$\bar{x}_1=$ (`r xb[1,-1]`) ,
$\bar{x}_2=$ (`r xb[2,-1]`) <br>


:::
::: {.column}

```{r}
#| echo: false
# Compute distances between each observation and each mean

d1 <- apply(df[,2:3], 1, function(x) round(sqrt(sum((x-xb[1,2:3])^2)),1))
d2 <- apply(df[,2:3], 1, function(x) round(sqrt(sum((x-xb[2,2:3])^2)),1))
df.km <- cbind(df, d1, d2)
df.km %>%
  kable("html", escape=F) %>%
  kable_styling("striped", position = "center", 
                font_size=24) %>% 
  row_spec(0, color = "white", background = "#3F9F7A") %>%
  column_spec(1:5, border_right=TRUE, width="1cm") %>%
  column_spec(4:5, color="#7570B3") 
```

:::

::::

## k-means clustering - algorithm [(5/8)]{.smallest}

:::: {.columns}
::: {.column}

<br><br><br><br>
Assign each observation to a cluster, based on which mean is closest.

:::
::: {.column}



```{r}
#| echo: false
df.km$cl <- ifelse(d1 < d2, 1, 2)
df.km$cl <- factor(df.km$cl)
kable(df.km) %>%
  kable_styling("striped", position = "center", 
                font_size=24) %>% 
  row_spec(0, color = "white", background = "#3F9F7A") %>%
  column_spec(1:6, border_right=TRUE, width="1cm") %>%
  column_spec(4:6, color="#7570B3")

```
:::
::::

## k-means clustering - algorithm [(6/8)]{.smallest}

:::: {.columns}
::: {.column}
Recompute means, and re-assign the cluster membership

```{r}
#| echo: false
xb <- df.km %>%
  group_by(cl) %>%
  summarise(x1=mean(x1), x2=mean(x2))
xb1 <- data.frame(x1=xb$x1[1], x2=xb$x2[1])
xb2 <- data.frame(x1=xb$x1[2], x2=xb$x2[2])
```

$\bar{x}_1=$ (`r round(xb[1,-1], 0)`) ,
$\bar{x}_2=$ (`r round(xb[2,-1], 0)`) <br>

```{r}
#| echo: false
d1 <- apply(df[,2:3], 1, function(x) round(sqrt(sum((x-xb[1,2:3])^2)),1))
d2 <- apply(df[,2:3], 1, function(x) round(sqrt(sum((x-xb[2,2:3])^2)),1))
df.km$d1 <- round(d1, 1)
df.km$d2 <- round(d2, 1)

df.km$cl <- ifelse(d1 < d2, 1, 2)
df.km$cl <- factor(df.km$cl)
kable(df.km) %>%
  kable_styling("striped", position = "center", 
                font_size=24) %>% 
  row_spec(0, color = "white", background = "#3F9F7A") %>%
  column_spec(1:6, border_right=TRUE, width="1cm") %>%
  column_spec(4:6, color="#7570B3")

```

:::
::: {.column}

```{r out.width="100%", fig.width=4, fig.height=4}
#| echo: false
mn <- data.frame(cl=xb$cl, lbl=c("m1", "m2"), 
                 x1=xb$x1, x2=xb$x2)
ggplot() + 
  geom_text(data=df.km, aes(x=x1, y=x2, label=lbl, color=cl)) + 
  geom_point(data=mn, aes(x=x1, y=x2, color=cl), 
             shape=3, size=3) + 
  scale_colour_brewer("", palette="Dark2") +
  xlab("") + ylab("") + theme_bw() +
  xlim(c(1,22)) + ylim(c(1,22)) +
  theme(aspect.ratio=1, legend.position="None") 
```

:::
::::

## k-means clustering - algorithm [(7/8)]{.smallest}

:::: {.columns}
::: {.column}
Recompute means, and re-assign the cluster membership

```{r}
#| echo: false
xb <- df.km %>%
  group_by(cl) %>%
  summarise(x1=mean(x1), x2=mean(x2))
xb1 <- data.frame(x1=xb$x1[1], x2=xb$x2[1])
xb2 <- data.frame(x1=xb$x1[2], x2=xb$x2[2])
```

$\bar{x}_1=$ (`r round(xb[1,-1], 0)`) ,
$\bar{x}_2=$ (`r round(xb[2,-1], 0)`) <br>

```{r}
#| echo: false
d1 <- apply(df[,2:3], 1, function(x) round(sqrt(sum((x-xb[1,2:3])^2)),1))
d2 <- apply(df[,2:3], 1, function(x) round(sqrt(sum((x-xb[2,2:3])^2)),1))
df.km$d1 <- round(d1, 1)
df.km$d2 <- round(d2, 1)

df.km$cl <- ifelse(d1 < d2, 1, 2)
df.km$cl <- factor(df.km$cl)
kable(df.km) %>%
  kable_styling("striped", position = "center", 
                font_size=24) %>% 
  row_spec(0, color = "white", background = "#3F9F7A") %>%
  column_spec(1:6, border_right=TRUE, width="1cm") %>%
  column_spec(4:6, color="#7570B3")

```
:::
::: {.column}

```{r out.width="100%", fig.width=4, fig.height=4}
#| echo: false
mn <- data.frame(cl=xb$cl, lbl=c("m1", "m2"), 
                 x1=xb$x1, x2=xb$x2)
ggplot() + 
  geom_text(data=df.km, aes(x=x1, y=x2, label=lbl, color=cl)) + 
  geom_point(data=mn, aes(x=x1, y=x2, color=cl), 
             shape=3, size=3) + 
  scale_colour_brewer("", palette="Dark2") +
  xlab("") + ylab("") + theme_bw() +
  xlim(c(1,22)) + ylim(c(1,22)) +
  theme(aspect.ratio=1, legend.position="None") 
```
:::
::::

## k-means clustering - algorithm [(8/8)]{.smallest}

:::: {.columns}
::: {.column}
Recompute means, and re-assign the cluster membership

```{r}
#| echo: false
xb <- df.km %>%
  group_by(cl) %>%
  summarise(x1=mean(x1), x2=mean(x2))
xb1 <- data.frame(x1=xb$x1[1], x2=xb$x2[1])
xb2 <- data.frame(x1=xb$x1[2], x2=xb$x2[2])
```

$\bar{x}_1=$ (`r round(xb[1,-1], 0)`) ,
$\bar{x}_2=$ (`r round(xb[2,-1], 0)`) <br>

```{r}
#| echo: false
d1 <- apply(df[,2:3], 1, function(x) round(sqrt(sum((x-xb[1,2:3])^2)),1))
d2 <- apply(df[,2:3], 1, function(x) round(sqrt(sum((x-xb[2,2:3])^2)),1))
df.km$d1 <- round(d1, 1)
df.km$d2 <- round(d2, 1)

df.km$cl <- ifelse(d1 < d2, 1, 2)
df.km$cl <- factor(df.km$cl)
kable(df.km) %>%
  kable_styling("striped", position = "center", 
                font_size=24) %>% 
  row_spec(0, color = "white", background = "#3F9F7A") %>%
  column_spec(1:6, border_right=TRUE, width="1cm") %>%
  column_spec(4:6, color="#7570B3")

```
:::
::: {.column}

```{r out.width="100%", fig.width=4, fig.height=4}
#| echo: false
mn <- data.frame(cl=xb$cl, lbl=c("m1", "m2"), 
                 x1=xb$x1, x2=xb$x2)
ggplot() + 
  geom_text(data=df.km, aes(x=x1, y=x2, label=lbl, color=cl)) + 
  geom_point(data=mn, aes(x=x1, y=x2, color=cl), 
             shape=3, size=3) + 
  scale_colour_brewer("", palette="Dark2") +
  xlab("") + ylab("") + theme_bw() +
  xlim(c(1,22)) + ylim(c(1,22)) +
  theme(aspect.ratio=1, legend.position="None") 
```
:::
::::

## Example: penguins

:::: {.columns}
::: {.column}
- We know there are three clusters, but generally [we don't know this]{.monash-orange2}. 
- Will $k=3$-means clustering see three?
- Fit for various values of $k$. Add cluster label to data.
- Examine [solution in plots]{.monash-blue2} of the data.
- Compute [cluster metrics]{.monash-blue2}.

- NOTE: `set.seed()` because results can depend on initialisation.
:::
::: {.column style="font-size: 70%;"}

```{r}
set.seed(712)
p_km3 <- kmeans(p_std[,2:5], 3)
p_std_km <- p_std |>
  mutate(cl = factor(p_km3$cluster))
```

```{r}
#| echo: false
#| eval: false
animate_xy(p_std_km[,2:5], col=p_std_km$cl)
render_gif(p_std_km[,2:5],
           grand_tour(),
           display_xy(col=p_std_km$cl),
           width=400,
           height=400,
           frames=100,
           gif_file = "../gifs/p_km.gif")
```

<center>
![](../gifs/p_km.gif){width=500}
</center>
:::
::::

## Choosing $k$ with cluster statistics [(1/2)]{.smallest}



- [within.cluster.ss]{.monash-blue2}: sum of distances within cluster. Want it to be [low]{.monash-orange2}, but always drops for each additional cluster so look for large drops.
- [WBRatio]{.monash-blue2}: average within/average between distances. Want it to be [low]{.monash-orange2}, but always drops for each additional cluster so look for large drops.
- [Hubert Gamma]{.monash-blue2}: (s+ - s-)/(s+ + s-) where $s+=$sum of number of within $<$ between, $s-=$ sum of number within $>$ between. Want this to be [high]{.monash-orange2}.
- [Dunn]{.monash-blue2}: ratio of (smallest distance between points from different clusters) to (maximum distance of points within any cluster). Want this to be [high]{.monash-orange2}.
- [Calinski-Harabasz Index]{.monash-blue2}: $\frac{\sum_{i=1}^p B_{ii}/(k-1)}{\sum_{i=1}^p W_{ii}/(n-k)}$. Want this to be [high]{.monash-orange2}.


## Choosing $k$ with cluster statistics [(2/2)]{.smallest}

:::: {.columns}
::: {.column width=70%}

```{r}
#| echo: false
#| fig-width: 8
#| fig-height: 6
#| out-width: 90%
set.seed(31)
f.km <- NULL; f.km.stats <- NULL
for (i in 2:10) {
  cl <- kmeans(p_std[,2:5], i, nstart=5)$cluster
  x <- cluster.stats(dist(p_std[,2:5]), cl)
  f.km <- cbind(f.km, cl)
  f.km.stats <- rbind(f.km.stats, 
                      c(x$within.cluster.ss, 
                        x$wb.ratio, x$ch,
                        x$pearsongamma, x$dunn, 
                        x$dunn2))
}
colnames(f.km.stats) <- c("within.cluster.ss",
                          "wb.ratio", "ch",
                          "pearsongamma", "dunn",
                          "dunn2")
f.km.stats <- data.frame(f.km.stats)
f.km.stats$cl <- 2:10
f.km.stats.m <- f.km.stats %>% 
  pivot_longer(`within.cluster.ss`:dunn2,
               names_to="stat", 
               values_to="value") |>
  mutate(stat = factor(stat,
                       levels=colnames(f.km.stats)))
ggplot(data=f.km.stats.m) + 
  geom_line(aes(x=cl, y=value)) + xlab("# clusters") + ylab("") +
  facet_wrap(~stat, ncol=3, scales = "free_y") + 
  theme_bw()
```
:::
::: {.column width=30%}
- Results are inconclusive. No agreement between metrics. 
- Not unusual. Stay tuned for nuisance variables and observations. 
:::

::::

## Hierarchical clustering {.transition-slide .center}

## Hierarchical clustering [1/4]{.smallest}

<br>

- [Agglomeration]{.monash-orange2}: Begin with all observations in singleton clusters. Sequentially [join]{.monash-orange2} points into clusters, until all are in one cluster. 
- [Divisive]{.monash-orange2}: Begin with all observtions in one cluster, adn sequentially [divide]{.monash-orange2} until all observations are in singleton clusters. 
- Produces a tree diagram illustrating the process, called a [dendrogram]{.monash-orange2}.

## Hierarchical clustering [2/4]{.smallest}

:::: {.columns}
::: {.column style="font-size: 80%;"}

$n\times n$ distance matrix `r set.seed(1)`

```{r}
#| echo: false
d <- as.matrix(round(dist(df[,2:3], diag=TRUE, upper=TRUE),1))
colnames(d) <- df$lbl
rownames(d) <- df$lbl
kable(d) %>%
  kable_styling("striped", position = "center", 
                row_label_position = "c", 
                font_size=18) %>%
  row_spec(0, color = "white", background = "#3F9F7A") %>%
  column_spec(1, color = "white", background = "#3F9F7A") %>%
  column_spec(1:12, border_right=TRUE, width="1cm") 
```

:::
::: {.column}

```{r out.width="60%", fig.width=4, fig.height=4}
#| echo: false
ggplot(data=df, aes(x1, x2)) + geom_text(aes(label=lbl)) + 
  xlab("") + ylab("") + theme_bw() + 
  theme(aspect.ratio=1) 
```

:::
::::

## Hierarchical clustering [3/4]{.smallest}

:::: {.columns}
::: {.column style="font-size: 80%;"}

```{r out.width="100%", fig.width=4, fig.height=4}
#| echo: false
ggplot(data=df, aes(x1, x2)) + geom_text(aes(label=lbl)) + 
  xlab("") + ylab("") + theme_bw() + 
  theme(aspect.ratio=1) 
```

:::
::: {.column}

```{r out.width="100%", fig.width=4, fig.height=4}
#| echo: false
df_hc <- hclust(as.dist(d), method="average")
ggdendrogram(df_hc, rotate = TRUE, size = 4)
```

:::
::::

## Hierarchical clustering [4/4]{.smallest}

:::: {.columns}
::: {.column}
```{r out.width="100%", fig.width=4, fig.height=4}
#| echo: false
df$cl11 <- factor(c(1,2,1,2,1,1,1,1,1,1,1,1))
ggplot(data=df, aes(x=x1, y=x2, colour=cl11)) +
  geom_text(aes(label=lbl)) + 
  scale_color_brewer("", palette="Dark2") +
  xlab("") + ylab("") + theme_bw() + 
  theme(aspect.ratio=1, legend.position="none") 
```

:::
::: {.column}
```{r out.width="100%", fig.width=4, fig.height=4}
#| echo: false
ggdendrogram(df_hc, rotate = TRUE, size = 4)
```
:::
::::


## Linkage

:::: {.columns}
::: {.column}
What is the [distance]{.monash-orange2} between the new [cluster (d,b)]{.monash-orange2} and all of the other observations?
:::
::: {.column}
Between points [in]{.monash-orange2} the cluster to points [not in]{.monash-orange2} the cluster. 

- [Single]{.monash-orange2}: minimum distance between points in the different clusters
- [Complete]{.monash-orange2}: maximum distance between points in the different clusters
- [Average]{.monash-orange2}: average of distances between points in the different clusters
- [Centroid]{.monash-orange2}: distances between the average of the different clusters
- [Wards]{.monash-orange2}: minimizes the total within-cluster variance
:::
::::

## Linkage

![](../images/linkage.png)


## Calculations with different linkage choices

:::: {.columns}
::: {.column style="font-size: 80%;"}

```{r}
#| echo: false
kable(d) %>%
  kable_styling("striped", position = "center", 
                row_label_position = "c", 
                font_size=20) %>%
  row_spec(0, color = "white", background = "#3F9F7A") %>%
  column_spec(1, color = "white", background = "#3F9F7A") %>%
  column_spec(1:12, border_right=TRUE, width="1cm") %>%
  column_spec(c(3,5), color="#D95F02") %>%
  row_spec(c(2,4), color="#D95F02") %>%
  column_spec(c(2,4), color="#7570B3") %>%
  row_spec(c(1,3), color="#7570B3")
``` 

Distance (b,d): <br>
Distance (a,c):

:::
::: {.column style="font-size: 80%;"}

```{r out.width="50%", fig.width=3, fig.height=3}
#| echo: false
df$cl11 <- factor(c(3,2,3,2,1,1,1,1,1,1,1,1))
ggplot(data=df, aes(x=x1, y=x2, colour=cl11)) +
  geom_text(aes(label=lbl)) + 
  scale_color_brewer("", palette="Dark2") +
  xlab("") + ylab("") + theme_bw() + 
  theme(aspect.ratio=1, legend.position="none") 
```

Linkage between (b,d) and (a,c)

Single:  <br>
Complete:  <br>
Average:  <br>
Centroid:  <br>


:::
::::

## Results from different linkage choices

```{r out.width="80%", fig.width=10, fig.height=6}
#| echo: false
df_hc1 <- hclust(as.dist(d), method="single")
p1 <- ggdendrogram(df_hc1, rotate = TRUE, size = 4) + ggtitle("single")
df_hc2 <- hclust(as.dist(d), method="complete")
p2 <- ggdendrogram(df_hc2, rotate = TRUE, size = 4) + ggtitle("complete")
df_hc3 <- hclust(as.dist(d), method="average")
p3 <- ggdendrogram(df_hc3, rotate = TRUE, size = 4) + ggtitle("average")
df_hc4 <- hclust(as.dist(d), method="centroid")
p4 <- ggdendrogram(df_hc4, rotate = TRUE, size = 4) + ggtitle("centroid")
df_hc5 <- hclust(as.dist(d), method="ward.D2")
p5 <- ggdendrogram(df_hc5, rotate = TRUE, size = 4) + ggtitle("wards")
p6 <- ggplot(data=df, aes(x=x1, y=x2)) +
  geom_text(aes(label=lbl)) + 
  xlab("") + ylab("") + theme_bw() + 
  theme(aspect.ratio=1, legend.position="none") 
grid.arrange(p1, p2, p3, p4, p5, p6, ncol=3)
```

## Dendrogram

:::: {.columns}
::: {.column}

- Each [leaf]{.monash-orange2} of the dendrogram represents one observation
- Leaves [fuse]{.monash-orange2} into branches and branches fuse, either with leaves or other branches.	
- Fusions [lower in the tree]{.monash-orange2} mean the groups of observations are more similar to each other. 

:::
::: {.column}
[Cut the tree]{.monash-orange2} to partition the data into $k$ clusters.

:::
::::

## Results from different linkage choices

*Model-in-the-data-space*

```{r out.width="90%", fig.width=10, fig.height=6}
#| echo: false
df$cl1 <- factor(cutree(df_hc1, 3))
df$cl2 <- factor(cutree(df_hc2, 3))
df$cl3 <- factor(cutree(df_hc3, 3))
df$cl4 <- factor(cutree(df_hc4, 3))
df$cl5 <- factor(cutree(df_hc5, 3))
p1 <- ggplot(data=df, aes(x=x1, y=x2, colour=cl1)) +
  geom_text(aes(label=lbl)) + 
  scale_color_brewer("", palette="Dark2") +
  xlab("") + ylab("") + theme_bw() + 
  theme(aspect.ratio=1, legend.position="none") + ggtitle("single")
p2 <- ggplot(data=df, aes(x=x1, y=x2, colour=cl2)) +
  geom_text(aes(label=lbl)) + 
  scale_color_brewer("", palette="Dark2") +
  xlab("") + ylab("") + theme_bw() + 
  theme(aspect.ratio=1, legend.position="none") + 
  ggtitle("complete")
p3 <- ggplot(data=df, aes(x=x1, y=x2, colour=cl3)) +
  geom_text(aes(label=lbl)) + 
  scale_color_brewer("", palette="Dark2") +
  xlab("") + ylab("") + theme_bw() + 
  theme(aspect.ratio=1, legend.position="none") + 
  ggtitle("average")
p4 <- ggplot(data=df, aes(x=x1, y=x2, colour=cl4)) +
  geom_text(aes(label=lbl)) + 
  scale_color_brewer("", palette="Dark2") +
  xlab("") + ylab("") + theme_bw() + 
  theme(aspect.ratio=1, legend.position="none") + 
  ggtitle("centroid")
p5 <- ggplot(data=df, aes(x=x1, y=x2, colour=cl5)) +
  geom_text(aes(label=lbl)) + 
  scale_color_brewer("", palette="Dark2") +
  xlab("") + ylab("") + theme_bw() + 
  theme(aspect.ratio=1, legend.position="none") + 
  ggtitle("wards")
grid.arrange(p1, p2, p3, p4, p5, ncol=3)
```

## Example: penguins

:::: {.columns}
::: {.column}
- We know there are three clusters, but generally [we don't know this]{.monash-orange2}. 
- Will $k=3$-means clustering see three?
- Fit for various values of $k$. Add cluster label to data.
- Examine [solution in plots]{.monash-blue2} of the data.
- Compute [cluster metrics]{.monash-blue2}.

- NOTE: No need for `set.seed()` because results are deterministic.
:::
::: {.column style="font-size: 70%;"}

```{r}
p_hc_w3 <- hclust(dist(p_std[,2:5]), method="ward.D2")
p_std_hc_w3 <- p_std |>
  mutate(cl = factor(cutree(p_hc_w3, 3)))
```

```{r}
#| echo: false
#| eval: false
animate_xy(p_std_hc_w3[,2:5], col=p_std_hc_w3$cl)
render_gif(p_std_hc_w3[,2:5],
           grand_tour(),
           display_xy(col=p_std_hc_w3$cl),
           width=400,
           height=400,
           frames=100,
           gif_file = "../gifs/p_hc_w3.gif")
```

<center>
![](../gifs/p_hc_w3.gif){width=500}
</center>
:::
::::

## Choosing $k$ with cluster statistics 

:::: {.columns}
::: {.column width=70%}

```{r}
#| echo: false
#| fig-width: 8
#| fig-height: 6
#| out-width: 90%
f.hc <- NULL; f.hc.stats <- NULL
for (i in 2:10) {
  cl <- cutree(p_hc_w3, i)
  x <- cluster.stats(dist(p_std[,2:5]), cl)
  f.hc <- cbind(f.hc, cl)
  f.hc.stats <- rbind(f.hc.stats, 
                      c(x$within.cluster.ss, 
                        x$wb.ratio, x$ch,
                        x$pearsongamma, x$dunn, 
                        x$dunn2))
}
colnames(f.hc.stats) <- c("within.cluster.ss",
                          "wb.ratio", "ch",
                          "pearsongamma", "dunn",
                          "dunn2")
f.hc.stats <- data.frame(f.hc.stats)
f.hc.stats$cl <- 2:10
f.hc.stats.m <- f.hc.stats %>% 
  pivot_longer(`within.cluster.ss`:dunn2,
               names_to="stat", 
               values_to="value") |>
  mutate(stat = factor(stat,
                       levels=colnames(f.hc.stats)))
ggplot(data=f.hc.stats.m) + 
  geom_line(aes(x=cl, y=value)) + xlab("# clusters") + ylab("") +
  facet_wrap(~stat, ncol=3, scales = "free_y") + 
  theme_bw()
```
:::
::: {.column width=30%}
- `within.cluster.ss` and `wb.ratio` suggest 3, and 5
- `pearsongamma` (Hubert) suggests 2-3
- `dunn`, `dunn2`, `ch` all 2?
:::

::::

<!--## There is no truth

You don't know what you don't know!

- Visualisation is really important to help decide in a solution is reasonable
- External validation
-->

## Next: Model-based clustering and self-organising maps {.transition-slide .center}


