Unsupervised Clustering: Five Major Families

Chapter 23 introduced K-means with a six-patient example, then ran PCA and K-means on the full synthetic four-variable table; this chapter starts from that K-means and surveys the other clustering families. Clustering is descriptive. A cluster is a group defined by an algorithm and chosen distance/likelihood assumptions, not a confirmed disease class. Before clustering, decide what counts as similar and whether rows represent people, visits, or time points.

1. Use features without the outcome

We use invented ages, BMIs, SBP, and pulse rates. We deliberately exclude patient_id and diabetes from the clustering inputs. The numeric values have different units, so standardization is important for distance-based methods.

library(tidymodels)
ehr <- readr::read_csv("data/synthetic_ehr.csv", show_col_types = FALSE)
x_raw <- ehr |> select(age, bmi, sbp, pulse)
prep_plan <- recipe(~ age + bmi + sbp + pulse, data = x_raw) |>
  step_impute_median(all_numeric_predictors()) |>
  step_normalize(all_numeric_predictors())
x <- as.matrix(bake(prep(prep_plan, training = x_raw), new_data = x_raw))
stopifnot(!anyNA(x))
dim(x)

[1] 600   4

round(colMeans(x), 2)

  age   bmi   sbp pulse 
    0     0     0     0 

The matrix has 600 rows and 4 columns, and no missing values remain. The column means are 0 because each variable was centered, which is what standardizing does.

Scope: This is a descriptive exercise, not a predictive validation. If a clustering result becomes an input to a supervised prediction task, the preprocessing and cluster assignment rules must be learned inside training folds and applied without access to the held-out set.

2. K-means: divide observations into k groups

This starts from the same standardized four-variable table (age, bmi, sbp, pulse) that chapter 23 clustered with k = 2, and now tries k = 3. Each point goes to a cluster whose centroid (average point) is closest. The algorithm tries to reduce within-cluster squared distance. You must specify k (section 9 shows how to choose it). It tends to favor relatively compact, approximately spherical groups and is sensitive to outliers.

set.seed(500)
km <- kmeans(x, centers = 3, nstart = 25)
table(km$cluster)

  1   2   3 
189 237 174 

round(km$centers, 2)

    age   bmi   sbp pulse
1  0.72  0.82  0.99  0.11
2 -0.99  0.03 -0.86 -0.13
3  0.57 -0.93  0.10  0.06

nstart=25 repeats initialization because K-means may settle in a local optimum. The cluster numbers are arbitrary IDs, not ranks. k=3 does not imply three biologically real patient types.

3. What is inside a kmeans() result

kmeans() returns a list. str() lists every part, so you can see what you can pull out and use later.

str(km)

List of 9
 $ cluster     : int [1:600] 3 2 3 2 1 2 1 1 2 1 ...
 $ centers     : num [1:3, 1:4] 0.7204 -0.9933 0.5704 0.8198 0.0299 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:3] "1" "2" "3"
  .. ..$ : chr [1:4] "age" "bmi" "sbp" "pulse"
 $ totss       : num 2396
 $ withinss    : num [1:3] 418 591 354
 $ tot.withinss: num 1363
 $ betweenss   : num 1033
 $ size        : int [1:3] 189 237 174
 $ iter        : int 3
 $ ifault      : int 0
 - attr(*, "class")= chr "kmeans"

• cluster: the group number (1 to k) of each of the 600 rows, in the original row order.
• centers: the average (centroid) of each variable in each group. The input was standardized, so these are in standard-deviation units, not years or mmHg.
• totss: total sum of squares, the spread of all points around the overall mean.
• withinss: one number per group, the spread of its points around its own centroid.
• tot.withinss: the sum of withinss. K-means tries to make this small, and it is the number the elbow plot uses.
• betweenss: the spread of the group centroids around the overall mean, the part of the spread the grouping explains.
• size: how many rows are in each group.
• iter: how many rounds the algorithm needed.
• ifault: a code for trouble; 0 means no problem was flagged.

The parts add up: totss equals tot.withinss plus betweenss.

c(totss = km$totss, within_plus_between = km$tot.withinss + km$betweenss)

              totss within_plus_between 
               2396                2396 

all.equal(km$totss, km$tot.withinss + km$betweenss)

[1] TRUE

round(km$betweenss / km$totss, 2)

[1] 0.43

all.equal() returns TRUE, so the identity holds. The last line is the share of the spread that the three groups explain. A small share means the groups separate the points only a little.

4. Reproducibility: set.seed() and nstart

K-means starts from randomly chosen centers, so two runs can end in different groupings. set.seed() makes the random choice repeatable, and nstart runs many random starts and keeps the best one (the lowest tot.withinss). Run the algorithm six times with one start each and no seed:

runs_no_seed <- replicate(6, kmeans(x, centers = 3, nstart = 1)$tot.withinss)
round(runs_no_seed, 1)

[1] 1389.9 1370.7 1372.3 1364.6 1364.6 1372.3

The six values usually differ, and they will not match this page when you run the code, because nothing fixes the random start. Now set the seed before each run, and use 25 starts:

run_seeded <- function() {
  set.seed(500)
  kmeans(x, centers = 3, nstart = 25)$tot.withinss
}
round(c(run_seeded(), run_seeded()), 1)

[1] 1362.9 1362.9

The two numbers are identical, so anyone with the same data and code gets your result. Always set a seed before random steps and report it. Use nstart of 10 or more so that the result does not depend on one lucky or unlucky start.

5. K-medoids / PAM: use representative observations

PAM (partitioning around medoids) chooses observed data points as cluster centers instead of mean vectors. This can be less sensitive to extreme points under some distances but is computationally heavier on large datasets.

pam_model <- cluster::pam(x, k = 3)
table(pam_model$clustering)

  1   2   3 
176 245 179 

pam_model$medoids

            age        bmi         sbp         pulse
[1,]  0.3949913 -0.9813892 -0.01619059 -0.3721749510
[2,] -0.9706702  0.3538013 -0.64292293  0.0009327693
[3,]  1.1828730  0.5569825  1.08059102 -0.0923441608

PAM puts 176, 245 and 179 rows in its three groups, and each row of medoids is one real record in standardized units (for example, the second medoid is a younger record, age -0.97, with lower SBP). Compare K-means and PAM qualitatively: do they group the same records? Which points change membership? They optimize different objectives.

6. Hierarchical clustering: build a tree of similarities

Hierarchical clustering builds nested groups at many levels. Agglomerative methods merge close groups from small to large; divisive methods split groups from large to small. Linkage options (single, complete, average, Ward) define what “two groups are close” means.

For a readable dendrogram, use 40 fabricated records, not all 600:

d <- dist(x[1:40, ], method = "euclidean")
hc <- hclust(d, method = "ward.D2")
plot(hc, labels = FALSE, main = "Hierarchical similarity (40 FAKE rows)")

↔ Swipe the plot sideways to see all of it.

groups <- cutree(hc, k = 3)
table(groups)

groups
 1  2  3 
18 14  8 

The branches show similarity under the chosen metric, not ancestry, progression, or medical severity.

7. Density-based: DBSCAN, HDBSCAN, OPTICS

Some datasets contain irregularly shaped groups and isolated points. Density-based algorithms group locally dense areas and may label sparse regions noise.

• DBSCAN: parameters eps (neighborhood radius) and minPts (minimum number of points in a neighborhood). Noise often receives cluster ID 0.
• HDBSCAN: studies clusters across density scales; can handle varying densities better in some settings.
• OPTICS: produces density-based ordering across neighborhood scales.

Optional R demonstration (requires package dbscan, not part of core Quarto rendering):

# install.packages("dbscan")
library(dbscan)
db <- dbscan(x, eps = 1.0, minPts = 8)
table(db$cluster)       # 0 indicates noise

Neither “noise” nor a rare cluster means abnormal disease; it may simply reflect the distance metric or sampling.

Example: the DBSCAN rule can be followed by hand on ten typed numbers, with no package. Take eps = 0.5 and minPts = 3. A point is a core point if at least 3 points (itself included) lie within 0.5 of it; a point that is not core but lies within 0.5 of a core point is a border point; the rest is noise.

pts <- c(1.0, 1.2, 1.4, 1.5, 5.0, 9.0, 9.1, 9.3, 9.2, 9.4)
eps <- 0.5
min_pts <- 3
near <- as.matrix(dist(pts)) <= eps          # TRUE where two points are within eps
n_near <- rowSums(near)
core <- n_near >= min_pts
near_core <- apply(near, 1, function(r) any(r & core))
tibble(point = pts, neighbors = n_near,
       role = ifelse(core, "core", ifelse(near_core, "border", "noise")))

# A tibble: 10 x 3
   point neighbors role 
   <dbl>     <dbl> <chr>
 1   1           4 core 
 2   1.2         4 core 
 3   1.4         4 core 
 4   1.5         4 core 
 5   5           1 noise
 6   9           5 core 
 7   9.1         5 core 
 8   9.3         5 core 
 9   9.2         5 core 
10   9.4         5 core 

The four points near 1 and the five points near 9 are all core points, so they form two dense groups. The lone point at 5.0 has only itself within 0.5, and no core point is near it, so it is noise and would get cluster ID 0. K-means would have forced it into one of the two groups.

8. Model-based, graph, and spectral clustering

• Gaussian mixture models (GMM): model data as a mixture of probability distributions; allow soft membership probabilities. Example R package: mclust.
• Latent-class analysis: categorical-data mixtures with conditional response models; example poLCA.
• Spectral clustering: uses the eigenvectors of a similarity graph before partitioning; example kernlab::specc().
• Community detection: partitions a network graph, such as a constructed patient-similarity network; example igraph. Graph construction and privacy matter a great deal here.

These algorithms assume different shapes or mathematical structures. Ask what kind of data and clusters are scientifically plausible.

Example: a Gaussian mixture gives each point a probability of belonging to each group instead of a hard label. Suppose two groups have means 3 and 7, both with standard deviation 1, and equal size. The membership probability of a point is its density under one group divided by the sum of both densities.

x0 <- c(4, 5, 6.5)
d_group1 <- dnorm(x0, mean = 3, sd = 1)
d_group2 <- dnorm(x0, mean = 7, sd = 1)
round(cbind(x = x0, group1 = d_group1 / (d_group1 + d_group2),
            group2 = d_group2 / (d_group1 + d_group2)), 3)

       x group1 group2
[1,] 4.0  0.982  0.018
[2,] 5.0  0.500  0.500
[3,] 6.5  0.002  0.998

A point at 4 belongs to group 1 with probability 0.98, a point at 6.5 to group 2 with 0.998, and a point at 5, exactly halfway, is 0.5 for each. K-means would give all three a single group number and hide that the point at 5 is a coin flip. The mclust package fits the means, spreads and sizes from data.

9. Choosing the number of clusters: diagnostics versus truth

Useful internal checks include within-cluster sum of squares (elbow), average silhouette width (how much closer a point is to its own group than to the nearest other group), gap statistic, stability under resampling, and (for mixture models) AIC/BIC. No single score confirms a clinically meaningful number of groups.

Elbow. Plot tot.withinss for several k. It always falls as k grows, so look for the k after which the drop becomes small.

k_values <- 2:6
set.seed(500)
wcss <- sapply(k_values, function(k) {
  kmeans(x, centers = k, nstart = 10)$tot.withinss
})
plot(k_values, wcss, type = "b",
     xlab = "Number of groups", ylab = "Within-cluster sum of squares")

↔ Swipe the plot sideways to see all of it.

A curve without a clear elbow does not justify inventing one.

Average silhouette width. For each row, cluster::silhouette() compares the distance to its own group with the distance to the nearest other group. The value runs from -1 to 1: near 1 means clearly in its group, near 0 means on the border, below 0 means probably in the wrong group. Higher average is better. It needs the cluster labels and a distance matrix.

d_all <- dist(x)
set.seed(500)
sil_avg <- sapply(k_values, function(k) {
  cl <- kmeans(x, centers = k, nstart = 10)$cluster
  mean(cluster::silhouette(cl, d_all)[, "sil_width"])
})
round(setNames(sil_avg, k_values), 3)

    2     3     4     5     6 
0.276 0.232 0.235 0.246 0.243 

The highest average is at k = 2 (0.276), and every value is far below 0.5, so no choice of k gives clearly separated groups. Rough guide (Kaufman and Rousseeuw): 0.26 to 0.50 is weak structure, and 0.25 or below means no substantial structure. The k = 3 value (0.232) is the lowest, which is one more reason not to read k = 3 as a real number of types.

Silhouette plot for one k. One bar per row, grouped by cluster and sorted, shows which groups are tight and which rows sit on a border. Here it is for the k = 3 model from section 2.

sil3 <- cluster::silhouette(km$cluster, d_all)
plot(sil3, col = "grey60", main = "", border = NA)

↔ Swipe the plot sideways to see all of it.

c(rows_below_zero = sum(sil3[, "sil_width"] < 0),
  smallest_width = round(min(sil3[, "sil_width"]), 3),
  largest_width = round(max(sil3[, "sil_width"]), 3))

rows_below_zero  smallest_width   largest_width 
         14.000          -0.059           0.459 

The number printed beside each group is its size and its average silhouette width, about 0.23 to 0.24 for all three. The bars are short (the largest width, printed after the plot, is about 0.46, so none reaches 0.5), which matches the low average. The first two numbers printed after the plot count the rows with a negative width (14 of 600 rows) and give the lowest width, about -0.06. Those rows are closer to another group than to their own, but only slightly, so they are rows on a border and not clear mistakes.

Gap statistic. cluster::clusGap() compares the within-group spread for each k with the spread you would get from random data with no groups. B is the number of random reference data sets; a small B keeps this fast, and 50 or more is better for real work.

set.seed(500)
gap <- cluster::clusGap(x, FUN = kmeans, K.max = 6, B = 20, nstart = 5)
round(gap$Tab[, c("gap", "SE.sim")], 3)

       gap SE.sim
[1,] 0.379  0.011
[2,] 0.343  0.010
[3,] 0.337  0.009
[4,] 0.340  0.011
[5,] 0.352  0.010
[6,] 0.350  0.009

cluster::maxSE(gap$Tab[, "gap"], gap$Tab[, "SE.sim"], method = "firstSEmax")

[1] 1

A larger gap means the grouping is further from what random data gives. Here gap is positive at every k (0.34 to 0.38), so the data are somewhat more clustered than random data. maxSE() with "firstSEmax" first finds where the gap stops rising, the first local maximum, and then returns the smallest k whose gap is within one standard error of that maximum. The gap already falls from 0.379 at k = 1 to 0.343 at k = 2, so the first local maximum is k = 1 and the function returns 1. So the gap statistic does not pick out any k of 2 or more, which agrees with the weak elbow and low silhouette values. Report the elbow and silhouette (and gap if you use it) together, and say plainly when they do not agree or are weak. The factoextra package has one-call versions: fviz_nbclust() draws these three checks and fviz_cluster() draws the cluster view in the next section.

Also compare groups against site/batch effects: clusters may reflect measurement process, not patient health.

10. Looking at the clusters: a PCA view

Four variables cannot be drawn on one flat plot. A cluster view projects the rows onto the first two principal components (PCA, introduced in chapter 23; chapter 25 goes deeper) and colors each point by its group. This is what fviz_cluster() draws.

pc <- prcomp(x)
view <- data.frame(PC1 = pc$x[, 1], PC2 = pc$x[, 2],
                   cluster = factor(km$cluster))
ggplot(view, aes(PC1, PC2, colour = cluster)) +
  geom_point(alpha = 0.6, size = 1) +
  labs(colour = "Cluster")

↔ Swipe the plot sideways to see all of it.

If the colors form separate islands, the groups are clearly apart; if they are neighboring slices of one cloud, as here, the groups are a convenient partition and not natural types.

11. Describing the clusters: a profile plot

A cluster profile answers “what is each group like?” Plot each group’s standardized mean for every variable, one line per group.

profile <- as.data.frame(km$centers) |>
  mutate(cluster = factor(row_number())) |>
  pivot_longer(-cluster, names_to = "variable", values_to = "mean_z") |>
  mutate(variable = factor(variable, levels = colnames(x)))
ggplot(profile, aes(variable, mean_z, colour = cluster, group = cluster)) +
  geom_line() +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = NULL, y = "Mean (standard-deviation units)", colour = "Cluster")

↔ Swipe the plot sideways to see all of it.

Zero (dashed) is the overall average; a point at 1 means the group’s mean is one standard deviation above it. Read each line as a description: for example, “group 2 is younger and has lower SBP”. These are descriptions of similarity groups, not diagnoses. The GGally package draws a related view for individual rows with ggparcoord().

12. Stability and downstream safety

To assess stability, rerun with different starts/subsamples and compare membership with methods such as adjusted Rand index, but understand that cluster labels may permute. To name a cluster “high risk,” independent verified outcomes and evaluation are needed. Clustering alone never labels a person as diseased.

Example: two k-means runs with different seeds can give the same groups under different numbers. A cross-table of the labels shows this.

set.seed(42)
km_b <- kmeans(x, centers = 3, nstart = 25)
table(first_run = km$cluster, second_run = km_b$cluster)

         second_run
first_run   1   2   3
        1   0 189   0
        2 235   0   2
        3   0   0 174

Each row has one large cell (189, 235 and 174 rows) and the rest are near zero, apart from 2 rows that moved. So the two runs found nearly the same groups and only the numbering changed (first-run group 1 is second-run group 2, 2 is 1, and 3 is 3); cells spread over a row would mean unstable groups. The adjusted Rand index summarizes this table as one number that ignores the numbering.

Runnable lab: the included script at scripts/06_unsupervised_comparison.R demonstrates K-means, PAM, hierarchical groups, PCA, anomaly scores, and association-rule calculations with fake data.

Try it

1. Fit K-means at k=2, k=3, and k=4. Are cluster sizes equal?

Show solution to question 1

Answer: No. The groups differ in size at every k, and asking for more groups does not make them equal.

Step by step

1. Fit each model with the same seed and nstart, and keep the size part of the result.

set.seed(500)
sizes <- lapply(2:4, function(k) kmeans(x, centers = k, nstart = 25)$size)
names(sizes) <- paste("k =", 2:4)
sizes

Output

$`k = 2`
[1] 307 293

$`k = 3`
[1] 189 235 176

$`k = 4`
[1] 161 155 150 134

2. Each list entry holds the group sizes in cluster-number order. The sizes within each k are unequal, and K-means has no rule that forces equal sizes. The k = 3 sizes (189, 235, 176) differ slightly from section 2 (189, 237, 174) because k = 2 runs first here, so the random draws are different.

2. What is the difference between a centroid and a medoid?

Show solution to question 2

Answer: A centroid is the average of the group’s points and need not be a real observation. A medoid is one of the actual observations, the one most central in its group.

Step by step

1. Take group 1 of the K-means fit (centroid) and group 1 of the PAM fit (medoid) from earlier sections.
2. Check whether each center equals a row of the data.

is_a_row <- function(center) any(apply(x, 1, function(r) all(abs(r - center) < 1e-8)))
c(centroid_is_a_row = is_a_row(km$centers[1, ]),
  medoid_is_a_row = is_a_row(pam_model$medoids[1, ]))

Output

centroid_is_a_row   medoid_is_a_row 
            FALSE              TRUE 

3. The centroid is not a row (FALSE), because averaging produces a point nobody was measured at. The medoid is a row (TRUE), so it can be shown as a real example patient.

3. Which method can identify irregularly shaped dense regions and explicitly mark noise?

Show solution to question 3

Answer: Density-based clustering: DBSCAN (and HDBSCAN or OPTICS).

Step by step

1. K-means and PAM assign every point to a group, and the groups tend to be compact blobs, so they cannot mark noise.
2. DBSCAN grows groups from dense neighborhoods (eps, minPts), so a group can take any shape, and points in no dense neighborhood stay unassigned.
3. The by-hand example in section 7 shows this: the point at 5.0 has the role “noise” while the other nine points form two dense groups.

4. Why could clustering separate clinic sites instead of medical states?

Show solution to question 4

Answer: Clusters only reflect whatever differs most between rows. If sites use different devices, units or patient populations, those differences can drive the groups even when health is the same.

Step by step

1. Cross the k = 3 groups with the clinic site.

table(cluster = km$cluster, site = ehr$clinic_site)

Output

       site
cluster Mobile North South
      1     69    62    58
      2     69    87    81
      3     62    51    61

2. If one site filled most of one cluster, the cluster would describe the site. Read the table row by row: here each cluster contains rows from all sites, so on this data the clusters do not follow the sites. 3. That check is the one to run on real data. A cluster that matches a site or device points to a measurement or workflow cause before any medical one.

5. Explain why a silhouette score does not prove diagnostic accuracy.

Show solution to question 5

Answer: The silhouette measures only how well separated the groups are in the chosen distance. It never uses a confirmed diagnosis, so a high value says nothing about whether the groups match a disease.

Step by step

1. The inputs to silhouette() are the cluster labels and the distance matrix. No outcome column is involved.
2. As a check, compare the diabetes rate across the three groups.

round(tapply(ehr$diabetes == "Yes", km$cluster, mean), 3)

Output

    1     2     3 
0.206 0.030 0.046 

3. Group 1 has a diabetes rate of 0.206 against 0.030 and 0.046 in the others, even though the groups were made without looking at diabetes. That is a pattern worth noticing, but it is not a diagnosis or a measure of accuracy: the silhouette says nothing about it, and confirming it needs verified outcomes and evaluation on new data.

6. Fit km2 <- kmeans(x, centers = 4, nstart = 25) after set.seed(1). Show that totss equals tot.withinss plus betweenss, and say which two parts size and withinss describe.

Show solution to question 6

Answer: all.equal() returns TRUE, so totss = tot.withinss + betweenss. size is the number of rows in each group and withinss is how spread out each group is around its centroid.

Step by step

1. Fit the model with the seed set first.
2. Compare totss with the sum of the other two parts, and print size and withinss.

set.seed(1)
km2 <- kmeans(x, centers = 4, nstart = 25)
all.equal(km2$totss, km2$tot.withinss + km2$betweenss)

Output

[1] TRUE

km2$size

Output

[1] 134 150 161 155

round(km2$withinss, 1)

Output

[1] 230.1 291.4 324.0 294.9

3. TRUE confirms the identity. Both size and withinss have one number per group (four each): group counts, and each group’s own spread. Their sum over groups is tot.withinss.

7. Run kmeans(x, centers = 3, nstart = 1) five times without a seed. Do the tot.withinss values agree? Then repeat with set.seed(1) before each run.

Show solution to question 7

Answer: Without a seed the values usually differ, because each run starts from random centers. With set.seed(1) before each run they are identical.

Step by step

1. Run five unseeded fits and print the totals. These are random, so your numbers will differ from the page.

round(replicate(5, kmeans(x, centers = 3, nstart = 1)$tot.withinss), 1)

Output

[1] 1372.6 1364.2 1364.7 1364.6 1364.2

2. Repeat with the seed set inside each repeat.

round(replicate(5, { set.seed(1); kmeans(x, centers = 3, nstart = 1)$tot.withinss }), 1)

Output

[1] 1362.9 1362.9 1362.9 1362.9 1362.9

3. The second line shows five equal values, so anyone with the same seed gets the same result. Setting a larger nstart makes unseeded runs agree more often, but only the seed guarantees it.

8. Compute the average silhouette width for k = 2 to 4 and name the k with the highest value. Is that value high enough to call the groups well separated?

Show solution to question 8

Answer: The highest average silhouette width is at k = 2 (see the printed values), and it is far below 0.5, so the groups are not well separated at any of these k.

Step by step

1. For each k, fit K-means, compute the silhouette of its labels against the distance matrix d_all, and average the widths.

set.seed(1)
sil_k <- sapply(2:4, function(k) {
  cl <- kmeans(x, centers = k, nstart = 10)$cluster
  mean(cluster::silhouette(cl, d_all)[, "sil_width"])
})
round(setNames(sil_k, 2:4), 3)

Output

    2     3     4 
0.276 0.232 0.235 

2. Pick the largest value: k = 2 with 0.276 (k = 3 gives 0.232 and k = 4 gives 0.235). 3. Compare with the guide in section 9: 0.26 to 0.50 is weak structure and 0.25 or below is no substantial structure. Even the best value is barely in the weak range, so the groups are a convenient partition and not clearly separated.

9. Make the cluster view (first two principal components) and the profile plot for your k = 4 model. Describe one group in words, without using a disease name.

Show solution to question 9

Answer: The groups overlap as neighboring slices of one cloud, so they are a convenient partition. One description: group 2’s average patient is younger with a fast pulse (0.95 standard deviations above average) and lower SBP, with BMI a little below average.

Step by step

1. Print the group centers: one row per group, in standard-deviation units.

round(km2$centers, 2)

Output

    age   bmi   sbp pulse
1  0.87 -0.80  0.39 -0.49
2 -0.60 -0.35 -0.62  0.95
3 -0.93  0.18 -0.77 -0.82
4  0.80  0.85  1.06  0.35

2. Project the rows onto the first two principal components and color by km2$cluster.

view4 <- data.frame(PC1 = pc$x[, 1], PC2 = pc$x[, 2], cluster = factor(km2$cluster))
ggplot(view4, aes(PC1, PC2, colour = cluster)) +
  geom_point(alpha = 0.6, size = 1)

↔ Swipe the plot sideways to see all of it.

3. Plot km2$centers, one line per group.

prof4 <- as.data.frame(km2$centers) |>
  mutate(cluster = factor(row_number())) |>
  pivot_longer(-cluster, names_to = "variable", values_to = "mean_z") |>
  mutate(variable = factor(variable, levels = colnames(x)))
ggplot(prof4, aes(variable, mean_z, colour = cluster, group = cluster)) +
  geom_line() +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = NULL, y = "Mean (standard-deviation units)", colour = "Cluster")

↔ Swipe the plot sideways to see all of it.

4. The cluster view shows no separate islands, so the four groups are not clearly apart. In the printed centers and the profile, read group 2’s row: age -0.60, sbp -0.62 and bmi -0.35 (all below the dashed zero line) and pulse +0.95. That gives “younger, lower SBP, higher pulse”, a description of similarity and not a diagnosis.

10. Look at the gap table in section 9: at which k is the gap largest, and what does maxSE() return? Does it agree with the silhouette values?

Show solution to question 10

Answer: The gap is largest at k = 1 (0.379), and maxSE() with "firstSEmax" returns 1. It agrees with the silhouette values in saying that no k of 2 or more is clearly supported.

Step by step

1. Print the gap table from section 9 and find the largest gap.

gap_tab <- round(gap$Tab[, c("gap", "SE.sim")], 3)
gap_tab

Output

       gap SE.sim
[1,] 0.379  0.011
[2,] 0.343  0.010
[3,] 0.337  0.009
[4,] 0.340  0.011
[5,] 0.352  0.010
[6,] 0.350  0.009

which.max(gap_tab[, "gap"])

Output

[1] 1

cluster::maxSE(gap$Tab[, "gap"], gap$Tab[, "SE.sim"], method = "firstSEmax")

Output

[1] 1

2. The largest gap is in row 1 (k = 1). The gap falls at k = 2, so the first local maximum is k = 1 and maxSE() returns 1. 3. The silhouette widths were all below 0.5 and weak, and the elbow had no clear bend. All three checks say the same: the data show no clear number of groups, and k = 3 or 4 is a choice made for convenience.