knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) options(datatable.print.topn = 3L)
We'll start by loading the necessary libraries:
library(data.table) library(ggplot2) library(gmeans)
To illustrate the purpose of the G-means algorithm, let's start with an adapted
k-means clustering example from the
tidymodels 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:
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.
k-means clustering requires specifying the number of clusters, k, beforehand.
To illustrate this, let's fit a k-means model with k = 3:
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:
tidy <- function(x, col.names = colnames(x$centers)) { if (is.null(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:
augment(kclust, points)
The tidy() function provides a per-cluster summary, displaying the cluster centers,
sizes, and within-cluster sum of squares:
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:
glance(kclust)
Using these helper functions, we can easily extract and manipulate the results of
k-means clustering for different values of k:
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
To enhance the visualization, let's add cluster centers:
p2 <- p1 + geom_point(data = clusters, size = 10, shape = "x") + labs(title = "k-means Clustering with Centers") p2
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:
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.
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.
Let's now apply the G-means algorithm to the same data:
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:
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:
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:
tidy(gclust)
The glance() function gives an overall summary of the model:
glance(gclust)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.