Short R Tutorial: Validating Sequence Analysis Typologies Using Parametric Bootstrap'

In WeightedCluster: Clustering of Weighted Data

Introduction

knitr::knit_hooks$set(time_it = local({
  now <- NULL
  function(before, options) {
    if (before) {
      # record the current time before each chunk
      now <<- Sys.time()
    } else {
      # calculate the time difference after a chunk
      res <- difftime(Sys.time(), now, units = "secs")
      # return a character string to show the time
      paste("Time for this code chunk to run:", round(res,
        2), "seconds")
    }
  }
}))
knitr::opts_chunk$set(dev = "png", dev.args = list(type = "cairo-png"), time_it=TRUE)

This document provides a very quick introduction to the R code needed to use parametric bootstraps for typology validation in sequence analysis. Readers interested in the methods and the exact interpretation of the results are referred to:

• Studer, M. (2021). Validating Sequence Analysis Typologies Using Parametric Bootstraps. Sociological Methodology 51(2), 290--318. https://doi.org/10.1177/00811750211014232

You are kindly asked to cite the above reference if you use the methods presented in this document.

Warning!! To avoid lengthy computations (and overloading the CRAN server), we restricted the number of bootstraps to 50. We recommend using a higher value (i.e., 1000).

Let's start by setting the seed for reproducible results.

set.seed(1)

Creating the State Sequence Object

For this example, we use the mvad dataset. Let's start with the creation of the state sequence object.

## Loading the TraMineR library
library(TraMineR)
## Loading the data
data(mvad)

## State properties
mvad.alphabet <- c("employment", "FE", "HE", "joblessness", "school", "training")
mvad.lab <- c("employment", "further education", "higher education", "joblessness", "school", "training")
mvad.shortlab <- c("EM","FE","HE","JL","SC","TR")

## Creating the state sequence object
mvad.seq <- seqdef(mvad, 17:86, alphabet = mvad.alphabet, states = mvad.shortlab, labels = mvad.lab, xtstep = 6)

Creating the typology

We will now create a typology using cluster analysis. Readers interested in more detail are referred to the WeightedCluster library manual (also available as a vignette), which goes into the details of the creation and computation of cluster quality measures.

We start by computing dissimilarities with the seqdist function using the Hamming distance. We then use Ward clustering to create a typology of the trajectories. For this step, we recommend the use of the fastcluster library, which considerably speed up the computations.

## Using fastcluster for hierarchical clustering
library(fastcluster)
## Distance computation
diss <- seqdist(mvad.seq, method="HAM")
## Hierarchical clustering
hc <- hclust(as.dist(diss), method="ward.D")

We can now compute several cluster quality indices using as.clustrange function from two to ten groups.

# Loading the WeightedCluster library
library(WeightedCluster)
# Computing cluster quality measures.
clustqual <- as.clustrange(hc, diss=diss, ncluster=10)
clustqual

Parametric Bootstrap

Parametric bootstrap aims to provide baseline values obtained by clustering similar but non-clustered data [@Studer2021]. This can be computed using the seqnullcqi function with the following parameters:

• R: number of bootstraps.
• model: The null model (see table below
• seqdist.args: list of arguments passed to seqdist (should be identical to first call to seqdist).
• hclust.method: hierarchical clustering method (should be identical to orginal clustering).
• kmedoid: If TRUE, use PAM (and the wcKMedRange function) instead of hierarchical clustering.
• parallel: If TRUE, use parallel computing to speed up the computations.

Combined Randomization

The following R code estimate expected values of the cluster quality indices when clustering similar sequences that are not clustered according to the "combined" model, using the Hamming distance and Ward hierarchical clustering. We set parallel=TRUE to use parallel computing. You can use progressbar=TRUE to show a progress bar and an estimation of the computation remaining time (not meaningful here within a document):

bcq.combined <- seqnullcqi(mvad.seq, clustqual, R=50, model="combined", seqdist.args=list(method="HAM"), hclust.method="ward.D", parallel = TRUE)

Once the parametric bootstrap is computed (may take a while...), the results are stored in the bcq.combined object. Printing the object (just by writing its name), already provides several information, the standardized cluster quality indices and the associated inconclusive intervals. Here, 2, 9 and 10 groups stand out.

bcq.combined

To get non-standardized values, use norm=FALSE. Notice that the ASW inconclusive intervals are well below the values recommended by Kaufman and Rousseeuw (over 0.5).

print(bcq.combined, norm=FALSE)

Several plots can then be used to inspect the results using the plot command and the type argument. First, one can look at the sequences generated by the null model by using type="seqdplot".

plot(bcq.combined, type="seqdplot")

The overall distribution of the CQI values can be plotted using type="density". In this case, one also needs to specify the CQI to be used. All tested number of groups are found to be significant. Any CQI computed by as.clustrange() can be used here. To show the density of the average silhouette width ("ASW"), one can use:

plot(bcq.combined, stat="ASW", type="density")

By using type="line", we plot the obtained and bootstrapped CQI values depending on the number of groups. Here again

plot(bcq.combined, stat="ASW", type="line")

Randomized Sequencing

To use another null model, one needs to change the model argument of the seqnullcqi function. The randomized sequencing keep the duration attached to each state, but randomizes the ordering of the spells. It can be used to uncover sequencing structure of the data.

bcq.seq <- seqnullcqi(mvad.seq, clustqual, R=50, model="sequencing", seqdist.args=list(method="HAM"), hclust.method="ward.D", parallel = TRUE)

We can then plot the results as before. Notice that solutions between 3 and 6 are below the critical line.

plot(bcq.seq, stat="ASW", type="line")

Randomized Duration

The randomized duration keeps the same ordering of the states, but randomizes the time spent in each spell. It can be used to uncover the duration-related structure of the data.

bcq.dur <- seqnullcqi(mvad.seq, clustqual, R=50, model="duration", seqdist.args=list(method="HAM"), hclust.method="ward.D", parallel = TRUE)

We can then plot the results as before. The solutions 3 and 4 groups solutions are below the "significance line". Otherwise, the ranking of the solutions is the same.

plot(bcq.dur, stat="ASW", type="line")

State Independence

The state independence null model generates sequence, position by position, independently of the previous state. This is a quite unrealistic assumption for longitudinal data, but a common one in statistical modeling.

bcq.stateindep <- seqnullcqi(mvad.seq, clustqual, R=50, model="stateindep", seqdist.args=list(method="HAM"), hclust.method="ward.D", parallel = TRUE)

Bootstrapped CQI values are extremely low compared to our clustering, meaning that we have a strong longitudinal structure (not surprising!).

plot(bcq.stateindep, stat="ASW", type="line")

First-order Markov Null Model

The first-order Markov null model generates sequences using time-invariant transition rates. As a result, the generated sequences are often quite similar to the observed ones. This model can uncover structure stemming from time-dependent transition rates.

bcq.Markov <- seqnullcqi(mvad.seq, clustqual, R=50, model="Markov", seqdist.args=list(method="HAM"), hclust.method="ward.D", parallel = TRUE)

plot(bcq.Markov, stat="ASW", type="line")

Choosing a Solution

The various null models lead to the same conclusions and ranking of the solutions. Solutions between 3 and 6 groups were not always above the critical lines (in the sequencing null model for instance), and can be avoided. We generally saw good clustering quality for a clustering in 9 groups. The solution is shown below.

seqdplot(mvad.seq, clustqual$clustering$cluster9, border=NA)

WeightedCluster documentation built on April 24, 2025, 3:01 a.m.