flat clustering

Flat clustering uses the concept of distance between samples in variable space to assign each sample in a dataset to one group, either with k-means, which divides the samples into a chosen number (k) of clusters around cluster means, or with dbscan, which forms clusters from regions where samples are densely packed and labels samples in sparse regions as noise, after which the cluster labels can be summarized and joined to other data. This analysis can be used to answer questions such as “Which of k clusters does each sample belong to?”, “At what value of k does adding another cluster stop substantially reducing the within-cluster variance?”, or “Which samples does dbscan label as noise at these settings?”

Do my samples fall into definable clusters? Sorting samples into groups is often the first step in making sense of a large set of measurements: fragments into kinds of meteorite, water samples into sources that share a chemistry, or plant extracts into chemotypes. Once samples have groups, those groups can be counted and summarized, compared with the statistical tests in the comparing groups chapter, or used to color an ordination plot such as a PCA. Some samples are interesting precisely because they fit no group at all. In the hierarchical clustering chapter we answered a related question by building a tree and then cutting it. This chapter sorts samples straight into groups without building a tree first, using two methods, k-means and dbscan, and shows how to decide how many groups to ask for and when a sample is better left out of every group.

kmeans

K-means is the first of the two methods. It sorts samples into however many clusters we ask for, placing each sample in the cluster whose center it sits closest to.

For the examples in this chapter we will use meteorite_fragments, which loads along with the rest of the course datasets. It describes 37 meteorite fragments recovered from four places where meteorite hunters search: the Atacama and Sahara deserts, the Nullarbor Plain in Australia, and the Antarctic ice fields. Each fragment was measured for its iron content (iron_percent), its density (density_g_per_cm3), and two more measurements we will come back to at the end of the chapter, nickel_percent and silica_percent. (These are simulated data, built with realistic values for each kind of meteorite.)

meteorite_fragments
## # A tibble: 37 × 6
##    fragment find_site         iron_percent density_g_per_cm3
##    <chr>    <chr>                    <dbl>             <dbl>
##  1 MF-01    Antarctic ice fi…         24.5              3.48
##  2 MF-02    Atacama Desert            23                3.45
##  3 MF-03    Antarctic ice fi…         37.5              4.08
##  4 MF-04    Nullarbor Plain           22.2              3.51
##  5 MF-05    Sahara Desert             85.8              7.8 
##  6 MF-06    Nullarbor Plain           19.6              3.39
##  7 MF-07    Nullarbor Plain           25                3.37
##  8 MF-08    Sahara Desert             89                7.83
##  9 MF-09    Sahara Desert             19.5              3.55
## 10 MF-10    Atacama Desert            91.9              7.67
## # ℹ 27 more rows
## # ℹ 2 more variables: nickel_percent <dbl>,
## #   silica_percent <dbl>

Before clustering anything, it is worth simply plotting two of the measurements against each other:

ggplot(meteorite_fragments, aes(x = iron_percent, y = density_g_per_cm3)) +
  geom_point(shape = 21, size = 4, fill = "grey70") +
  theme_classic() +
  labs(x = "iron (%)", y = "density (g/cm³)")

Iron content and density of 37 meteorite fragments. A scatter plot with iron content on the horizontal axis and density on the vertical axis; the points fall into three separate clumps, at low, middle and high iron content and density, with two points lying between the clumps.

Three groups stand out right away, plus a couple of fragments that do not sit comfortably in any of them. K-means does formally what the eye just did. Let’s try it using runMatrixAnalysis. We set analysis = "kmeans", name the columns to cluster on, and give the number of clusters we want to the parameters argument as n_clusters (the same way we asked for a number of groups when cutting a tree in the hierarchical clustering chapter). Below we ask for three clusters, and the next section shows how to arrive at that number. (With parameters left out altogether, desktop R opens a small application for choosing the number of clusters with a slider. That application cannot run in a browser, so in a markdown document, in a WebR cell, or in an escape room, the number of clusters always has to be given explicitly.)

meteorite_clusters <- runMatrixAnalysis(
    data = meteorite_fragments,
    analysis = c("kmeans"),
    parameters = c(n_clusters = 3),
    columns_w_values_for_single_analyte = c("iron_percent", "density_g_per_cm3"),
    columns_w_sample_ID_info = c("fragment", "find_site")
)

meteorite_clusters
## # A tibble: 37 × 6
##    sample_unique_ID  fragment find_site cluster iron_percent
##    <chr>             <chr>    <chr>     <chr>          <dbl>
##  1 MF-01_Antarctic … MF-01    Antarcti… cluste…         24.5
##  2 MF-02_Atacama De… MF-02    Atacama … cluste…         23  
##  3 MF-03_Antarctic … MF-03    Antarcti… cluste…         37.5
##  4 MF-04_Nullarbor … MF-04    Nullarbo… cluste…         22.2
##  5 MF-05_Sahara Des… MF-05    Sahara D… cluste…         85.8
##  6 MF-06_Nullarbor … MF-06    Nullarbo… cluste…         19.6
##  7 MF-07_Nullarbor … MF-07    Nullarbo… cluste…         25  
##  8 MF-08_Sahara Des… MF-08    Sahara D… cluste…         89  
##  9 MF-09_Sahara Des… MF-09    Sahara D… cluste…         19.5
## 10 MF-10_Atacama De… MF-10    Atacama … cluste…         91.9
## # ℹ 27 more rows
## # ℹ 1 more variable: density_g_per_cm3 <dbl>

The output has one row per fragment: its ID columns, the cluster k-means assigned it to, and the two measurements it was clustered on. As in the previous chapters, runMatrixAnalysis scales each column before it clusters (scale_variance = TRUE is the default). That matters here: iron runs from about 20 to over 90 percent while density spans only about 3 to 8 g/cm³, so without scaling, the iron percentages would drive the clusters almost on their own. We can plot the results and color them according to the group that k-means suggested. We can also highlight groups using geom_mark_ellipse. Note that it is recommended to specify both fill and label for geom_mark_ellipse:

ggplot() +
  geom_point(
    data = meteorite_clusters,
    aes(x = iron_percent, y = density_g_per_cm3, fill = cluster), shape = 21, size = 4, alpha = 0.6
  ) +
  geom_mark_ellipse(
    data = meteorite_clusters,
    aes(x = iron_percent, y = density_g_per_cm3, label = cluster, fill = cluster), alpha = 0.2
  ) +
  theme_classic() +
  labs(x = "iron (%)", y = "density (g/cm³)") +
  scale_fill_manual(values = discrete_palette)

The meteorite fragments colored by k-means cluster, with an ellipse drawn around each cluster. The three clumps from the previous plot are each one cluster; the two in-between fragments have each been placed in a neighboring cluster.

The three clusters match three real kinds of meteorite. Stony meteorites are mostly silicate rock with a little metal mixed in, iron meteorites are almost solid iron-nickel metal (and so are more than twice as dense), and the rarer stony-iron meteorites are roughly half and half. Notice too what happened to the two fragments that sat between the clumps: k-means put each of them into a cluster anyway. K-means has no way to say “this one does not belong to any group”, and that is one of the things the dbscan method below does differently.

concept check

Fill in the blank with the analysis that sorts the fragments into a set number of clusters. Press Run.

meteorite_clusters <- runMatrixAnalysis(
  data = meteorite_fragments,
  analysis = "____",
  parameters = c(n_clusters = 3),
  columns_w_values_for_single_analyte = c("iron_percent", "density_g_per_cm3"),
  columns_w_sample_ID_info = c("fragment", "find_site")
)

ggplot(meteorite_clusters, aes(x = iron_percent, y = density_g_per_cm3, fill = cluster)) +
  geom_point(shape = 21, size = 4) +
  theme_classic()

Self-check: runMatrixAnalysis scales each column before k-means runs. Why does that matter for these two columns?

Scaling makes every cluster end up with the same number of fragments Iron spans a much wider range of numbers than density, so without scaling iron would decide the clusters almost on its own k-means cannot run on columns that contain decimal numbers

choosing k

Above we simply asserted that three clusters was the right number. How would we know? Fewer clusters are generally easier to interpret, but too few leave a lot of variance inside each cluster, so the goal is a balance between the two. The elbow method, illustrated in the image at the top of this chapter, is a way of finding that balance. In that image, the same points are clustered into one, two, three, and more clusters (panel A), the within-cluster variance is measured for each number of clusters (B), and those values are plotted against the number of clusters (C). The best balance sits at the “elbow”, the point where the line bends most sharply (the point with the smallest angle theta). The elbow point indicates the number of clusters we can use to get the smallest within-cluster variance before diminishing returns set in. We can ask runMatrixAnalysis for the numbers it needs by setting parameters = "elbow":

meteorite_elbow <- runMatrixAnalysis(
    data = meteorite_fragments,
    analysis = c("kmeans"),
    parameters = "elbow",
    columns_w_values_for_single_analyte = c("iron_percent", "density_g_per_cm3"),
    columns_w_sample_ID_info = c("fragment", "find_site")
)

meteorite_elbow
## # A tibble: 10 × 2
##        k within_cluster_variance
##    <int>                   <dbl>
##  1     1                  72    
##  2     2                  10.0  
##  3     3                   1.62 
##  4     4                   0.732
##  5     5                   0.367
##  6     6                   0.280
##  7     7                   0.218
##  8     8                   0.166
##  9     9                   0.141
## 10    10                   0.118

The elbow follows the same scale_variance setting as the clustering, so the elbow on the plot always belongs to the clusters that could then be requested by changing the parameters argument to c(n_clusters = 3), or whatever number is chosen. This gives us one row for each number of clusters it tried: k, and the within_cluster_variance that number of clusters leaves behind. Notice that the within-cluster variance always falls as we add clusters. At the extreme, where every sample is its own cluster, there is no within-cluster variance left at all. What we are after is the bend:

ggplot(meteorite_elbow, aes(x = k, y = within_cluster_variance)) +
  geom_line() +
  geom_point(size = 3) +
  theme_classic() +
  labs(x = "number of clusters (k)", y = "within-cluster variance")

Elbow plot for the meteorite fragments. Within-cluster variance falls steeply from one to two to three clusters and then levels off, so the bend in the line is at three clusters.

Reading left to right: going from one cluster to two, and from two to three, each buys a large reduction in within-cluster variance, and after that the returns collapse. The line has clearly flattened by three clusters, which is where the three we used above came from. Here the elbow is sharp, because the fragments really do come in three kinds. With real data it is often gentler, and then choosing k is a judgment call, not a calculation. It is an important skill to be able to recognize when the elbow is sharp (and so the number of groups is really a property of the data) and when it is gentle (and so the number of groups is largely a choice). Even a sharp elbow belongs to the columns that were clustered, not to the samples alone: the same samples clustered on a different set of measurements can bend just as sharply at a different number of groups, and both answers can be right.

By default runMatrixAnalysis tries every number of clusters from one up to ten. To try a different range, pass the maximum as a second value:

runMatrixAnalysis(
    data = meteorite_fragments,
    analysis = c("kmeans"),
    parameters = c("elbow", 6),
    columns_w_values_for_single_analyte = c("iron_percent", "density_g_per_cm3"),
    columns_w_sample_ID_info = c("fragment", "find_site")
)

concept check

Fill in the blank so that runMatrixAnalysis returns the within-cluster variance for each number of clusters instead of the clusters themselves. Press Run and find the bend.

meteorite_elbow <- runMatrixAnalysis(
  data = meteorite_fragments,
  analysis = "kmeans",
  parameters = "____",
  columns_w_values_for_single_analyte = c("iron_percent", "density_g_per_cm3"),
  columns_w_sample_ID_info = c("fragment", "find_site")
)

ggplot(meteorite_elbow, aes(x = k, y = within_cluster_variance)) +
  geom_line() +
  geom_point(size = 3) +
  theme_classic()

Self-check: the within-cluster variance is lowest at k = 10. Why is 10 not the number of clusters to choose?

Because runMatrixAnalysis cannot make more than six clusters Because an even number of clusters is never correct Because the variance always falls as clusters are added, so the lowest value says nothing; the bend, where extra clusters stop helping much, is what to read

dbscan

There is another method to define clusters that we call dbscan. In this method, not all points are necessarily assigned to a cluster, and we define clusters according to a set of parameters, instead of simply defining the number of clusters. dbscan grows a cluster outwards from points that have at least “k” neighbors within a “threshold” distance of them, and a point that is not close enough to any cluster is left unassigned. In interactive mode, runMatrixAnalysis() will again load an interactive means of selecting these two parameters. In the context of a markdown document, escape room, or sandbox, simply provide “k” and “threshold” to the parameters argument. Note that the threshold is measured on the scaled columns, not in percent or g/cm³:

meteorite_dbscan <- runMatrixAnalysis(
    data = meteorite_fragments,
    analysis = c("dbscan"),
    parameters = c(4, 0.3),
    columns_w_values_for_single_analyte = c("iron_percent", "density_g_per_cm3"),
    columns_w_sample_ID_info = c("fragment", "find_site")
)

We can make the plot in the same way, but please note that to get geom_mark_ellipse to omit the ellipse for NAs it needs data without NAs:

ggplot() +
  geom_point(
    data = meteorite_dbscan,
    aes(x = iron_percent, y = density_g_per_cm3, fill = cluster), shape = 21, size = 4, alpha = 0.6
  ) +
  geom_mark_ellipse(
    data = drop_na(meteorite_dbscan),
    aes(x = iron_percent, y = density_g_per_cm3, label = cluster, fill = cluster), alpha = 0.2
  ) +
  theme_classic() +
  labs(x = "iron (%)", y = "density (g/cm³)") +
  scale_fill_manual(values = discrete_palette)

The meteorite fragments colored by dbscan cluster. The same three clumps form three clusters, but the two in-between fragments are left unassigned (NA) and have no ellipse.

dbscan finds the same three groups, but this time the two in-between fragments are left out, with NA as their cluster. For a meteorite hunter that is the more useful answer: a fragment that fits none of the known kinds is exactly the one worth a closer look.

concept check

Fill in the blank with the function that tests for missing values, so that only the fragments dbscan left unassigned are kept. Press Run.

meteorite_dbscan <- runMatrixAnalysis(
  data = meteorite_fragments,
  analysis = "dbscan",
  parameters = c(4, 0.3),
  columns_w_values_for_single_analyte = c("iron_percent", "density_g_per_cm3"),
  columns_w_sample_ID_info = c("fragment", "find_site")
)

meteorite_dbscan %>%
  filter(____(cluster))

Self-check: k-means was run on the same two fragments earlier in the chapter. What did it do with them?

It placed each of them in one of the three clusters anyway It left them unassigned, just as dbscan did It put the two of them together in a fourth cluster of their own

summarize by cluster

One more important point: when using k-means or dbscan, we can use the clusters as groupings for summary statistics. The simplest summary is a count: how many fragments ended up in each cluster? Group by the cluster column and use n(), which counts the rows in each group:

meteorite_dbscan %>%
  group_by(cluster) %>%
  summarize(n = n())
## # A tibble: 4 × 2
##   cluster       n
##   <chr>     <int>
## 1 cluster_1    14
## 2 cluster_2    12
## 3 cluster_3     9
## 4 <NA>          2

Each real cluster gets a row with its size, and the two fragments dbscan left unassigned are counted together in an NA row. A count like this is often the question itself: which kind of meteorite turned up most often, or how big is the largest group?

concept check

Using meteorite_dbscan from the previous concept check, fill in the blank with the function that counts the rows in each group. Press Run.

meteorite_dbscan %>%
  group_by(cluster) %>%
  summarize(n = ____())

Self-check: what does the row with NA in the cluster column count?

Fragments that were missing an iron or density measurement Fragments that dbscan did not assign to any cluster A fourth, smaller cluster that dbscan found

joining data

We can also ask how the two measurements we did not cluster on, nickel and silica, differ among the clusters. The clustering output only carries the ID columns and the columns we clustered on, so nickel and silica have to be brought back from meteorite_fragments. Combining two tables like this is called a join. left_join(x, y, by = ...) keeps every row of x, finds the row of y with the same value in the by column(s), and adds that row’s other columns alongside. Each fragment has its own ID in the fragment column of both tables, so we match on that, and take only the two columns we need from meteorite_fragments:

meteorite_dbscan <- left_join(
  meteorite_dbscan,
  select(meteorite_fragments, fragment, nickel_percent, silica_percent),
  by = "fragment"
)

Selecting first matters: any column that is in both tables but not in by comes back twice, renamed with .x and .y endings, and code that refers to it by its old name stops working. Naming by is worth it too: left out, left_join matches on every column the two tables share, which is not always what was intended. (The networks chapter used full_join(), which works the same way but also keeps rows of y that have no match in x.) With nickel and silica back in place, we can summarize them by cluster:

meteorite_dbscan %>%
  select(cluster, nickel_percent, silica_percent) %>%
  pivot_longer(cols = c(2,3), names_to = "measurement", values_to = "percent") %>%
  drop_na() %>%
  group_by(cluster, measurement) -> meteorite_dbscan_clean

plot_2 <- ggplot() + 
  geom_col(
    data = summarize(
      meteorite_dbscan_clean,
      mean = mean(percent), sd = sd(percent)
    ),
    aes(x = mean, y = cluster, fill = cluster),
    color = "black", alpha = 0.6
  ) +
  geom_errorbar(
    data = summarize(
      meteorite_dbscan_clean,
      mean = mean(percent), sd = sd(percent)
    ),
    aes(y = cluster, xmin = mean-sd, xmax = mean+sd),
    color = "black", width = 0.2
  ) +
  facet_grid(.~measurement, scales = "free_x") + theme_bw() +
  theme(panel.spacing = unit(1.5, "lines")) +
  geom_jitter(
    data = meteorite_dbscan_clean,
    aes(x = percent, y = cluster, fill = cluster), width = 0, height = 0.05,
    shape = 21
  ) +
  labs(x = "percent") +
  scale_fill_manual(values = discrete_palette) 

plot_1 <- ggplot() +
  geom_point(
    data = meteorite_dbscan,
    aes(x = iron_percent, y = density_g_per_cm3, fill = cluster), shape = 21, size = 4, alpha = 0.6
  ) +
  geom_mark_ellipse(
    data = drop_na(meteorite_dbscan),
    aes(x = iron_percent, y = density_g_per_cm3, label = cluster, fill = cluster), alpha = 0.2
  ) +
  theme_classic() +
  labs(x = "iron (%)", y = "density (g/cm³)") +
  scale_fill_manual(values = discrete_palette)

plot_1 + plot_2

Left: the meteorite fragments colored by dbscan cluster. Right: mean nickel and silica content in each cluster, shown as horizontal bars with the clusters on the vertical axis, with error bars showing one standard deviation and the individual fragments overlaid. Nickel rises and silica falls from the low-iron cluster to the high-iron cluster.

The clusters were found from iron and density alone, yet they also separate cleanly on nickel and silica: the iron-rich cluster carries the most nickel (iron meteorites are an iron-nickel alloy) and almost no silica, while the stony cluster is the reverse. When a grouping found from some measurements also explains measurements it never saw, that is good evidence the groups are real rather than an artifact of the method.

concept check

Fill in the blank with the column that identifies a fragment in both tables, so that each fragment’s nickel content is matched up with its cluster. Press Run.

meteorite_dbscan %>%
  select(fragment, cluster) %>%
  left_join(meteorite_fragments, by = "____") %>%
  select(fragment, cluster, nickel_percent)

Self-check: what does the by argument tell left_join?

Which column to use to match each row of the first table to a row of the second Which columns to bring across from the second table Which column to sort the joined table by

Finally, nothing limits clustering to two columns: any number of columns can be listed in columns_w_values_for_single_analyte and k-means and dbscan work exactly the same way. The only difficulty is drawing the result, since a scatter plot shows two measurements at a time. One common answer is to cluster on all of the columns and then color a principal components plot by cluster, using the PCA from the dimensional reduction chapter.

chapter summary

Download this summary as a one-page PDF

Choosing a flat clustering method. This decision flowchart starts from a set of samples to sort into groups and ends at a runMatrixAnalysis analysis and its parameters setting. Yellow boxes are questions about the samples and green boxes are recommended settings.

Flat clustering sorts samples straight into groups. K-means places every sample in the cluster whose center it sits closest to, so it needs the number of clusters determined ahead of time and can never leave a sample out of a cluster. Dbscan grows clusters outwards from points with enough close neighbors and can leave a sample that fits no cluster as NA. For k-means, plot the within-cluster variance against k and examine whether there is a bend in the plot (the elbow). A sharp elbow means the number of groups is a property of the data, and a gentle one means that the data are not obviously structured into clusters. Both methods scale the data first by default. Once samples have clusters, the cluster column works like any other grouping column.

Function or setting	What it does	What to vary
`runMatrixAnalysis(analysis = "kmeans", parameters = c(n_clusters = 3))`	sorts every sample into the number of clusters asked for	`n_clusters`
`analysis = "kmeans", parameters = "elbow"`	returns `k` and `within_cluster_variance` for one to ten clusters, for an elbow plot	`c("elbow", 6)` sets a different maximum
`analysis = "dbscan", parameters = c(4, 0.3)`	grows clusters from samples with at least k neighbors within a threshold distance, and leaves the rest as `NA`	k (the number of neighbors), then the threshold (measured on the scaled columns)
`columns_w_values_for_single_analyte`	names the columns to cluster on	any number of columns
`scale_variance`	scales each column before clustering (`TRUE` by default)	leave on when columns are in different units
`group_by(cluster) %>% summarize(n = n())`	counts the samples in each cluster (the `NA` row counts unassigned samples, not a cluster)	other summary statistics
`left_join(x, y, by = "fragment")`	brings back columns that were not clustered on	`by`, and `select()` the needed columns of `y` first
`geom_mark_ellipse(aes(fill = cluster, label = cluster))`	draws an ellipse around each cluster	use `drop_na()` on dbscan output to skip the `NA` samples

exercises

Where the Waters Divide

Flat clustering · k-means + choosing k + scaling

An order of alchemists bled a valley dry to feed their great sorting engine, and stole the stone that kept its springs sweet. Climb down through the failing waterworks, sort as its makers sorted to work the levers, and carry the Wellheart home before the flood closes over the engine floor.

Nobody Works Alone

Flat clustering · k-means, choosing k + scaling · two players

Sixty robberies in three years, and nobody can say how many thieving crews are behind them. Play as a pair: one of you works the crime scenes, the other follows the stolen goods to their resellers. Your counts will not agree, and you will have to talk each other through to the end.

Integrated Bioanalytics

flat clustering

kmeans

concept check

choosing k

concept check

dbscan

concept check

summarize by cluster

concept check

joining data

concept check

chapter summary

exercises

further reading