---
title: "Why G-means?"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Why G-means?}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
options(datatable.print.topn = 3L)
```

We'll start by loading the necessary libraries:

```{r setup, message = FALSE}
library(data.table)
library(ggplot2)
library(gmeans)
```

## Overview

To illustrate the purpose of the G-means algorithm, let's start with an adapted
k-means clustering example from the
[tidymodels](https://www.tidymodels.org/learn/statistics/k-means/) website.
This example shows the challenge of determining the number of clusters (`k`) in
clustering analysis. Throughout this vignette, we will use `data.table` for data
manipulation, and the custom `tidy`, `augment`, and `glance` functions for handling
model output, inspired by the `broom` package functionality.

We begin by generating some random two-dimensional data that naturally forms three
clusters. Each cluster's data comes from a different multivariate Gaussian
distribution with unique means:

```{r, dpi = 300}
set.seed(27)

centers <- data.table(
  cluster = factor(1:3),
  num_points = c(100, 150, 50),
  x1 = c(5, 0, -3),
  x2 = c(-1, 1, -2)
)
points <- centers[, .(
  x1 = rnorm(num_points, mean = x1),
  x2 = rnorm(num_points, mean = x2)
), by = cluster]

ggplot(points, aes(x1, x2, color = cluster)) +
  geom_point(alpha = 0.3)
```

In this simple example, we know that there are three clusters.
However, in real-world scenarios, the number of clusters is often unknown and must
be determined as part of the analysis.

## Challenges with k-means

k-means clustering requires specifying the number of clusters, `k`, beforehand.
To illustrate this, let's fit a k-means model with `k = 3`:

```{r}
points <- points[, cluster := NULL]
kclust <- kmeans(points, centers = 3)
kclust
```

Here, we fit the k-means model with the correct number of clusters because we know the
true structure of the data. However, this knowledge is often not available in practice.

To explore the effect of different `k`, we can fit k-means models with varying numbers
of clusters and visualize the results. To make handling the k-means output easier, we
define `tidy`, `augment`, and `glance` functions that mimic the functionality of the
`broom` package:

```{r}
`%||%` <- function(x, y) if (!is.null(x)) x else y

tidy <- function(x, col.names = colnames(x$centers)) {
  col.names <- col.names %||% paste0("x", seq_len(ncol(x$centers)))
  dt <- as.data.table(x$centers)
  setnames(dt, col.names)
  dt[, let(
    size = x$size,
    withinss = x$withinss,
    cluster = factor(seq_len(.N))
  )][]
}

augment <- function(x, data) {
  if (inherits(data, "matrix") && is.null(colnames(data))) {
    colnames(data) <- paste0("X", seq_len(ncol(data)))
  }
  dt <- as.data.table(data)
  dt[, .cluster := as.factor(x$cluster)][]
}

glance <- function(x) {
  as.data.table(x[c("totss", "tot.withinss", "betweenss", "iter")])
}
```

The `augment()` function adds the cluster assignments to the original dataset,
allowing us to see how each data point is classified:

```{r}
augment(kclust, points)
```

The `tidy()` function provides a per-cluster summary, displaying the cluster centers,
sizes, and within-cluster sum of squares:

```{r}
tidy(kclust)
```

To obtain a single-row summary with overall metrics such as total sum of squares and
the number of iterations, use the `glance()` function:

```{r}
glance(kclust)
```

Using these helper functions, we can easily extract and manipulate the results of
k-means clustering for different values of `k`:

```{r, dpi = 300}
kclusts <- data.table(k = 1:9)
kclusts[, kclust := lapply(k, \(x) kmeans(points, x))]
kclusts[, let(
  tidied = lapply(kclust, tidy),
  glanced = lapply(kclust, glance),
  augmented = lapply(kclust, augment, points)
)]

clusters <- kclusts[, .(k, rbindlist(tidied))]
assignments <- kclusts[, .(k, rbindlist(augmented))]
clusterings <- kclusts[, .(k, rbindlist(glanced))]

p1 <- ggplot(assignments, aes(x = x1, y = x2)) +
  geom_point(aes(color = .cluster), alpha = 0.8) +
  facet_wrap(~k) +
  labs(title = "k-means Clustering Results with Different Values of k")
p1
```

## Visualizing cluster centers

To enhance the visualization, let's add cluster centers:

```{r}
p2 <- p1 + geom_point(data = clusters, size = 10, shape = "x") +
  labs(title = "k-means Clustering with Centers")
p2
```

## Evaluating clustering performance

Finally, we can look at how the total within-cluster sum of squares (WSS) changes with
different values of `k`. This helps us see how well the data is being clustered as `k`
increases:

```{r, dpi = 300}
ggplot(clusterings, aes(k, tot.withinss)) +
  geom_line() +
  geom_point() +
  labs(
    title = "Total Within-Cluster Sum of Squares vs. Number of Clusters (k)",
    x = "Number of Clusters (k)",
    y = "Total Within-Cluster Sum of Squares"
  )
```

In general, the WSS decreases as the number of clusters k increases, which is expected
since having more clusters usually results in a better fit. However, we often look for
a point in the plot where the decrease in WSS starts to slow down, creating a
noticeable "elbow". This elbow suggests that adding more clusters beyond this
point offers little improvement, indicating a good number of clusters. In our example,
this bend is around `k = 3`, suggesting that three clusters capture the main structure
of the data effectively.

## Motivation for G-means

As seen from the plots, choosing the right number of clusters is not straightforward.
We could use metrics like WSS to help decide, but these methods can be subjective and
prone to error. This is where the G-means algorithm comes in: it automatically
determines the number of clusters by assessing the data distribution within each
cluster.

By using statistical hypothesis testing (the Anderson-Darling test in our
implementation), G-means provides a more robust and automated way to find the "correct"
number of clusters. In the next section, we'll see how to use G-means and explore its
benefits over traditional k-means clustering.

## G-means

Let's now apply the G-means algorithm to the same data:

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

gmeans(points)
```

As expected from our previous analysis, G-means identifies 3 clusters, aligning with
the elbow point observed in the WSS plot.

Next, let's explore how G-means performs on a different dataset:

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

x <- as.matrix(iris[, -5])
gclust <- gmeans(x)
```

Since `gmeans()` uses `stats::kmeans()` under the hood, we can use our previously
defined helper functions to analyze the clustering results.

The `augment()` function adds cluster assignments to the original dataset for easy
plotting:

```{r, dpi = 300}
augment(gclust, x) |>
  ggplot(aes(x = Petal.Length, y = Petal.Width)) +
  geom_point(aes(color = .cluster))
```

The `tidy()` function provides a summary of each cluster:

```{r}
tidy(gclust)
```

The `glance()` function gives an overall summary of the model:

```{r}
glance(gclust)
```