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.2689417

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

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.2337931

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.5694224

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

Choosing k, dissimilarity, and method

cluster_choice() sweeps any combination of k, dissimilarity, and method in a single call and returns one row per configuration with silhouette, mean within-cluster distance, and cluster-size balance. Pass a vector to any axis to sweep it; dissimilarity = "all" and method = "all" expand to every supported option.

ch <- cluster_choice(net, k = 2:4,
                      method = c("pam", "ward.D2", "complete", "average"))
ch
#> Cluster Choice (sweep: k x method)
#> 
#>  k method   silhouette within_dist sizes      ratio   best    
#>  2 pam      0.558      19.500      [108, 418]   3.870         
#>  3 pam      0.280      17.379      [94, 260]    2.766         
#>  4 pam      0.158      16.538      [94, 158]    1.681         
#>  2 ward.D2  0.588      19.650      [91, 435]    4.780         
#>  3 ward.D2  0.569      18.389      [8, 435]    54.375         
#>  4 ward.D2  0.513      17.601      [8, 435]    54.375         
#>  2 complete 0.810      22.609      [8, 518]    64.750         
#>  3 complete 0.802      22.360      [1, 518]   518.000         
#>  4 complete 0.781      22.265      [1, 518]   518.000         
#>  2 average  0.824      23.400      [4, 522]   130.500 <-- best
#>  3 average  0.816      23.153      [1, 522]   522.000         
#>  4 average  0.763      21.879      [1, 515]   515.000

plot() accepts an explicit type = so you can pick the chart that suits the sweep shape ("lines", "bars", "heatmap", "tradeoff", "facet", or the default "auto").

plot(ch, type = "lines")

To compare several dissimilarities at a fixed k, pass a vector and use the bars plot type. dissimilarity = "all" expands to every supported metric (hamming, osa, lv, dl, lcs, qgram, cosine, jaccard, jw); some metrics require denser sequences than others, so for this short-sequence sample we name a curated subset directly. abbrev = TRUE shortens the metric names on the tick labels.

ch_d <- cluster_choice(net, k = 2,
                        dissimilarity = c("hamming", "lv", "lcs"),
                        method = "ward.D2")
ch_d
#> Cluster Choice (sweep: dissimilarity)
#> 
#>  dissimilarity silhouette within_dist sizes      ratio best    
#>  hamming       0.588      19.650      [91, 435]  4.780 <-- best
#>  lv            0.584      15.375      [110, 416] 3.782         
#>  lcs           0.514      17.544      [126, 400] 3.175

plot(ch_d, type = "bars", abbrev = TRUE)

A high silhouette can come from one big cluster and a tiny outlier cluster. The tradeoff type plots silhouette against the cluster-size ratio so you can see both axes at once.

plot(ch_d, type = "tradeoff", abbrev = TRUE)

Higher silhouette indicates better-defined clusters. The <-- best marker in the printed table is the silhouette-max row; summary(ch)$best returns it programmatically. We select ward.D2 with 2 clusters here.

clust <- build_clusters(net, k = 2, method = "ward.D2", seed = 42)
summary(clust)
#> Sequence Clustering Summary
#>   Method:        ward.D2 
#>   Dissimilarity: hamming 
#>   Silhouette:    0.5877 
#> 
#> Per-cluster statistics:
#>  cluster size mean_within_dist
#>        1  435         13.14836
#>        2   91         50.72747
#>   cluster size mean_within_dist
#> 1       1  435         13.14836
#> 2       2   91         50.72747

Validating the choice with cluster_diagnostics()

Once a clustering is fit, cluster_diagnostics() returns a uniform diagnostic surface that works for both distance-based fits (net_clustering) and model-based fits (net_mmm, net_mmm_clustering). The print method shows per-cluster size, mean within-distance, and silhouette for distance-based fits; AvePP, mixing share, and per-cluster classification error for MMM fits.

diag <- cluster_diagnostics(clust)
diag
#> Cluster Diagnostics (distance) [ward.D2 / hamming]
#>   Sequences: 526  |  Clusters: 2
#>   Quality: silhouette = 0.588
#> 
#>   Cluster  N            Mean within-dist  Silhouette
#>   1        435 (82.7%)  13.148            0.727
#>   2        91 (17.3%)   50.727            -0.076

plot(diag, type = "silhouette")

as.data.frame(diag) returns the per-cluster table for downstream use.

From build_clusters() to per-cluster networks

build_clusters() returns clustering only; pass it to build_network() to get one network per cluster as a netobject_group (group networks). The shortcut cluster_network() (below) does both steps in one call.

clust <- build_clusters(net, k = 2, method = "ward.D2")
cluster_net <- build_network(clust)
cluster_net
#> Group Networks (2 clusters via ward.D2 / hamming)
#> 
#>   Group      Nodes  Edges  Weights         N
#>   Cluster 1  3      9      [0.127, 0.634]  435 (82.7%)
#>   Cluster 2  3      9      [0.109, 0.571]  91 (17.3%)

plot(cluster_net)

Shortcut: cluster_network() (distance-based)

Combines build_clusters() + build_network() into one call.

## `cograph::cluster_network()` also exists with a different signature
## (matrix aggregation); qualify with `Nestimate::` to avoid masking.
grp_dist <- Nestimate::cluster_network(net, k = 2, cluster_by = "ward.D2")
grp_dist
#> Group Networks (2 clusters via ward.D2 / hamming)
#> 
#>   Group      Nodes  Edges  Weights         N
#>   Cluster 1  3      9      [0.127, 0.634]  435 (82.7%)
#>   Cluster 2  3      9      [0.109, 0.571]  91 (17.3%)

Shortcut: cluster_mmm() (model-based)

Model-based equivalent: fits a mixture of Markov models and returns per-cluster networks.

grp_mmm <- cluster_mmm(net, k = 2)
grp_mmm
#> Group Networks (2 clusters from MMM)
#> 
#>   Group      Nodes  Edges  Weights         N
#>   Cluster 1  3      9      [0.090, 0.565]  192 (36.5%)
#>   Cluster 2  3      9      [0.144, 0.647]  334 (63.5%)

Both return a netobject_group — a list of netobjects with clustering info accessible via attr(, "clustering"):

# Access cluster assignments
attr(grp_dist, "clustering")$assignments[1:10]
#>  [1] 1 1 2 2 1 1 1 2 1 1
attr(grp_mmm, "clustering")$assignments[1:10]
#>  [1] 1 2 1 2 2 2 1 1 2 2

# Access individual cluster networks
grp_dist[[1]]$weights[1:3, 1:3]
#>               Directive Evaluative Metacognitive
#> Directive     0.6339144  0.1806656     0.1854200
#> Evaluative    0.4807407  0.3925926     0.1266667
#> Metacognitive 0.4261294  0.3406593     0.2332112

Comparing Clusters

We can compare which transition probabilities differ significantly among clusters using permutation testing:

comparison <- permutation(grp_dist, iter = 100)