flat clustering

“Do my samples fall into definable clusters?”

kmeans

One of the questions we’ve been asking is “which of my samples are most closely related?”. We’ve been answering that question using clustering. However, now that we know how to run principal components analyses, we can use another approach. This alternative approach is called k-means, and can help us decide how to assign our data into clusters. It is generally desirable to have a small number of clusters, however, this must be balanced by not having the variance within each cluster be too big. To strike this balance point, the elbow method is used. For it, we must first determine the maximum within-group variance at each possible number of clusters. An illustration of this is shown in A below:

One we know within-group variances, we find the “elbow” point - the point with minimum angle theta - thus picking the outcome with a good balance of cluster number and within-cluster variance (illustrated above in B and C.)

Let’s try k-means using runMatrixAnalysis. For this example, let’s run it on the PCA projection of the alaska lakes data set. We set analysis = "kmeans" and hand the number of clusters we want to the parameters argument. Below we ask for five; the next section shows how to arrive at that number yourself. (If you leave parameters out altogether, desktop R will open a small application that lets you choose 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, always say what you want explicitly.) Notice the argument scale_variance = FALSE. It is explained in the section on scaling below, and for now it just tells runMatrixAnalysis to cluster the PCA scores exactly as they are.

alaska_lake_data %>%
  select(-element_type) %>%
  pivot_wider(names_from = "element", values_from = "mg_per_L") -> alaska_lake_data_wide

alaska_lake_data_pca <- runMatrixAnalysis(
    data = alaska_lake_data_wide,
    analysis = c("pca"),
    columns_w_values_for_single_analyte = colnames(alaska_lake_data_wide)[3:dim(alaska_lake_data_wide)[2]],
    columns_w_sample_ID_info = c("lake", "park")
)

alaska_lake_data_pca_clusters <- runMatrixAnalysis(
    data = alaska_lake_data_pca,
    analysis = c("kmeans"),
    parameters = c(5),
    columns_w_values_for_single_analyte = c("Dim.1", "Dim.2"),
    columns_w_sample_ID_info = "sample_unique_ID",
    scale_variance = FALSE
)

alaska_lake_data_pca_clusters <- left_join(alaska_lake_data_pca_clusters, alaska_lake_data_pca)

We can plot the results and color them according to the group that kmeans 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:

alaska_lake_data_pca_clusters$cluster <- factor(alaska_lake_data_pca_clusters$cluster)
ggplot() +
  geom_point(
    data = alaska_lake_data_pca_clusters,
    aes(x = Dim.1, y = Dim.2, fill = cluster), shape = 21, size = 5, alpha = 0.6
  ) +
  geom_mark_ellipse(
    data = alaska_lake_data_pca_clusters,
    aes(x = Dim.1, y = Dim.2, label = cluster, fill = cluster), alpha = 0.2
  ) +
  theme_classic() +
  coord_cartesian(xlim = c(-7,12), ylim = c(-4,5)) +
  scale_fill_manual(values = discrete_palette)

choosing k

Above we simply asserted that five clusters was the right number. How would we know? This is what the elbow method, illustrated at the top of this chapter, is for, and we can ask runMatrixAnalysis for the numbers it needs by setting parameters = "elbow":

alaska_lake_elbow <- runMatrixAnalysis(
    data = alaska_lake_data_pca,
    analysis = c("kmeans"),
    parameters = "elbow",
    columns_w_values_for_single_analyte = c("Dim.1", "Dim.2"),
    columns_w_sample_ID_info = "sample_unique_ID",
    scale_variance = FALSE
)

alaska_lake_elbow
## # A tibble: 10 × 2
##        k within_cluster_variance
##    <int>                   <dbl>
##  1     1                  166.  
##  2     2                   63.0 
##  3     3                   28.9 
##  4     4                   13.9 
##  5     5                    6.84
##  6     6                    4.47
##  7     7                    3.25
##  8     8                    2.10
##  9     9                    1.39
## 10    10                    1.11

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 — in the limit, where every sample is its own cluster, there is no within-cluster variance left at all — so the smallest value is never the answer. What we are after is the bend:

ggplot(alaska_lake_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")

Reading left to right: the first few clusters each buy us a large reduction in within-cluster variance, and then the returns collapse. Here the line has clearly flattened by about five clusters — which is where the five we used above came from. Note that this is a judgement call, not a calculation: four would also be defensible for these data, and part of the skill is recognising 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 your choice).

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 = alaska_lake_data_pca,
    analysis = c("kmeans"),
    parameters = c("elbow", 6),
    columns_w_values_for_single_analyte = c("Dim.1", "Dim.2"),
    columns_w_sample_ID_info = "sample_unique_ID",
    scale_variance = FALSE
)

scaling

K-means works from distances, so a column measured in large numbers has more say over the clusters than a column measured in small ones. Suppose one column of a table is in grams and another is in kilograms. Without any correction, the grams column would decide nearly every cluster by itself, simply because its numbers are a thousand times larger.

For this reason runMatrixAnalysis centers each column and divides it by its standard deviation before it clusters, so that every column carries the same weight. This is the default for every analysis, and it is controlled by the scale_variance argument. For k-means, parameters = "elbow" follows the same setting, so the elbow you read always belongs to the clusters you then ask for.

Leave scale_variance alone whenever your columns are in different units, or just on very different scales. Turn it off (scale_variance = FALSE) when the columns are already comparable and their relative sizes mean something. That is the situation in the examples above: the two PCA scores are on one scale, and the first axis is larger than the second because it explains more of the variation, which is information we want the clusters to use. Scaling would stretch the second axis to the same width as the first and quietly change the answer.

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 clusteres, as in kmeans. In interactive mode, runMatrixAnalysis() will again load an interactive means of selecting parameters for defining dbscan clusters (“k”, and “threshold”). In the context of markdown document, simply provide “k” and “threshold” to the parameters argument:

alaska_lake_data_pca <- runMatrixAnalysis(
    data = alaska_lake_data_wide,
    analysis = c("pca"),
    columns_w_values_for_single_analyte = colnames(alaska_lake_data_wide)[3:dim(alaska_lake_data_wide)[2]],
    columns_w_sample_ID_info = c("lake", "park")
) 

alaska_lake_data_pca_clusters <- runMatrixAnalysis(
    data = alaska_lake_data_pca,
    analysis = c("dbscan"),
    parameters = c(4, 0.45),
    columns_w_values_for_single_analyte = c("Dim.1", "Dim.2"),
    columns_w_sample_ID_info = "sample_unique_ID"
)

alaska_lake_data_pca_clusters <- left_join(alaska_lake_data_pca_clusters, alaska_lake_data_pca)

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:

alaska_lake_data_pca_clusters$cluster <- factor(alaska_lake_data_pca_clusters$cluster)
ggplot() +
  geom_point(
    data = alaska_lake_data_pca_clusters,
    aes(x = Dim.1, y = Dim.2, fill = cluster), shape = 21, size = 5, alpha = 0.6
  ) +
  geom_mark_ellipse(
    data = drop_na(alaska_lake_data_pca_clusters),
    aes(x = Dim.1, y = Dim.2, label = cluster, fill = cluster), alpha = 0.2
  ) +
  theme_classic() +
  coord_cartesian(xlim = c(-7,12), ylim = c(-4,5)) +
  scale_fill_manual(values = discrete_palette)

summarize by cluster

One more important point: when using kmeans or dbscan, we can use the clusters as groupings for summary statistics. For example, suppose we want to see the differences in abundances of certain chemicals among the clusters:

alaska_lake_data_pca <- runMatrixAnalysis(
  data = alaska_lake_data_wide,
  analysis = c("pca"),
  columns_w_values_for_single_analyte = colnames(alaska_lake_data_wide)[3:dim(alaska_lake_data_wide)[2]],
  columns_w_sample_ID_info = c("lake", "park")
)

alaska_lake_data_pca_clusters <- runMatrixAnalysis(
  data = alaska_lake_data_pca,
  analysis = c("dbscan"),
  parameters = c(4, 0.45),
  columns_w_values_for_single_analyte = c("Dim.1", "Dim.2"),
  columns_w_sample_ID_info = "sample_unique_ID"
) 

alaska_lake_data_pca_clusters <- left_join(alaska_lake_data_pca_clusters, alaska_lake_data_pca)

alaska_lake_data_pca_clusters %>%
  select(cluster, S, Ca) %>%
  pivot_longer(cols = c(2,3), names_to = "analyte", values_to = "mg_per_L") %>%
  drop_na() %>%
  group_by(cluster, analyte) -> alaska_lake_data_pca_clusters_clean

plot_2 <- ggplot() + 
  geom_col(
    data = summarize(
      alaska_lake_data_pca_clusters_clean,
      mean = mean(mg_per_L), sd = sd(mg_per_L)
    ),
    aes(x = cluster, y = mean, fill = cluster),
    color = "black", alpha = 0.6
  ) +
  geom_errorbar(
    data = summarize(
      alaska_lake_data_pca_clusters_clean,
      mean = mean(mg_per_L), sd = sd(mg_per_L)
    ),
    aes(x = cluster, ymin = mean-sd, ymax = mean+sd, fill = cluster),
    color = "black", alpha = 0.6, width = 0.5, size = 1
  ) +
  facet_grid(.~analyte) + theme_bw() +
  geom_jitter(
    data = alaska_lake_data_pca_clusters_clean,
    aes(x = cluster, y = mg_per_L, fill = cluster), width = 0.05,
    shape = 21
  ) +
  scale_fill_manual(values = discrete_palette) 

plot_1<- ggplot() +
  geom_point(
    data = alaska_lake_data_pca_clusters,
    aes(x = Dim.1, y = Dim.2, fill = cluster), shape = 21, size = 5, alpha = 0.6
  ) +
  geom_mark_ellipse(
    data = drop_na(alaska_lake_data_pca_clusters),
    aes(x = Dim.1, y = Dim.2, label = cluster, fill = cluster), alpha = 0.2
  ) +
  theme_classic() + coord_cartesian(xlim = c(-7,12), ylim = c(-4,5)) +
  scale_fill_manual(values = discrete_palette)

plot_1 + plot_2

Integrated Bioanalytics