49 Clustering in rtemis

cluster() fits a clustering with a configuration from one of the setup_*() functions listed by available_clustering(). It returns a HardClustering, with one cluster label per case, or a SoftClustering, which adds a matrix of memberships. This chapter clusters the patients of the diabetes data by their nine numeric baseline measurements, chooses the number of clusters, compares methods, and describes the clusters by the outcome that the clustering did not use.

available_clustering()

           CMeans: Fuzzy c-means clustering
           DBSCAN: DBSCAN density-based clustering
              GMM: Gaussian mixture model clustering
           HOPACH: Hierarchical Ordered Partitioning and Collapsing Hybrid
           HardCL: Hard competitive learning
           KMeans: K-means clustering
        NeuralGas: Neural Gas clustering
              PAM: Partitioning Around Medoids
             PAMK: Partitioning Around Medoids with estimation of the number of clusters
  SpectralLaplace: Spectral clustering with an exponential kernel, which decays more slowly than the Gaussian and so keeps more weight on distant pairs
    SpectralLocal: Spectral clustering with a Gaussian kernel whose width is set per case from the distance to its seventh nearest neighbor, which lets one clustering hold groups of differing density
      SpectralRBF: Spectral clustering with a Gaussian kernel of one shared width

49.1 Data

The nine numeric features are standardized, so that each contributes equally to the Euclidean distance (Section 48.1):

db <- read_diabetes()
x <- as.data.table(scale(db[, !c("sex", "progression")]))
ex <- setup_SerialExecution(seed = 2026L)

49.2 k-means and the number of clusters

k-means with \(k\) from 2 to 6, with the average silhouette width of each partition computed by cluster::silhouette():

d <- dist(x)
asw <- rbindlist(lapply(2:6, function(k) {
  fit <- cluster(x, config = setup_KMeans(k = k), execution_config = ex, verbosity = 0L)
  sil <- cluster::silhouette(as.integer(fit@clusters), d)
  data.table(k = k, average_silhouette = mean(sil[, "sil_width"]))
}))
asw

       k average_silhouette
   <int>              <num>
1:     2             0.2535
2:     3             0.1712
3:     4             0.1504
4:     5             0.1393
5:     6             0.1399

draw_line(
  asw[["k"]], list(`Average silhouette width` = asw[["average_silhouette"]]),
  xlab = "Number of clusters", ylab = "Average silhouette width"
)

The average silhouette width is highest at \(k\) = 2, at 0.25. Kaufman and Rousseeuw (1990) describe an average silhouette width of 0.25 or less as indicating no substantial structure, and 0.26 to 0.50 as weak structure that could be artificial. The values here are at or below that boundary for every \(k\): the patients form a continuum rather than separated groups, as the embeddings of Section 47.3 suggested.

setup_PAMK() fits partitioning around medoids for a range of \(k\) and selects \(k\) by a criterion, here the average silhouette width:

fit_pamk <- cluster(
  x,
  config = setup_PAMK(krange = 2:6, criterion = setup_ASWCriterion()),
  execution_config = ex,
  verbosity = 0L
)
table(fit_pamk@clusters)


  1   2 
227 215

49.3 A Gaussian mixture

setup_GMM() fits Gaussian mixture models with the mclust package and, with k = NULL, selects the number of components and the covariance model by BIC:

fit_gmm <- cluster(x, config = setup_GMM(), execution_config = ex, verbosity = 0L)
fit_gmm@k

[1] 6

table(fit_gmm@clusters)


  1   2   3   4   5   6 
 46 132  21  78  79  86

The soft clustering gives each patient a posterior probability for each component. The median of each patient’s largest posterior probability is 1.00, and 1.8% of patients have no component with a posterior probability of 0.8 or more.

BIC and the silhouette select different numbers of clusters. BIC rewards a better description of the density, which several overlapping normal components can provide for a single continuous cloud of cases; the silhouette rewards separated groups.

49.4 Agreement between methods

The adjusted Rand index compares two partitions of the same cases (Section 48.7):

fit_km3 <- cluster(x, config = setup_KMeans(k = 3L), execution_config = ex, verbosity = 0L)
fit_pam3 <- cluster(x, config = setup_PAM(k = 3L), execution_config = ex, verbosity = 0L)
fit_sp3 <- cluster(x, config = setup_SpectralRBF(k = 3L), execution_config = ex, verbosity = 0L)
c(
  kmeans_vs_pam = mclust::adjustedRandIndex(fit_km3@clusters, fit_pam3@clusters),
  kmeans_vs_spectral = mclust::adjustedRandIndex(fit_km3@clusters, fit_sp3@clusters),
  kmeans_vs_gmm = mclust::adjustedRandIndex(fit_km3@clusters, fit_gmm@clusters)
)

     kmeans_vs_pam kmeans_vs_spectral      kmeans_vs_gmm 
            0.3813             0.5679             0.1647

Methods with the same number of clusters agree partly. In data without separated groups, each method divides the continuum according to its own objective, and the partitions differ accordingly.

49.5 Density-based clustering

setup_DBSCAN() requires the neighborhood radius eps. The distance of each case to its 5th nearest neighbor shows the scale of local distances:

knn_dist <- sort(dbscan::kNNdist(x, k = 5L))
draw_line(
  seq_along(knn_dist), list(`5th nearest neighbor distance` = knn_dist),
  points = FALSE, xlab = "Case (sorted)", ylab = "Distance"
)

dbscan_fits <- rbindlist(lapply(c(0.5, 1.5, 2.5), function(eps) {
  fit <- cluster(x, config = setup_DBSCAN(eps = eps, min_points = 5L), execution_config = ex, verbosity = 0L)
  data.table(
    eps = eps,
    n_clusters = fit@k,
    noise = sum(fit@clusters == 0L)
  )
}))
dbscan_fits

     eps n_clusters noise
   <num>      <int> <int>
1:   0.5          0   442
2:   1.5          1   199
3:   2.5          1     4

With eps = 0.5, smaller than almost every nearest-neighbor distance in nine standardized dimensions, every patient is noise. With eps = 1.5, DBSCAN finds a single cluster of 243 patients and labels the rest noise; with eps = 2.5, a single cluster containing all but 4. At none of the three radii does it find more than one dense region.

49.6 Describing clusters by an external variable

The k-means clusters with \(k = 3\) can be described by variables not used to form them, such as disease progression:

profile <- data.table(
  cluster = fit_km3@clusters,
  progression = db[["progression"]],
  bmi = db[["bmi"]],
  s5 = db[["s5"]]
)[, .(n = .N, progression = mean(progression), bmi = mean(bmi), s5 = mean(s5)), by = cluster]
setorder(profile, cluster)
profile

   cluster     n progression   bmi    s5
     <int> <int>       <num> <num> <num>
1:       1   137       178.7 28.48 4.761
2:       2   133       184.7 28.10 5.054
3:       3   172       105.8 23.36 4.227

pca_x <- prcomp(x)
draw_scatter(
  pca_x[["x"]][, 1], pca_x[["x"]][, 2],
  group = factor(paste("Cluster", fit_km3@clusters)),
  xlab = "PC1", ylab = "PC2",
  point_alpha = 0.6
)

The clusters differ in mean progression, since they differ in BMI and s5, which predict it. In the space of the first two principal components, the three clusters are adjacent regions of one cloud. They are a description of how these patients’ measurements vary, under the chosen features and distance, not evidence of three distinct types of patient (Section 48.7).