--- title: "Sequence Clustering" author: - name: Mohammed Saqr url: https://saqr.me affiliation: University of Eastern Finland - name: Sonsoles López-Pernas url: https://sonsoles.me affiliation: University of Eastern Finland - name: Kamila Misiejuk url: https://kamilamisiejuk.com affiliation: FernUniversität in Hagen, Germany output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Sequence Clustering} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5, fig.alt = "Visualization" ) ``` ## Introduction Clusters represent subgroups within the data that share similar patterns. Such patterns may reflect similar temporal dynamics (when we are analyzing sequence data, for example) or relationships between variables (as is the case in psychological networks). Units within the same cluster are more similar to each other, while units in different clusters differ more substantially. In this vignette, we demonstrate how to perform clustering on sequence data using `Nestimate`. ## Data To illustrate clustering, we will use the `human_long` dataset, which contains 10,796 coded human interactions from 429 human-AI pair programming sessions across 34 projects. Each row represents a single interaction event with a timestamp, session identifier, and action label. ```{r data} library(cograph) library(Nestimate) data("human_long") head(human_long) ``` We can build a transition network using this dataset using `build_network`. We need to determine the `actor` (`session_id`), the `action` (`cluster`), and the `time` (`timestamp`). We will use the overall network object as the starting point to find subgroups since it structures the raw data into the appropriate units of analysis to perform clustering. ```{r} net <- build_network(human_long, method = "tna", action = "cluster", actor = "session_id", time = "timestamp") ``` ## Dissimilarity-based Clustering Dissimilarity-based clustering groups units of analysis (in our case, sessions, since that is what we provided as `actor`) by directly comparing their observed sequences. In our case, each session is represented by its sequence of actions, and similarity between sessions is defined using a distance metric that quantifies how different two sequences are. To implement this method using `Nestimate`, we can use the `build_clusters()` function, which takes either raw sequence data or a network object such as the `net` object that we estimated (which also contains the original sequences in `$data`): ```{r cluster-basic} clust <- build_clusters(net, k = 3) clust ``` The default clustering mechanism uses **Hamming distance** (number of positions where sequences differ) with **PAM** (Partitioning Around Medoids). The result contains the cluster `assignments` (which cluster each session belongs to), the cluster `sizes`, and a `silhouette` score that reflects the quality of the clustering (higher values indicate better separation between clusters), among other useful information. ```{r cluster-components} # Cluster assignments (first 20 sessions) head(clust$assignments, 20) # Cluster sizes clust$sizes # Silhouette score (clustering quality: higher is better) clust$silhouette ``` ### 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. ```{r cluster-plot, fig.alt = "Silhouette plot showing cluster quality"} plot(clust, type = "silhouette") ``` The MDS (multidimensional scaling) plot projects the distance matrix to 2D, showing cluster separation. ```{r cluster-mds, fig.alt = "MDS plot 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: | Metric | Description | Best for | |--------|-------------|----------| | `hamming` | Positions where sequences differ | Equal-length sequences | | `lv` | Levenshtein (edit distance) | Variable-length, insertions/deletions | | `osa` | Optimal string alignment | Edit distance + transpositions | | `dl` | Damerau-Levenshtein | Full edit + adjacent transpositions | | `lcs` | Longest common subsequence | Preserving order, ignoring gaps | | `qgram` | Q-gram frequency difference | Pattern-based similarity | | `cosine` | Cosine of q-gram vectors | Normalized pattern similarity | | `jaccard` | Jaccard index of q-grams | Set-based pattern overlap | | `jw` | Jaro-Winkler | Short strings, typo detection | 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: ```{r cluster-metrics} # Levenshtein distance (allows insertions/deletions) clust_lv <- build_clusters(net, k = 3, dissimilarity = "lv") clust_lv$silhouette # Longest common subsequence clust_lcs <- build_clusters(net, k = 3, dissimilarity = "lcs") clust_lcs$silhouette ``` Some distance metrics may take additional arguments. For example, the Hamming distance accepts **temporal weighting** to emphasize earlier or later positions: ```{r cluster-weighted} # Emphasize earlier positions (higher lambda = faster decay) clust_weighted <- build_clusters(net, k = 3, dissimilarity = "hamming", weighted = TRUE, lambda = 0.5) clust_weighted$silhouette ``` ### 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`. ```{r cluster-methods} # Ward's method (minimizes within-cluster variance) clust_ward <- build_clusters(net, k = 3, method = "ward.D2") clust_ward$silhouette # Complete linkage clust_complete <- build_clusters(net, k = 3, method = "complete") clust_complete$silhouette ``` ### 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. ```{r choose-k} ch <- cluster_choice(net, k = 2:4, method = c("pam", "ward.D2", "complete", "average")) ch ``` `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"`). ```{r choose-k-plot, fig.alt = "Silhouette across k for each clustering method"} 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. ```{r choose-d} ch_d <- cluster_choice(net, k = 2, dissimilarity = c("hamming", "lv", "lcs"), method = "ward.D2") ch_d ``` ```{r choose-d-plot, fig.alt = "Silhouette per dissimilarity at k = 2"} 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. ```{r choose-tradeoff, fig.alt = "Quality vs cluster-size balance"} 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. ```{r} clust <- build_clusters(net, k = 2, method = "ward.D2", seed = 42) summary(clust) ``` ### 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. ```{r cluster-diagnostics} diag <- cluster_diagnostics(clust) diag ``` ```{r cluster-diagnostics-plot, fig.alt = "Per-observation silhouette by cluster"} plot(diag, type = "silhouette") ``` `as.data.frame(diag)` returns the per-cluster table for downstream use. ## Mixture Markov Models Instead of clustering sequences based on how similar they are to one another, we can cluster them together based on their transition dynamics. Mixture Markov models (MMM) fit separate Markov models, and sequences are assigned to the cluster whose transition structure best matches their observed behavior. To implement MMM, we can use `build_mmm()`, which returns a `net_mmm` object with full model details: ```{r} mmm_fit <- build_mmm(net, k = 2) summary(mmm_fit) ``` The `net_mmm` object contains posterior probabilities, model fit statistics (BIC, AIC, ICL), and per-cluster transition matrices in `$models`. `cluster_diagnostics()` works on it too, with MMM-specific columns (mixing share, average posterior probability, per-cluster classification error): ```{r mmm-diagnostics} diag_mmm <- cluster_diagnostics(mmm_fit) diag_mmm ``` The `posterior` plot type shows the distribution of each sequence's max posterior probability, faceted by cluster. Tight distributions near 1 indicate clean assignment; spread-out distributions indicate fuzzy clusters. ```{r mmm-diagnostics-plot, fig.alt = "Posterior certainty per MMM cluster"} plot(diag_mmm, type = "posterior") ``` ## Building Networks per Cluster Once sequences are grouped, the goal is usually to compare *how transitions behave within each group*. A single transition network estimated on all sequences averages across the groups and hides the contrast. Estimating one network per cluster preserves each group's distinct dynamics — same nodes, different edge weights. The result is a `netobject_group`: a list of `netobject`s sharing the same node set, one per cluster. Every group network supports the same downstream operations as a single `netobject` — plotting via `plot()`, resampling via `bootstrap_network()`, group-level edge comparison via `permutation()`. Two paths produce the same `netobject_group`: 1. **Manual** — fit clustering with `build_clusters()`, then pass the result to `build_network()` to build one network per cluster. 2. **Shortcut** — `cluster_network()` and `cluster_mmm()` collapse both steps into one call. The distance-based shortcut (`cluster_network()`) groups sessions by sequence similarity. The model-based shortcut (`cluster_mmm()`) groups them by which Markov model best explains each sequence's transitions. | Function | Clustering method | Returns | |----------|------------------|---------| | `cluster_network()` | Distance-based (Hamming, LCS, etc.) | `netobject_group` | | `cluster_mmm()` | Model-based (MMM) | `netobject_group` | ### 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. ```{r cluster-networks} clust <- build_clusters(net, k = 2, method = "ward.D2") cluster_net <- build_network(clust) cluster_net ``` ```{r cluster-networks-plot, fig.alt = "Per-cluster transition networks"} plot(cluster_net) ``` ### Shortcut: `cluster_network()` (distance-based) Combines `build_clusters()` + `build_network()` into one call. ```{r cluster-network-shortcut} ## `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 ``` ### Shortcut: `cluster_mmm()` (model-based) Model-based equivalent: fits a mixture of Markov models and returns per-cluster networks. ```{r cluster-mmm} grp_mmm <- cluster_mmm(net, k = 2) grp_mmm ``` Both return a `netobject_group` — a list of `netobject`s with clustering info accessible via `attr(, "clustering")`: ```{r} # Access cluster assignments attr(grp_dist, "clustering")$assignments[1:10] attr(grp_mmm, "clustering")$assignments[1:10] # Access individual cluster networks grp_dist[[1]]$weights[1:3, 1:3] ``` ### Comparing Clusters We can compare which transition probabilities differ significantly among clusters using permutation testing: ```{r} comparison <- permutation(grp_dist, iter = 100) ``` ## Workflow Summary The end-to-end path from raw sequence data to validated per-cluster networks: ```{r workflow, eval = FALSE} # 1. Build the network from long-format data net <- build_network(human_long, method = "tna", actor = "session_id", action = "cluster", time = "timestamp") # 2. Sweep the parameter space ch <- cluster_choice(net, k = 2:5, dissimilarity = c("hamming", "lcs", "cosine"), method = c("pam", "ward.D2")) plot(ch, type = "facet", abbrev = TRUE) # 3. Pick a configuration and fit clust <- build_clusters(net, k = 2, dissimilarity = "hamming", method = "ward.D2") # 4. Validate diag <- cluster_diagnostics(clust) diag plot(diag, type = "silhouette") # 5. Build per-cluster networks grp <- build_network(clust) # 6. Optional: model-based second opinion mmm <- build_mmm(net, k = 2) plot(cluster_diagnostics(mmm), type = "posterior") ```