Merge pull request #1516 from exploratory-io/fix/issue-36106

kei51e · web-flow · commit 3fb1477baae6 · 2026-06-02T18:47:26.000-07:00
Add silhouette analysis to K-means clustering
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: exploratory
 Type: Package
 Title: R package for Exploratory
-Version: 15.1.12
+Version: 15.1.13
 Date: 2026-05-30
 Authors@R: c(person("Hideaki", "Hayashi", email = "hideaki@exploratory.io", role = c("aut", "cre")), person("Hide", "Kojima", email = "hide@exploratory.io", role = c("aut")), person("Kan", "Nishida", email = "kan@exploratory.io", role = c("aut")), person("Kei", "Saito", email = "kei@exploratory.io", role = c("aut")), person("Yosuke", "Yasuda", email = "double.y.919.quick@gmail.com", role = c("aut")))
 URL: https://github.com/exploratory-io/exploratory_func
diff --git a/R/kmeans.R b/R/kmeans.R
@@ -40,7 +40,8 @@ iterate_silhouette <- function(df, max_centers = 10,
                                algorithm = "Hartigan-Wong",
                                trace = FALSE,
                                normalize_data = TRUE,
-                               seed = NULL) {
+                               seed = NULL,
+                               dist = NULL) {
   mat <- as_numeric_matrix_(df, columns = colnames(df))
   if (normalize_data) {
     mat <- scale(mat)
@@ -54,7 +55,7 @@ iterate_silhouette <- function(df, max_centers = 10,
     return(data.frame(center = integer(0), avg_silhouette = numeric(0),
                       min_silhouette = numeric(0), pct_negative = numeric(0)))
   }
-  d <- stats::dist(mat)
+  d <- if (is.null(dist)) stats::dist(mat) else dist
   purrr::map_dfr(seq(2, upper), function(k) {
     if (!is.null(seed)) {
       set.seed(seed)
@@ -78,6 +79,45 @@ iterate_silhouette <- function(df, max_centers = 10,
   })
 }
 
+# Compute per-row silhouette widths for an already-built clustering.
+#
+# cluster_ids: integer cluster assignment vector (length n, aligned to mat rows).
+#              For K-Means this is x$kmeans$cluster (integer 1..k).
+# mat:         normalized numeric matrix used for clustering (n rows).
+# d:           optional precomputed stats::dist(mat) to reuse. Computed when NULL.
+#
+# Returns a tibble with n rows and columns silhouette_score, nearest_cluster, cluster_width.
+# Returns all-NA columns (no error) when silhouette is undefined
+# (fewer than 2 clusters, or fewer than 2 distinct points).
+compute_silhouette_per_row <- function(cluster_ids, mat, d = NULL) {
+  n <- length(cluster_ids)
+  na_result <- tibble::tibble(
+    silhouette_score = rep(NA_real_, n),
+    nearest_cluster = rep(NA_integer_, n),
+    cluster_width = rep(NA_real_, n)
+  )
+  ids <- as.integer(cluster_ids)
+  if (length(unique(ids)) < 2 || nrow(unique(mat)) < 2) {
+    return(na_result)
+  }
+  if (is.null(d)) {
+    d <- stats::dist(mat)
+  }
+  sil <- cluster::silhouette(ids, d)
+  if (!inherits(sil, "silhouette")) {
+    # silhouette() can return NA (not a matrix) for degenerate input.
+    return(na_result)
+  }
+  widths <- as.numeric(sil[, "sil_width"])
+  clusters <- as.integer(sil[, "cluster"])
+  clus_avg <- tapply(widths, clusters, mean, na.rm = TRUE)
+  tibble::tibble(
+    silhouette_score = widths,
+    nearest_cluster = as.integer(sil[, "neighbor"]),
+    cluster_width = as.numeric(clus_avg[as.character(clusters)])
+  )
+}
+
 #' analytics function for K-means view
 #' @export
 exp_kmeans <- function(df, ...,
@@ -154,6 +194,23 @@ exp_kmeans <- function(df, ...,
                    !!!rlang::syms(selected_cols))
   ret <- dplyr::ungroup(ret) # ungroup once so that the following mutate with purrr::map2 works.
 
+  # Build the normalized clustering matrix ONCE; the single dist() is the only
+  # O(n^2) cost and is shared with iterate_silhouette() (decision (a), #36106).
+  # Only build the shared dist for the ungrouped (UI) case; grouped models build
+  # their own per-group matrix below.
+  is_single_model <- length(ret$model) == 1
+  shared_sil_mat <- NULL
+  shared_sil_dist <- NULL
+  if (is_single_model) {
+    kmeans_df_shared <- df %>% dplyr::select(!!!rlang::syms(selected_cols))
+    shared_sil_mat <- as_numeric_matrix_(kmeans_df_shared, columns = colnames(kmeans_df_shared))
+    if (normalize_data) {
+      shared_sil_mat <- scale(shared_sil_mat)
+      shared_sil_mat[is.nan(shared_sil_mat)] <- 0
+    }
+    shared_sil_dist <- stats::dist(shared_sil_mat)
+  }
+
   # Compute elbow or silhouette results depending on optimal_method.
   elbow_result <- NULL
   silhouette_result <- NULL
@@ -176,13 +233,26 @@ exp_kmeans <- function(df, ...,
                                             algorithm = algorithm,
                                             trace = trace,
                                             normalize_data = normalize_data,
-                                            seed = NULL) # Seed is already set once at the top of exp_kmeans(). Skip it.
+                                            seed = NULL, # Seed is already set once at the top of exp_kmeans(). Skip it.
+                                            dist = shared_sil_dist) # Reuse the single dist when available.
   }
 
   ret <- ret %>% dplyr::mutate(model = purrr::imap(model, function(x, idx) {
     tryCatch({
       y <- kmeans_model_df$model[[idx]]
       x$kmeans <- y # Might need to be more careful on guaranteeing x and y are from same group, but we are not supporting group_by on UI at this point.
+      # Always attach per-row silhouette for the chosen clustering (#36106).
+      if (!is.null(shared_sil_mat)) {
+        x$silhouette <- compute_silhouette_per_row(y$cluster, shared_sil_mat, shared_sil_dist)
+      } else {
+        # Grouped case: build the matrix from this group's own data.
+        gmat <- as_numeric_matrix_(x$df[, selected_cols, drop = FALSE], columns = selected_cols)
+        if (normalize_data) {
+          gmat <- scale(gmat)
+          gmat[is.nan(gmat)] <- 0
+        }
+        x$silhouette <- compute_silhouette_per_row(y$cluster, gmat)
+      }
       x$sampled_nrow <- sampled_nrow
       x$excluded_nrow <- excluded_nrow
       if (!is.null(elbow_result)) {
diff --git a/R/prcomp.R b/R/prcomp.R
@@ -175,6 +175,23 @@ tidy.prcomp_exploratory <- function(x, type="variances", n_sample=NULL, pretty.n
   }
   else if (type == "summary") { # This is only for kmeans case. TODO: We might want to separate PCA code and k-means code.
     res <- broom::tidy(x$kmeans)
+    if (!is.null(x$silhouette)) {
+      # Per-cluster silhouette aggregates keyed by the same cluster labels as broom::tidy(kmeans).
+      sil_summary <- x$silhouette %>%
+        dplyr::mutate(.cluster_key = as.character(x$kmeans$cluster)) %>%
+        dplyr::group_by(.cluster_key) %>%
+        dplyr::summarise(
+          # Guard the degenerate all-NA case so it yields NA (not NaN/Inf) consistently.
+          avg_silhouette = if (all(is.na(silhouette_score))) NA_real_ else mean(silhouette_score, na.rm = TRUE),
+          min_silhouette = if (all(is.na(silhouette_score))) NA_real_ else min(silhouette_score, na.rm = TRUE),
+          pct_negative = if (all(is.na(silhouette_score))) NA_real_ else mean(silhouette_score < 0, na.rm = TRUE),
+          .groups = "drop"
+        )
+      res <- res %>%
+        dplyr::mutate(.cluster_key = as.character(cluster)) %>%
+        dplyr::left_join(sil_summary, by = ".cluster_key") %>%
+        dplyr::select(-.cluster_key)
+    }
     if (with_excluded_rows) {
       res <- res %>% tibble::add_row(size=x$excluded_nrow)
     }
@@ -195,6 +212,11 @@ tidy.prcomp_exploratory <- function(x, type="variances", n_sample=NULL, pretty.n
       res <- res %>% dplyr::mutate(dplyr::across(dplyr::all_of(column_names), exploratory::normalize))
     }
 
+    if (type == "data" && !is.null(x$silhouette)) {
+      # Bind per-row silhouette (aligned positionally to x$df, same as the cluster column above).
+      res <- res %>% dplyr::bind_cols(x$silhouette)
+    }
+
     if (!is.null(n_sample)) { # default is no sampling.
       # limit n_sample so that no more dots are created than the max that can be plotted on scatter plot, which is 5000.
       n_sample <- min(n_sample, floor(5000 / length(column_names)))
diff --git a/tests/testthat/test_kmeans_1.R b/tests/testthat/test_kmeans_1.R
@@ -138,3 +138,86 @@ test_that("exp_kmeans elbow_method_mode backward compatibility: logical TRUE/FAL
 #  model_df <- exp_kmeans(df, cyl, mpg, hp, elbow_method_mode=TRUE, max_centers=3)
 #  model_df %>% tidyr::unnest(model)
 #})
+
+test_that("compute_silhouette_per_row returns aligned per-row columns", {
+  set.seed(1)
+  mat <- as.matrix(iris[, 1:4])
+  mat <- scale(mat)
+  km <- stats::kmeans(mat, centers = 3)
+  res <- compute_silhouette_per_row(km$cluster, mat)
+  expect_equal(nrow(res), nrow(mat))
+  expect_setequal(colnames(res), c("silhouette_score", "nearest_cluster", "cluster_width"))
+  expect_true(all(res$silhouette_score >= -1 & res$silhouette_score <= 1))
+  expect_true(all(res$nearest_cluster %in% unique(km$cluster)))
+  # cluster_width is the per-cluster mean of silhouette_score, broadcast to rows.
+  expected_avg <- tapply(res$silhouette_score, km$cluster, mean)
+  for (cl in unique(km$cluster)) {
+    expect_equal(unique(res$cluster_width[km$cluster == cl]), unname(expected_avg[as.character(cl)]))
+  }
+  # Passing precomputed dist gives identical result to recomputing it.
+  d <- stats::dist(mat)
+  res_d <- compute_silhouette_per_row(km$cluster, mat, d)
+  expect_equal(res, res_d)
+})
+
+test_that("compute_silhouette_per_row returns all-NA for degenerate input (no error)", {
+  mat <- matrix(rep(1, 20), ncol = 2)            # all identical points
+  res <- compute_silhouette_per_row(rep(1L, 10), mat)  # single cluster
+  expect_equal(nrow(res), 10)
+  expect_true(all(is.na(res$silhouette_score)))
+  expect_true(all(is.na(res$nearest_cluster)))
+  expect_true(all(is.na(res$cluster_width)))
+})
+
+test_that("iterate_silhouette accepts a precomputed dist and matches recompute", {
+  set.seed(1)
+  df <- iris[, 1:4]
+  mat <- scale(as_numeric_matrix_(df, columns = colnames(df)))
+  d <- stats::dist(mat)
+  set.seed(1); a <- iterate_silhouette(df, max_centers = 4, normalize_data = TRUE, seed = 1)
+  set.seed(1); b <- iterate_silhouette(df, max_centers = 4, normalize_data = TRUE, seed = 1, dist = d)
+  expect_equal(a, b)
+})
+
+test_that("exp_kmeans attaches per-row silhouette to each model (all elbow modes)", {
+  df <- mtcars
+  for (mode in list("none", "elbow", "silhouette")) {
+    model_df <- exp_kmeans(df, cyl, mpg, hp, centers = 3, elbow_method_mode = mode, max_nrow = 30)
+    model <- model_df$model[[1]]
+    expect_true(!is.null(model$silhouette))
+    expect_equal(nrow(model$silhouette), nrow(model$df))
+    expect_setequal(colnames(model$silhouette),
+                    c("silhouette_score", "nearest_cluster", "cluster_width"))
+  }
+})
+
+test_that("tidy_rowwise summary includes per-cluster silhouette aggregates", {
+  df <- mtcars
+  model_df <- exp_kmeans(df, cyl, mpg, hp, centers = 3, max_nrow = 30)
+  res <- model_df %>% tidy_rowwise(model, type = "summary")
+  expect_true(all(c("avg_silhouette", "min_silhouette", "pct_negative") %in% colnames(res)))
+  non_excluded <- res[!is.na(res$cluster), ]
+  expect_true(all(non_excluded$pct_negative >= 0 & non_excluded$pct_negative <= 1))
+  expect_true(all(non_excluded$avg_silhouette >= non_excluded$min_silhouette))
+})
+
+test_that("tidy_rowwise data includes per-row silhouette columns", {
+  df <- mtcars
+  model_df <- exp_kmeans(df, cyl, mpg, hp, centers = 3, max_nrow = 30)
+  res <- model_df %>% tidy_rowwise(model, type = "data")
+  expect_true(all(c("silhouette_score", "nearest_cluster", "cluster_width") %in% colnames(res)))
+  expect_true(all(res$silhouette_score >= -1 & res$silhouette_score <= 1))
+  # gathered_data must NOT carry the per-row silhouette columns (charts select their own cols).
+  g <- model_df %>% tidy_rowwise(model, type = "gathered_data")
+  expect_false("silhouette_score" %in% colnames(g))
+})
+
+test_that("exp_kmeans with strange column name still yields silhouette columns", {
+  df <- mtcars
+  df <- df %>% dplyr::rename(`Cy l !#$%` = cyl)
+  model_df <- exp_kmeans(df, `Cy l !#$%`, mpg, hp, centers = 3, max_nrow = 30)
+  res <- model_df %>% tidy_rowwise(model, type = "data")
+  expect_true("silhouette_score" %in% colnames(res))
+  smry <- model_df %>% tidy_rowwise(model, type = "summary")
+  expect_true("avg_silhouette" %in% colnames(smry))
+})
