Sequence Clustering

Visualizing Clusters

The silhouette plot shows how well each sequence fits its assigned cluster. Values near 1 indicate good fit; values near 0 suggest the sequence is between clusters; negative values indicate possible misclassification.

plot(clust, type = "silhouette")

The MDS (multidimensional scaling) plot projects the distance matrix to 2D, showing cluster separation.

plot(clust, type = "mds")

Distance Metrics

A distance metric defines how (dis)similarity between sequences is measured. In other words, it quantifies how different two sequences are from each other. Nestimate currently supports 9 distance metrics for comparing sequences:

Different metrics may produce different clustering results. You need to choose this based on your research question:

• Hamming: compares sequences position by position (best for aligned sequences of equal length).
• Edit distances (lv, osa, dl): allow insertions and deletions (best when sequences may be shifted or vary in length).
• LCS: captures shared subsequences (best when overall patterns matter more than exact alignment).

We can specify which distance metric we want to use through the dissimilarity argument:

# Levenshtein distance (allows insertions/deletions)
clust_lv <- build_clusters(net, k = 3, dissimilarity = "lv")
clust_lv$silhouette
#> [1] 0.5277661

# Longest common subsequence
clust_lcs <- build_clusters(net, k = 3, dissimilarity = "lcs")
clust_lcs$silhouette
#> [1] 0.227412

Some distance metrics may take additional arguments. For example, the Hamming distance accepts temporal weighting to emphasize earlier or later positions:

# Emphasize earlier positions (higher lambda = faster decay)
clust_weighted <- build_clusters(net, 
                               k = 3,
                               dissimilarity = "hamming",
                               weighted = TRUE,
                               lambda = 0.5)
clust_weighted$silhouette
#> [1] 0.2545125

Clustering Methods

By default, Nestimate uses PAM (Partitioning Around Medoids) to form clusters, which assigns each sequence to the cluster represented by the most central sequence (the medoid). Besides PAM, Nestimate supports hierarchical clustering methods, which build clusters step by step by progressively merging similar units into a tree-like structure (a dendrogram):

• ward.D2 (“Ward’s Method, Squared Distances”): Minimizes the increase in within-cluster variance using squared distances. Typically produces compact, well-separated clusters.
• ward.D (“Ward’s Method”): An alternative implementation of Ward’s approach using a different distance formulation. Similar behavior, but results may vary slightly.
• complete (“Complete Linkage”): Defines the distance between clusters as the maximum distance between their members. Produces tight, compact clusters.
• average (“Average Linkage”): Uses the average distance between all pairs of points across clusters. Provides a balance between compactness and flexibility.
• single (“Single Linkage”): Uses the minimum distance between points in two clusters. Can capture chain-like structures but may lead to - loosely connected clusters.
• mcquitty (“McQuitty’s Method” / “WPGMA”): A weighted version of average linkage that gives equal weight to clusters regardless of size.
• centroid (“Centroid Linkage”): Defines cluster distance based on the distance between cluster centroids (means). Can produce intuitive groupings but may introduce inconsistencies in the hierarchy.

To use any of these methods instead of PAM, we need to provide the method argument to build_clusters.

# Ward's method (minimizes within-cluster variance)
clust_ward <- build_clusters(net, k = 3, method = "ward.D2")
clust_ward$silhouette
#> [1] 0.411989

# Complete linkage
clust_complete <- build_clusters(net, k = 3, method = "complete")
clust_complete$silhouette
#> [1] 0.6411682

Choosing k (Number of Clusters)

To choose the right clustering solution and method, we need to compare the silhouette scores across different k values and clustering methods (and also distance metrics if we want):

methods <- c("pam", "ward.D2", "complete", "average")

silhouettes <- lapply(methods, function(m) {
  sapply(2:4, function(k) {
    build_clusters(net, k = k, method = m, seed = 42)$silhouette
  })
})

names(silhouettes) <- methods

silhouettes
#> $pam
#> [1] 0.3802219 0.3901787 0.3375285
#> 
#> $ward.D2
#> [1] 0.8221011 0.4119890 0.4202447
#> 
#> $complete
#> [1] 0.8221011 0.6411682 0.6190995
#> 
#> $average
#> [1] 0.8221011 0.6411682 0.6190995

methods <- names(silhouettes)
colors <- rainbow(length(methods))

plot(2:4, silhouettes[[1]], type = "b", pch = 19, col = colors[1],
     xlab = "Number of clusters (k)",
     ylab = "Average silhouette width",
     ylim = c(0, 1),
     main = "Choosing k")

for (i in 2:length(methods)) {
  lines(2:4, silhouettes[[i]], type = "b", pch = 19, col = colors[i])
}

legend("topright", legend = methods, col = colors, lty = 1, pch = 19)

Higher silhouette scores indicate better-defined clusters. Look for an “elbow” or maximum. Here we select ward.D2 with 2 clusters, which yields a reasonable silhouette width.

clust <- build_clusters(net, k = 2, method = "ward.D2", seed = 42)

summary(clust)
#> Sequence Clustering Summary
#>   Method:        ward.D2 
#>   Dissimilarity: hamming 
#>   Silhouette:    0.8221 
#> 
#> Per-cluster statistics:
#>  cluster size mean_within_dist
#>        1   97          17.7049
#>        2    1           0.0000