26  Clustering for biologists

26.1 Introduction

Now that we have a way of quantifying multivariate differences between samples as distance, we can use those measures to explore patterns in data. One application of distance metrics is to use them to organize data into clusters or groups consisting of samples that are more similar to each other than they are to other samples. In other words, the within-group distances should be smaller than the between-group distances.

The figure below illustrates this concept with an imaginary example. The clusters shown cleanly separate red, blue, and green points. For example, the average distance between the red points is smaller than the average distance between the red points and the blue points.

Figure 26.1: Hypothetical dataset containing three distinct clusters of observations in two-dimensional space.

26.1.1 Types of clustering

There are many methods for clustering observations. Many methods are hierarchical, with groups being nested within other groups. Many biological phenomena are also hierarchical in nature, so these methods are naturally quite useful. Other methods are non-hierarchical, where groups are not constrained to be within other groups.

Most hierarchical methods produce non-overlapping groups. Hierarchical methods are easier to compute in general. However, some authors (see Legendre and Legendre (2012)) consider hierarchical cluster methods to be “classification” rather than true clustering. The biological reasoning is clear: the hierarchy produced by a clustering algorithm is analogous to the hierarchies we see in phylogeny or Linnaean taxonomy.

Figure 26.2: In nonhierarchical clustering, groups are not nested within each other. In hierarchical clustering, samples are partitioned into nested groups.

Clustering methods also differ in their manner of group construction. Many methods are agglomerative, which means that groups are formed by combining samples. Other methods are divisive, which means that all samples start in a single group which is then broken up. Divisive techniques are also called partitioning techniques.

Figure 26.3: Agglomerative clustering starts with all samples completely ungrouped, then puts groups together by combining samples. Divisive clustering or partitioning starts with all samples in the same group, then divides the group until some stopping rule is reached.

Examples of each class of clustering method are shown below.

Table 26.1: Common approaches to clustering and their methods of forming groups.

Clustering approach	Group formation	Common methods
Partitioning (nonhierarchical)	Observations are assigned and reassigned among a specified number of clusters	K-means clustering
Hierarchical agglomerative	Individual observations or clusters are progressively merged into larger clusters	Ward’s method (Ward 1963); single linkage; complete linkage
Hierarchical divisive	All observations begin in one cluster, which is progressively divided into smaller clusters	Association analysis (Ludwig et al. 1988; Legendre and Legendre 2012); TWINSPAN (Hill 1979; Roleček et al. 2009)

26.2 K-means clustering

26.2.1 Introduction to K-means clustering

K-means clustering is a method for dividing a dataset with n observations, with p variables, into K categories such that each observation belongs to the cluster with the most similar mean. The “mean” for each category (i.e., group) is really a position in a p-dimensional hyperspace called a centroid. In this context, centroids are calculated as the point at the central coordinate of each dimension. The mean is most often used, but the median can be a useful alternative.

There are several algorithms for k-means clustering. One of the most common works like this:

1. Select the desired number of clusters k.
2. Make an initial guess of the mean of each cluster. I.e., Start with a K random positions in the p-dimensional hyperspace. These positions are known as centroids.
3. Calculate the Euclidean distance between each sample and each centroid.
4. Assign each sample to the category corresponding to the nearest centroid.
5. Update the centroids for each cluster by calculating the new centroid for the samples assigned to each cluster.
6. Repeat steps 3, 4, and 5 until cluster assignments no longer change.

This algorithm is not guaranteed to find the best possible solution, but it quite often finds decent solutions. Distance metrics other than the Euclidean can be used, but Euclidean is most common. Transforming, scaling, and/or standardizing variables so that they are all on the same scale can improve the results.

The figure below shows the effects of K-means clustering on a pretend dataset with 2 variables with \(K = 3\). The group assignments were determined because groups with those centroids had the smallest possible set of Euclidean distances between the points and the group centroids (right panel). Centroids are shown with a thick outline.

Figure 26.4: Illustration of K-means clustering in two-dimensional space. The left panel shows the original observations, while the right panel shows their assignment to three clusters. Large, black-outlined symbols indicate cluster centroids, and dashed lines connect individual observations to their respective centroids.

How many clusters are appropriate? That’s a good question to ask, and sometimes a tricky question to answer. There are statistical ways to assess how many groups are appropriate…but, you always need to use your biological judgement when interpreting groups. Sometimes you may be justified on biological grounds using 1 group more or less than the numerical summaries indicate.

One common way to assess the optimal K is to try several K and calculate the within-cluster sum of squared errors (WSS) for all points. The optimal value of K is often assumed to be the value at which a plot of WSS vs. K starts to bend. The example elbow plot below illustrates this method for the data above.

# generate a random dataset with 3 clusters
set.seed(1234)
x <- c(rnorm(20, 5, 2), rnorm(10, 12, 2))
y <- c(rnorm(10, 3, 2), rnorm(10, 10, 2), rnorm(10, 14, 2))
dx <- data.frame(x,y)

# perform k-means clustering with k from 2 to 8
kvec <- 2:8
nk <- length(kvec)
k.list <- vector("list", nk)
for(i in 1:nk){
  k.list[[i]] <- kmeans(dx, kvec[i])
}

# gather WSS for each k
wss <- sapply(k.list, function(x){x$tot.withinss})

# plot WSS vs. k
par(mfrow=c(1,1), mar=c(5,5,1,1))
plot(kvec, wss, type="o",
     lwd=3, pch=1.5,
     xlab="Clusters (k)",
     ylab="Total WSS")

Figure 26.5: Elbow plot showing the relationship between the number of clusters (K) and total within-cluster sum of squares (WSS) in K-means clustering. WSS decreases as the number of clusters increases, with a pronounced reduction up to approximately four clusters, suggesting that K = 4 may be a reasonable choice.

Notice that total WSS drops off much more slowly above 3 clusters. That is, going from 3 to 4 clusters reduces WSS much less than going from 2 to 3 clusters. The result above could reasonably justify 3 or 4 clusters. The figure below shows how the data are divided into 3 and 4 clusters.

# assemble values needed for plots
## get clusters
dx$k3 <- k.list[[which(kvec == 3)]]$cluster
dx$k4 <- k.list[[which(kvec == 4)]]$cluster

## coordinates of group centroids
km3 <- k.list[[which(kvec == 3)]]
km4 <- k.list[[which(kvec == 4)]]

dx$cenx3 <- km3$centers[dx$k3,1]
dx$ceny3 <- km3$centers[dx$k3,2]
dx$cenx4 <- km4$centers[dx$k4,1]
dx$ceny4 <- km4$centers[dx$k4,2]

## colors for each group
dx$col3 <- rainbow(3)[dx$k3]
dx$col4 <- rainbow(4)[dx$k4]

## symbols for each group
dx$pch3 <- (15:17)[dx$k3]
dx$pch4 <- (15:18)[dx$k4]

# make plot
par(mfrow=c(1,2), mar=c(5, 5, 2, 1))
plot(x,y, pch=dx$pch3, col=dx$col3, cex=1.4,
    xlab="Variable 1", ylab="Variable 2",
    main=expression(italic(k)==3))
segments(dx$cenx3, dx$ceny3, dx$x, dx$y,
    col=dx$col3, lwd=2, lty=2)
points(dx$cenx3, dx$ceny3, pch=21, cex=2,
    bg="yellow", col=dx$col3, lwd=4)

plot(x,y, pch=dx$pch4, col=dx$col4, cex=1.4,
    xlab="Variable 1", ylab="Variable 2",
    main=expression(italic(k)==4))
segments(dx$cenx4, dx$ceny4, dx$x, dx$y,
    col=dx$col4, lwd=2, lty=2)
points(dx$cenx4, dx$ceny4, pch=21, cex=2,
    bg="yellow", col=dx$col4, lwd=4)

Figure 26.6: Comparison of K-means clustering solutions with three clusters (left) and four clusters (right) for the same hypothetical dataset. Increasing K from 3 to 4 divides the upper-right group into two smaller clusters, illustrating how the number of clusters influences group assignments and centroid locations.

Most of the group assignments stayed the same, but the algorithm split the blue group when K increased from 3 to 4. Accepting the 3 cluster or 4 cluster solution is as much of a biological decision as it is a statistical one. When you get to this point, you need to think about whether the groups make sense or not. Notice that K-means clustering does not necessarily result in groups of equal sizes. This may or may not be a problem depending on the structure of your dataset and the question you are trying to answer.

26.2.2 Example with NEON data

Let’s revisit the NEON small mammal data. Download the individual measurement dataset and the taxonomic dataset and put them in your working directory.

tax <- read.csv("data-neon-taxa.csv")
ind <- read.csv("data-neon-ind-meas.csv")

We are going to use these data to try out a few common clustering techniques. First we need to organize the data. Cut the measurement data down to just the columns we need: individual, species, and the four standard measurements.

use.cols <- c("tagID", "abb", "totalLength",
              "tailLength", "hindfootLength", "earLength" )
dx <- ind[,use.cols]

# only want individuals with no missing measurements
dx <- dx[complete.cases(dx),]

# pick Peromyscus maniculatus (Perman) and
# P. leucopus (Perleu). These species are
# notoriously difficult to tell apart in
# the field. Let's see if we can do it with clustering.
use.spp <- c("Perleu", "Perman")
dx <- dx[dx$abb %in% use.spp,]

# assign colors and shapes for plotting later
cols <- c("tomato", "royalblue")
pchs <- 1:2
names(cols) <- use.spp
names(pchs) <- use.spp
dx$col <- cols[dx$abb]
dx$pch <- pchs[dx$abb]

The Eastern Deer Mouse (Peromyscus maniculatus) and the White-Footed Mouse (P. leucopus) are two morphologically similar species of mice that are widespread throughout much of North America. Among mammalogists, they are notoriously difficult to tell apart in the field. This is a problem, because they are both habitat generalists and co-occur in many places (Light et al. (2021)).

Figure 26.7: The Eastern Deer Mouse (Peromyscus maniculatus) and the White-Footed Mouse (P. leucopus). Both photos by the author.

Several field characters or suites of characters have been shown to differentiate them, with some success (Rich et al. (1996)), but regional variation in both species confounds these methods. Some studies have found that over a third of species in museum collections may have been misidentified (Light et al. (2021)). For this example analysis, we are going to assume that all individuals in the dataset have been identified correctly.

Let’s apply K-means clustering to see if we can identify any groups of similar individuals. If the species can be distinguished morphologically, then the clusters may recover the species. We expect to find 2 clusters, but there may be more. We will test \(K\) from 2 through 8.

Unlike most of the multivariate methods in this guide, we won’t need the package vegan: kmeans() is a base function.

# numbers of groups to test
kvec <- 2:8

# set up a list to catch results
nk <- length(kvec)
k.list <- vector("list", nk)

# columns with measurements
meas.cols <- c("totalLength", "tailLength",
               "hindfootLength","earLength")

# fit clusters in a group
for(i in 1:nk){
    k.list[[i]] <- kmeans(dx[,meas.cols], kvec[i])
}

# gather within-group sums of squares (wss)
# from each K for scree plot
wss <- sapply(k.list, function(x){x$tot.withinss})

# plot wss vs. K
par(mar=c(5,5,1,1))
plot(kvec, wss, type="o", lwd=3, pch=1.5,
    xlab="Clusters (k)",
    ylab="Total WSS")

Figure 26.8: Elbow plot showing the relationship between the number of clusters (K) and total within-cluster sum of squares (WSS) for K-means clustering of morphological measurements. WSS decreases as K increases, but the absence of a clearly defined elbow makes it difficult to identify an optimal number of clusters.

Unlike our first example (Figure 26.5), this curve doesn’t have a clear break point. It is definitely leveling out by 4 to 6 clusters, but there isn’t as clean of an answer as the plot for the data that I made up earlier. Let’s explore the potential clusters in an ordination plot. This is a way to visualize multivariate data that we will cover in Chapter 27. For now, just understand that points that are closer together are more similar to each other.

library(vegan)

# quick pca on scaled and centered data
pca.dx <- scale(dx[,meas.cols])
p1 <- pca(dx[,meas.cols])
summary(p1)

Call:
pca(X = dx[, meas.cols]) 

Partitioning of variance:
              Inertia Proportion
Total           459.1          1
Unconstrained   459.1          1

Eigenvalues, and their contribution to the variance 

Importance of components:
                           PC1      PC2      PC3      PC4
Eigenvalue            414.9543 38.29408 4.158315 1.693792
Proportion Explained    0.9038  0.08341 0.009058 0.003689
Cumulative Proportion   0.9038  0.98725 0.996311 1.000000

First, plot the PCA and symbolize the points by species.

xy <- scores(p1, 1:4, "sites")

par(mar=c(5,5,1,1))
plot(xy[,1], xy[,2], pch=dx$pch, col=dx$col, cex=1.4,
     xlab="PC1 (90.4%)", ylab="PC2 (8.3%)")
legend("topright", legend=use.spp,
       col=cols, pch=pchs, cex=1.4, bty="n")

Figure 26.9: Principal components analysis (PCA) of four morphological measurements from Peromyscus leucopus and P. maniculatus collected at NEON sites. The first two principal components explain 98.7% of the total variation. Although the species differ somewhat along PC1, their distributions overlap substantially, indicating that the morphological measurements do not clearly separate the two species.

There is some separation between the species along PC1, with P. leucopus having larger PC1 values and P. maniculatus having smaller scores on that axis. There is still plenty of overlap though, which makes sense given the difficulty in field identification. We also can’t rule out that some P. leucopus were misidentified as P. maniculatus and vice versa (or, for that matter, that individuals of other species were misidentified as one of P. leucopus or P. maniculatus).

Let’s compare the “known” species to the the assignments under 4, 5, and 6 cluster scenarios.

# extract clusters
dx$k4 <- k.list[[3]]$cluster # [[1]] == k = 2
dx$k5 <- k.list[[4]]$cluster
dx$k6 <- k.list[[5]]$cluster

# assign symbology
dx$pch4 <- (0:3)[dx$k4]
dx$pch5 <- (0:4)[dx$k5]
dx$pch6 <- (0:5)[dx$k6]

dx$col4 <- rainbow(4)[dx$k4]
dx$col5 <- rainbow(5)[dx$k5]
dx$col6 <- rainbow(6)[dx$k6]

# get measurement scores
zz <- scores(p1, display="species", scaling=2)
zz <- 0.1*zz

par(mar=c(5,5,1,1), mfrow=c(2,2), xpd=NA)
# plot with species
plot(xy[,1], xy[,2], pch=dx$pch, col=dx$col, cex=1.4,
     xlab="PC1 (90.4%)", ylab="PC2 (8.3%)")
legend("topright", legend=use.spp,
       col=cols, pch=pchs, cex=1.4, bty="n")
sc <- 1.03
segments(0, 0, zz[,1], zz[,2], lwd=2)
text(sc*zz[,1], sc*zz[,2], rownames(zz), adj=0)

# plot with 4 clusters
plot(xy[,1], xy[,2], pch=dx$pch4, col=dx$col4, cex=1.4,
     xlab="PC1 (90.4%)", ylab="PC2 (8.3%)")
legend("topright", legend=paste("Group", 1:4),
       col=rainbow(4), pch=0:3, bty="n")
segments(0, 0, zz[,1], zz[,2], lwd=2)
text(sc*zz[,1], sc*zz[,2], rownames(zz), adj=0)

# plot with 5 clusters
plot(xy[,1], xy[,2], pch=dx$pch5, col=dx$col5, cex=1.4,
     xlab="PC1 (90.4%)", ylab="PC2 (8.3%)")
legend("topright", legend=paste("Group", 1:5),
       col=rainbow(5), pch=0:4, bty="n")
segments(0, 0, zz[,1], zz[,2], lwd=2)
text(sc*zz[,1], sc*zz[,2], rownames(zz), adj=0)

# plot with 6 clusters
plot(xy[,1], xy[,2], pch=dx$pch6, col=dx$col6, cex=1.4,
     xlab="PC1 (90.4%)", ylab="PC2 (8.3%)")
legend("topright", legend=paste("Group", 1:6),
       col=rainbow(6), pch=0:5, bty="n")
segments(0, 0, zz[,1], zz[,2], lwd=2)
text(sc*zz[,1], sc*zz[,2], rownames(zz), adj=0)

Figure 26.10: Principal components analysis (PCA) of four morphological measurements from Peromyscus leucopus and P. maniculatus, showing species identities (upper left) and K-means clustering solutions with four (upper right), five (lower left), and six (lower right) clusters. Arrows indicate the contributions of the four morphological measurements to the principal components. Increasing the number of clusters progressively subdivides observations along PC1, which explains 90.4% of the total variation, rather than recovering the two species as distinct groups.

The PCA let’s us see something important about the clusters. First, all of the biplot vectors (black lines) go off to the right. This tells us that PC1 is basically a gradient of overall size. That’s not inherently a problem; overall size could be a factor that helps differentiate two species. Second, the clusters that the K-means algorithm discovered are essentially splitting the individuals up along that gradient of size. The top left panel of Figure 26.10 shows that P. leucopus and P. maniculatus do separate along PC1 somewhat, but they also separate somewhat along PC2.

Just for fun, let’s redraw the figure to show PC1 and PC3 instead of PC1 and PC2. This might allow us to see some better separation.

# get measurement scores
zz <- scores(p1, choices=c(1,3), display="species", scaling=2)
zz <- 0.1*zz

par(mar=c(5,5,1,1), mfrow=c(2,2), xpd=NA)
# plot with species
plot(xy[,1], xy[,3], pch=dx$pch, col=dx$col, cex=1.4,
     xlab="PC1 (90.4%)", ylab="PC3 (4.2%)")
legend("topright", legend=use.spp,
       col=cols, pch=pchs, cex=1.4, bty="n")
sc <- 1.03
segments(0, 0, zz[,1], zz[,2], lwd=2)
text(sc*zz[,1], sc*zz[,2], rownames(zz), adj=0)

# plot with 4 clusters
plot(xy[,1], xy[,3], pch=dx$pch4, col=dx$col4, cex=1.4,
     xlab="PC1 (90.4%)", ylab="PC3 (4.2%)")
legend("topright", legend=paste("Group", 1:4),
       col=rainbow(4), pch=0:3, bty="n")
segments(0, 0, zz[,1], zz[,2], lwd=2)
text(sc*zz[,1], sc*zz[,2], rownames(zz), adj=0)

# plot with 5 clusters
plot(xy[,1], xy[,3], pch=dx$pch5, col=dx$col5, cex=1.4,
     xlab="PC1 (90.4%)", ylab="PC3 (4.2%)")
legend("topright", legend=paste("Group", 1:5),
       col=rainbow(5), pch=0:4, bty="n")
segments(0, 0, zz[,1], zz[,2], lwd=2)
text(sc*zz[,1], sc*zz[,2], rownames(zz), adj=0)

# plot with 6 clusters
plot(xy[,1], xy[,3], pch=dx$pch6, col=dx$col6, cex=1.4,
     xlab="PC1 (90.4%)", ylab="PC3 (4.2%)")
legend("topright", legend=paste("Group", 1:6),
       col=rainbow(6), pch=0:5, bty="n")
segments(0, 0, zz[,1], zz[,2], lwd=2)
text(sc*zz[,1], sc*zz[,2], rownames(zz), adj=0)

Figure 26.11: Principal components analysis (PCA) of four morphological measurements from Peromyscus leucopus and P. maniculatus, showing species identities (upper left) and K-means clustering solutions with four (upper right), five (lower left), and six (lower right) clusters. PC1 explains 90.4% of the total variation, while PC3 explains 4.2%. Increasing the number of clusters progressively subdivides observations along PC1, with little apparent separation along PC3. The clustering primarily reflects variation in overall body size rather than species identity.

Nope, PC1 is still dominating the figure. Let’s check PC2 and PC3:

# get measurement scores
zz <- scores(p1, choices=c(2,3), display="species", scaling=2)
zz <- 0.1*zz

par(mar=c(5,5,1,1), mfrow=c(2,2), xpd=NA)
# plot with species
plot(xy[,2], xy[,3], pch=dx$pch, col=dx$col, cex=1.4,
     xlab="PC2 (8.3%)", ylab="PC3 (4.2%)")
legend("topright", legend=use.spp,
       col=cols, pch=pchs, cex=1.4, bty="n")
sc <- 1.03
segments(0, 0, zz[,1], zz[,2], lwd=2)
text(sc*zz[,1], sc*zz[,2], rownames(zz), adj=0)

# plot with 4 clusters
plot(xy[,2], xy[,3], pch=dx$pch4, col=dx$col4, cex=1.4,
     xlab="PC2 (8.3%)", ylab="PC3 (4.2%)")
legend("topright", legend=paste("Group", 1:4),
       col=rainbow(4), pch=0:3, bty="n")
segments(0, 0, zz[,1], zz[,2], lwd=2)
text(sc*zz[,1], sc*zz[,2], rownames(zz), adj=0)

# plot with 5 clusters
plot(xy[,2], xy[,3], pch=dx$pch5, col=dx$col5, cex=1.4,
     xlab="PC2 (8.3%)", ylab="PC3 (4.2%)")
legend("topright", legend=paste("Group", 1:5),
       col=rainbow(5), pch=0:4, bty="n")
segments(0, 0, zz[,1], zz[,2], lwd=2)
text(sc*zz[,1], sc*zz[,2], rownames(zz), adj=0)

# plot with 6 clusters
plot(xy[,2], xy[,3], pch=dx$pch6, col=dx$col6, cex=1.4,
     xlab="PC2 (8.3%)", ylab="PC3 (4.2%)")
legend("topright", legend=paste("Group", 1:6),
       col=rainbow(6), pch=0:5, bty="n")
segments(0, 0, zz[,1], zz[,2], lwd=2)
text(sc*zz[,1], sc*zz[,2], rownames(zz), adj=0)

Figure 26.12: Principal components analysis (PCA) of four morphological measurements from Peromyscus leucopus and P. maniculatus, showing species identities (upper left) and K-means clustering solutions with four (upper right), five (lower left), and six (lower right) clusters. PC2 explains 8.3% of the total variation, while PC3 explains 4.2%. Neither species nor the K-means clusters show clear separation along these axes, indicating that the clustering primarily reflects variation along PC1 rather than differences in the secondary components of morphological variation.

That is somehow worse. Okay, let’s check the 2 cluster case. We can omit the biplot vectors because they are essentially redundant with PC1.

dx$k2 <- k.list[[1]]$cluster
table(dx$abb, dx$k2, useNA="ifany")

        
           1   2
  Perleu 389 248
  Perman 160 548

dx$col2 <- rainbow(2)[dx$k2]
dx$pch2 <- (0:1)[dx$k2]

par(mar=c(5,5,1,1), mfrow=c(2,2))
# plot with species
plot(xy[,1], xy[,2], pch=dx$pch, col=dx$col, cex=1.4,
     xlab="PC1 (90.4%)", ylab="PC2 (8.3%)")
legend("topright", legend=use.spp,
       col=cols, pch=pchs, cex=1.4, bty="n")

# plot with 2 clusters vs. PC1 and PC2
plot(xy[,1], xy[,2], pch=dx$pch2, col=dx$col2, cex=1.4,
     xlab="PC1 (90.4%)", ylab="PC2 (8.3%)")
legend("topright", legend=paste("Group", 1:2),
       col=rainbow(2), pch=0:1, bty="n")

# plot with 2 clusters vs. PC1 and PC3
plot(xy[,1], xy[,3], pch=dx$pch2, col=dx$col2, cex=1.4,
     xlab="PC1 (90.4%)", ylab="PC3 (4.2%)")
legend("topright", legend=paste("Group", 1:2),
       col=rainbow(2), pch=0:1, bty="n")

# plot with 2 clusters vs. PC2 and PC3
plot(xy[,2], xy[,3], pch=dx$pch2, col=dx$col2, cex=1.4,
     xlab="PC2 (8.3%)", ylab="PC3 (4.2%)")
legend("topright", legend=paste("Group", 1:2),
       col=rainbow(2), pch=0:1, bty="n")

Figure 26.13: Principal components analysis (PCA) of four morphological measurements from Peromyscus leucopus and P. maniculatus, comparing species identities (upper left) with a K-means clustering solution containing two clusters (upper right and lower panels). The two clusters are separated primarily along PC1, which explains 90.4% of the total variation, with little separation along PC2 or PC3. Although the clustering partially distinguishes the two species, it primarily reflects variation in overall body size rather than species identity.

So what did we learn? First, the clustering is overwhelmingly driven by overall size, and the clustering mirrored that. This makes sense mathematically. K-means minimizes within-cluster sums of squared Euclidean distances. When nearly all the variation lies along one dimension, the algorithm can reduce that sum of squares most efficiently by dividing observations along that dimension. In other words, the algorithm worked exactly as intended, but it answers a different question than the one we hoped it would.

Second, the clustering was not entirely useless. The contingency table showed that the 2 cluster solution was correct about 69.7% of the time (937 out of 1,345 individuals correctly assigned). However, the accuracy is asymmetric: 77.4% of P. maniculatus are in Cluster 2, but only 61.1% of P. leucopus are in Cluster 1. So there is a real association between the clusters and species, even though the clusters don’t correspond particularly well to species boundaries. A \(\chi^2\) test on the frequency table will confirm this lack of independence of cluster and species.

Third, the analysis detected considerable variation in body size that must be associated with age, sex, geography, or some factor (or combination of factors) other than species identity. Indeed, both species exhibit considerable geographic variation, which we will explore in Chapter 27.

This leads us to an important question: if we already know the species identities, would a supervised classification method have been more appropriate? We will explore one such method, linear discriminant analysis, in the next chapter.

26.3 Hierarchical clustering

26.3.1 Introduction

Hierarchical clustering works by constructing groups of samples, and then combining samples and groups. Most hierarchical clustering algorithms are agglomerative. The general method is as follows.

1. Initialize the clustering. Each observation begins as the sole member of a cluster with one member.
2. Calculate distances. Calculate a distance matrix describing the dissimilarity between every pair of observations. The choice of distance measure (or dissimilarity, for reviewer 2) is important, so think carefully about that.
3. Find the closest clusters. Identify the two clusters separated by the smallest distance. Initially, each cluster contains only one observation.
4. Merge the clusters. Combine the two closest clusters into a single, larger cluster.
5. Update the distances. Calculate the distance between the newly formed cluster and each of the remaining clusters (Figure 26.14).
6. Repeat. Continue merging clusters and updating distances until all observations belong to a single cluster.

The important question is how to calculate distances between clusters containing multiple observations. Several methods are available, and different methods can produce different clustering results. Legendre and Legendre (2012) note that many hierarchical agglomerative clustering methods are actually special cases of a general method (Lance and Williams (1966), Lance and Williams (1967)). Rather than recalculating distances from all the individual observations after each merger, this method uses distances already calculated in previous steps.

In Figure 26.14, clusters b and c have just been merged into a new cluster, bc. We already know the distances between clusters a and b, between a and c, and between b and c. The Lance–Williams equation uses these distances, along with the numbers of observations in each cluster, to calculate the new distance \(d\left(\mathbf{a},\mathbf{bc}\right)\).

Figure 26.14: In combinatorial hierarchical agglomerative clustering methods, the distances between already-existing clusters, along with their sample sizes, are used to determine the next cluster. Here cluster bc is already formed from combining clusters b and c. Cluster a is being considered for addition to cluster bc to form cluster abc. The distance between a and bc, \(d\left(\mathbf{a},\mathbf{bc}\right)\), is a function of the distances between a and b, a and c, b and c, and their respective sample sizes.

Different clustering methods use different rules for this calculation. For example, single linkage, also known as nearest neighbor clustering, defines the distance between clusters using their closest pair of observations, whereas Ward’s method chooses mergers that minimize the increase in within-cluster variation.

26.3.2 Chaining

One potential pitfall in hierarchical clustering is a phenomenon called chaining. This is a situation where new samples are added to a pre-existing cluster one at a time, producing a “chain” of mostly pairwise connections going down one side of the dendrogram. Strong chaining in a dendrogram may indicate a continuous gradient, not a clear partition between groups. Chaining can also arise from outliers, noise, or observations that bridge otherwise distinct groups. Single-linkage agglomeration is the most vulnerable to chaining; Ward’s method and complete linkage are much less susceptible.

Figure 26.15: Hypothetical hierarchical clustering dendrogram illustrating the difference between relatively balanced branching and possible chaining. Samples B–E on the left form two pairs that subsequently merge into a larger cluster. On the right, samples G–M exhibit a more sequential pattern of mergers, with individual observations repeatedly joining an existing cluster. This pattern, known as chaining, may reflect continuous variation rather than clearly separated groups.

Figure 26.15 shows a dendrogram of a hypothetical clustering result with evidence of possible chaining. Note that on the left side (samples B, C, D, and E), the samples are partitioned into groups with both sides of most splits having a subsequent split. On the right (samples G, H, I, J, K, L, and M), each sample is joined to the tree one at a time. This indicates that variation among these samples may be driven by a continuous gradient rather than some clear grouping structure.

26.3.3 Hierarchical clustering methods

26.3.3.1 Hierarchical clustering with single linkage

One of the simplest hierarchical clustering methods is to combine samples based on the shortest pairwise dissimilarities. This is called single linkage, or sometimes nearest neighbor clustering. “Nearest neighbor” can refer to other methods as well, so we’ll refer to this clustering method as “single linkage”. The basic idea is to merge the two clusters separated by the smallest distance between any pair of observations, one from each cluster. Because the distance between clusters depends only on their closest observations, individual samples can be added to existing clusters one at a time, producing the phenomenon of chaining.

26.3.3.1.1 Example using NEON data

Let’s revisit the NEON small mammal data. Download the site-level capture dataset and the site metadata dataset and put them in your working directory.

cap <- read.csv("data-neon-captures-site.csv")
met <- read.csv("data-neon-metadata-site.csv")

Let’s use hierarchical clustering to look for structure among the sites in terms of their small mammal fauna. We will convert the abundances of each species to presence-absence (1 and 0) and check for empty sites.

# organize data frame
dx <- cap
rownames(dx) <- dx$site
dx$X <- NULL
dx$site <- NULL

# convert capture numbers to presence/absence
dx[dx>0] <- 1

# verify that all sites have at least 1 species
# should be TRUE
all(rowSums(dx) > 0)

[1] TRUE

Data frame dx has 46 rows (sites) and 134 variables (species presence). Let’s try a single-linkage clustering to look for structure among the sites.

library(vegan)

# use Bray-Curtis (Sørensen) dissimilarity
D <- vegdist(dx)

# execute clustering
clust1 <- hclust(D, method="single")

Now plot the dendrogram

plot(clust1)

Figure 26.16: Hierarchical agglomerative clustering of 46 NEON terrestrial field sites using single linkage (nearest-neighbor clustering) and Sørensen dissimilarity calculated from species presence/absence data. The height of each merger represents the dissimilarity between the two clusters being joined. Several geographically related sites cluster together at relatively low dissimilarities, while the dendrogram also exhibits chaining, with individual sites or small groups joining larger clusters sequentially.

We can superimpose that tree on the actual geography using package terra. Grab the Natural Earth Data 110 m countries dataset and put it in your working directory. Them import it using terra.

library(terra)

world <- vect("ne_110m_admin_0_countries.shp")

# crop to North America
na <- crop(world, ext(-160, -60, 15, 75))

# Create spatial points
coords <- met[,c("lat", "lon")]
pts <- vect(coords, geom = c("lon", "lat"), crs = "EPSG:4326")

# Project to a suitable planar CRS before computing centroids
pts.pr <- project(pts, "EPSG:5070")  # CONUS Albers; change if needed
na.proj <- project(na, "EPSG:5070")

xy <- crds(pts.pr)

# Initialize coordinates for all nodes
n <- nrow(xy)
nodes <- matrix(NA_real_, nrow = n - 1, ncol = 2)

# Plot sampling locations
par(mar=c(5,5,1,1))
plot(na.proj)
plot(pts.pr, pch = 19, add=TRUE)

# Draw connections for each merger
for (i in seq_len(n - 1)) {

  children <- clust1$merge[i, ]

  # Get coordinates for each child
  get_xy <- function(j) {
    if (j < 0) xy[-j, ] else nodes[j, ]
  }

  a <- get_xy(children[1])
  b <- get_xy(children[2])

  # New node at midpoint of child nodes
  nodes[i, ] <- (a + b) / 2

  # Draw branches
  segments(a[1], a[2], nodes[i, 1], nodes[i, 2], col = "gray50")
  segments(b[1], b[2], nodes[i, 1], nodes[i, 2], col = "gray50")
}

points(xy, pch = 19)

Figure 26.17: Geographic representation of a single linkage tree connecting 46 NEON terrestrial field sites based on Sørensen dissimilarity in species presence/absence. Points indicate sampling locations, and lines connect pairs of sites selected by the single linkage tree algorithm. The connections represent similarities in species composition rather than geographic proximity, although many nearby sites are linked, suggesting a geographic component to community similarity.

The results are pretty interesting. For the most part, nearby sites are linked directly together. That make sense, because sites that are close together should have more similar environmental conditions, and overlap the same species distributions.

There are other ways to explore the results. A dendrogram simplifies a complex dissimilarity matrix into a hierarchical structure. One way to assess how much information is preserved is the cophenetic correlation coefficient.

The cophenetic distance between two sites is the height at which they first join the same cluster in the dendrogram. We can check the correlation between that distance and the original Sørensen dissimilarity matrix. Values closer to 1 indicate that the dendrogram preserves the original dissimilarity relationships relatively well.

# D is your original Sørensen dissimilarity matrix
# clust1 is your hclust result

cd <- cophenetic(clust1)

cor(D, cd)

[1] 0.7754411

A cophenetic r of 0.775 is pretty good. We can visualize this on a scatterplot.

plot(as.vector(D), as.vector(cd),
     xlab = "Sørensen dissimilarity",
     ylab = "Cophenetic distance")
abline(lm(as.vector(cd) ~ as.vector(D)), col = "red")

Figure 26.18: Relationship between original Sørensen dissimilarities and cophenetic distances from single-linkage hierarchical clustering of 46 NEON terrestrial field sites. Each point represents a pair of sites, with the cophenetic distance indicating the height at which the two sites first join the same cluster in the dendrogram. The red line shows the fitted linear regression. The positive relationship indicates that the dendrogram preserves some of the original dissimilarity structure, although considerable scatter suggests imperfect representation.

This is particularly useful for comparing different hierarchical clustering methods. Just keep in mind that cophenetic correlation measures fidelity to the original dissimilarities, not whether the resulting groups are biologically meaningful.

Another way to interpret the clustering result is to identify clusters by cutting the dendrogram. The dendrogram represents a hierarchy, but you can choose a particular level to define discrete groups.

For example, suppose you want four clusters:

groups <- cutree(clust1, k = 4)

plot(clust1)
rect.hclust(clust1, k = 4, border = 2:5)

Figure 26.19: Single-linkage hierarchical clustering of 46 NEON terrestrial field sites based on Sørensen dissimilarity in species presence/absence. Colored rectangles indicate four groups identified by cutting the dendrogram at approximately 0.64. Three relatively small clusters contain sites in Puerto Rico, the Pacific coastal mountains, and the Rocky Mountains, while the fourth contains the remaining 36 sites. The highly unequal cluster sizes illustrate how single-linkage clustering can produce large, heterogeneous groups.

Four groups puts sites GUAN and LAJA in one group; SJER, SOAP, TEAK, ABBY, and WREF in another; YELL, NIWO, and RMNP in a third, and every other site in the last group. Geographically this makes sense: GUAN and LAJA are both in southwestern Puerto Rico. The second group consists of sites in the Pacific coastal mountain regions of California and Washington. The third consists of mountainous sites in the Rocky Mountains of Wyoming and Colorado. The lumping of the other 36 sites is a bit concerning, though. That group includes sites in Alaska and throughout the continental US east of the Rockies. Let’s see what happens if we allow for 7 groups instead of 4.

groups <- cutree(clust1, k = 7)

plot(clust1)
rect.hclust(clust1, k = 7, border = 2:8)

Figure 26.20: Single-linkage hierarchical clustering of 46 NEON terrestrial field sites based on Sørensen dissimilarity in species presence/absence. Colored rectangles identify six groups obtained by cutting the dendrogram at approximately 0.60. Compared with the four-group solution, the Pacific coastal mountain sites are divided into two clusters, and the Arctic site BARR forms its own cluster. However, most sites remain together in one large group, illustrating the tendency of single-linkage clustering to produce unevenly sized clusters.

This time we still have a “most of the continental US” group, but a few were peeled off into their own groups. Now BARR (Utqiaġvik, Alaska; arctic tundra) is in its own group. JORN (near Las Cruces, New Mexico) and SRER (near Tucson, Arizona) ar both in the arid southwest. Let’s try 12 groups.

groups <- cutree(clust1, k = 12)

plot(clust1)
rect.hclust(clust1, k = 12, border = 2:13)

Figure 26.21: Single-linkage hierarchical clustering of 46 NEON terrestrial field sites based on Sørensen dissimilarity in species presence/absence. Colored rectangles identify 12 groups obtained by cutting the dendrogram at approximately 0.49. Compared with the four- and six-group solutions, the 12-group solution reveals more geographically and ecologically recognizable clusters, including sites in Alaska, the Rocky Mountains, the Pacific Northwest, the southwestern deserts, and the southeastern United States.

Now we’re getting somewhere! Jumping up to 12 groups gives us a bit more separation within the “rest of the US” group. Let’s see what these groups might be:

sort(groups)

ABBY WREF BARR BART BLAN DCFS GRSM HARV MLBS NOGP ORNL SCBI SERC STEI TREE UNDE 
   1    1    2    3    3    3    3    3    3    3    3    3    3    3    3    3 
WOOD BONA DEJU HEAL TOOL CLBJ CPER KONA KONZ OAES STER UKFS DELA DSNY JERC LENO 
   3    4    4    4    4    5    5    5    5    5    5    5    6    6    6    6 
OSBS TALL GUAN LAJA JORN SRER MOAB ONAQ NIWO RMNP YELL SJER SOAP TEAK 
   6    6    7    7    8    8    9    9   10   10   10   11   11   12 

Table 26.2: Geographic and ecological characterization of 12 clusters of NEON terrestrial field sites based on Sørensen dissimilarities in species presence/absence.

Group	Sites	Description
1	ABBY, WREF	Pacific Northwest; Cascade Range forests
2	BARR	Arctic tundra of northern Alaska
3	BART, BLAN, DCFS, GRSM, HARV, MLBS, NOGP, ORNL, SCBI, SERC, STEI, TREE, UNDE, WOOD	Primarily eastern and central North American temperate forests, with some grassland and transitional sites
4	BONA, DEJU, HEAL, TOOL	Interior Alaska; boreal forests and tundra
5	CLBJ, CPER, KONA, KONZ, OAES, STER, UKFS	Central and southern Great Plains; grasslands, savannas, and transitional habitats
6	DELA, DSNY, JERC, LENO, OSBS, TALL	Southeastern United States; coastal plain, wetlands, and subtropical forests
7	GUAN, LAJA	Southwestern Puerto Rico; tropical dry forests
8	JORN, SRER	Southwestern United States; desert grasslands and shrublands
9	MOAB, ONAQ	Interior western United States; arid shrublands and desert environments
10	NIWO, RMNP, YELL	Rocky Mountains; montane and subalpine environments
11	SJER, SOAP	California Sierra Nevada; foothill and montane forests
12	TEAK	California Sierra Nevada; higher-elevation coniferous forest

Some of these descriptions are more ecologically coherent than others. A few interesting features stand out. Groups 3 and 5 encompass sites spread across substantial geographic distances, suggesting that similar communities can occur across broad regions. Groups 2 and 12 contain single sites, potentially indicating distinctive species compositions. And Groups 11 and 12 separate TEAK from SJER and SOAP, despite all three being in the Sierra Nevada, which might reflect differences in elevation or vegetation.

Let’s plot the sites on a map, symbolized by their group.

# colors and symbols for points
cols <- hcl.colors(12, palette = "Dark 3")
use.cols <- cols[groups]
pchs <- rep(c(0, 15:18, 4), times=2)
use.pch <- pchs[groups]

# Plot sampling locations
par(mar=c(5,5,1,1))
plot(na.proj)
plot(pts.pr, pch=46, add=TRUE)

# Draw connections for each merger
for (i in seq_len(n - 1)) {

  children <- clust1$merge[i, ]

  # Get coordinates for each child
  get_xy <- function(j) {
    if (j < 0) xy[-j, ] else nodes[j, ]
  }

  a <- get_xy(children[1])
  b <- get_xy(children[2])

  # New node at midpoint of child nodes
  nodes[i, ] <- (a + b) / 2

  # Draw branches
  segments(a[1], a[2], nodes[i, 1], nodes[i, 2], col = "gray50")
  segments(b[1], b[2], nodes[i, 1], nodes[i, 2], col = "gray50")
}

points(xy, col=use.cols, cex=1.5, pch = use.pch)
legend("bottomleft", legend=paste("Group", 1:12),
       pch=pchs, col=cols, cex=1.3)

Figure 26.22: Geographic distribution of 12 clusters of NEON terrestrial field sites identified using single-linkage hierarchical clustering and Sørensen dissimilarity in species presence/absence. Colors and symbols distinguish clusters, while gray lines illustrate the hierarchical relationships between sites projected onto geographic space. Many clusters correspond to recognizable geographic or ecological regions, although some geographically distant sites belong to the same cluster and nearby sites sometimes belong to different clusters.

We can see from the map that while the groups are certainly associated with geography, they are not totally determined by it. Some of the groups appear to reflect biogeographic regions rather than geographic proximity alone. The Alaskan sites, Rocky Mountain sites, southwestern deserts, and southeastern sites form fairly recognizable groups. Meanwhile, some geographically close sites belong to different clusters, presumably because their habitats or species assemblages differ.

There are two particularly interesting interpretations:

• Geographic distance-decay: Nearby sites tend to share more species, producing geographically coherent clusters. We could check on that with a scatterplot of Sørensen dissimilarity and geographic distance. Hey, why not? (Figure 26.23).
• Environmental filtering: Sites with similar environments may share species even when separated by substantial distances. Conversely, nearby sites in different habitats may have quite different assemblages.

Of course, the map alone can’t distinguish those mechanisms.

Here is the plot of Sørensen dissimilarity vs. geographic distance.

dgeo <- distance(pts)/1000
par(mar=c(5,5,1,1))
plot(dgeo, D,
     ylab="Sørensen dissimilarity",
     xlab="Geographic distance (km)",
     log="x")

Figure 26.23: Relationship between geographic distance and Sørensen dissimilarity in species composition among 46 NEON terrestrial field sites. Each point represents a pair of sites, with geographic distance (km, logarithmic scale) on the X-axis and Sørensen dissimilarity on the Y-axis. Species composition generally becomes more dissimilar as geographic distance increases, suggesting a pattern of distance decay in community similarity. However, considerable variation exists among site pairs, indicating that geographic distance alone does not explain differences in species composition.

And the correlation:

cor(D, log(dgeo))

[1] 0.7197855

The last thing I’ll point out is that the apparent strength of the geographic pattern (and its ecological sensibleness) depends on where we cut the dendrogram. With four groups, most sites fell into one enormous cluster. At twelve groups, much more recognizable geographic structure emerges. Neither solution is necessarily more correct.

26.3.3.2 Hierarchical clustering with Ward’s method

Ward’s method is another hierarchical agglomerative technique, which combines clusters to minimize the increase in within-cluster variation, measured by the sum of squared deviations from cluster means. In other words, Ward’s method merges the two clusters that will produce the smallest increase in overall within-cluster variation, while single linkage merges the two clusters containing the most similar pair of observations.

For an analogy, think of buying clothes. The single-linkage approach is to buy a new shirt that is most similar to at least one shirt you already own. The Ward’s method approach is to buy a new shirt that makes the smallest increase in the overall variation among the shirts in your wardrobe.

26.3.3.2.1 Example using NEON data

Coming soon!

26.3.3.3 Other methods for hierarchical clustering

There are other approaches to hierarchical clustering, which we will not cover here, that measure the relationships between clusters differently than single-linkage or Ward’s method. All of these are available in the R function hclust().

Some of these methods weight clusters by their number of observations (WPGMA and WPGMC), and the others are unweighted (UPGMA and UPGMC). They also differ in how updated centroid is calculated. In WPGMA and UPGMA the new centroid is the mean of the previous centroids, while in WPGMC and UPGMC in whether they use unweighted or weighted pairwise distances, and in whether they calculate centroids using arithmetic means or medians. Of these, UPGMA is probably the most commonly used.

• UPGMA: Unweighted pair-group method using arithmetic averages. The distance between two clusters is the average of all pairwise distances between observations in the two clusters, giving each observation equal weight regardless of cluster size. Use argument method="average" in hclust().
• WPGMA: Weighted pair-group method using arithmetic averages. Similar to UPGMA, but when clusters are merged, each of the two original clusters contributes equally to the new distance calculations, regardless of how many observations each contains. Use argument method="mcquitty" in hclust().
• UPGMC: Unweighted pair-group method using centroids. The distance between two clusters is calculated from their centroids (multivariate means), with each observation contributing equally to the centroid of its cluster. Use argument method="centroid" in hclust().
• WPGMC: Weighted pair-group method using centroids. Similar to UPGMC, but the centroid of a newly merged cluster is calculated by averaging the centroids of the two original clusters, giving each cluster equal weight regardless of its size. Use argument method="median" in hclust().

Hill MO. 1979. A FORTRAN program for arranging multivariate data in an ordered two-way table by classification of the individuals and attributes. TWINSPAN.

Lance GN, Williams WT. 1967. A general theory of classificatory sorting strategies: 1. Hierarchical systems. Comput. J. 9(4):373–380.

Lance G, Williams W. 1966. A generalized sorting strategy for computer classifications. Nature 212(5058):218–218.

Legendre P, Legendre L. 2012. Numerical ecology. Amsterdam: Elsevier.

Light JE, Siciliano-Martina L, Dohlanik EG, Hafner DJ, Lawing AM, Greenbaum IF. 2021. Morphological differentiation of peromyscus leucopus and p. Maniculatus in east texas. Therya 12(2):369–387.

Ludwig JA, Reynolds JF, QUARTET L, Reynolds J. 1988. Statistical ecology: a primer in methods and computing. John Wiley & Sons.

Rich SM, Kilpatrick CW, Shippee JL, Crowell KL. 1996. Morphological differentiation and identification of peromyscus leucopus and p. Maniculatus in northeastern north america. Journal of Mammalogy 77(4):985–991.

Roleček J, Tichý L, Zelený D, Chytrý M. 2009. Modified TWINSPAN classification in which the hierarchy respects cluster heterogeneity. Journal of Vegetation Science 20(4):596–602. https://doi.org/10.1111/j.1654-1103.2009.01062.x

Ward JH. 1963. Hierarchical grouping to optimize an objective function. Journal of the American Statistical Association 58(301):236–244. https://doi.org/10.1080/01621459.1963.10500845