Clustering school-to-work trajectories with DBHC

The data

mvad ships with TraMineR. It follows 712 young people in Northern Ireland month by month from September 1993 to June 1999, recording what each was doing in each of 70 months: employment, further education, higher education, joblessness, school, or training. It comes from a cohort study by McVicar and Anyadike-Danes (2002) on predicting successful and unsuccessful transitions from school to work.

This suits DBHC well. The observations are discrete and unordered, the sequences are long enough that the order of states carries most of the information, and the data set carries background variables that were never shown to the algorithm — so the clusters can be checked against something external.

data("mvad", package = "TraMineR")

mvad.alphabet <- c("employment", "FE", "HE", "joblessness", "school", "training")
mvad.scodes <- c("EM", "FE", "HE", "JL", "SC", "TR")
mvad.seq <- seqdef(mvad[, 17:86], id = mvad$id, alphabet = mvad.alphabet,
                   states = mvad.scodes, cnames = 1:70)
dim(mvad.seq)
#> [1] 712  70
alphabet(mvad.seq)
#> [1] "EM" "FE" "HE" "JL" "SC" "TR"

Columns 17 to 86 are the months from September 1993 onward, the convention TraMineR’s own examples use; it drops the two summer months at the start of the survey. Months are numbered 1 to 70 rather than labelled Sep.93 and so on, because seqHMM wants time indices it can coerce to numeric.

Spelling out alphabet is worth the extra line. Without it seqdef() takes the alphabet from whichever states happen to occur and matches states to them in order, so on a subset where some state never appears the short codes would silently attach to the wrong states. Passing both keeps the mapping fixed no matter which sequences you hand it.

Two things to note. mvad also has a weight column from the survey design; hmm.clust() has no weights argument, so every sequence counts once here. And the six state codes are the alphabet of the mixture: each cluster model emits these six symbols from its own hidden states.

Running the algorithm

hmm.clust() searches for the number of clusters and for the number of hidden states in each cluster model on its own, so the whole run is one call. It is randomised — it picks seeds at random and restarts the EM fits — so set a seed if you want to be able to reproduce a partition.

set.seed(2024)
elapsed <- system.time(res <- hmm.clust(mvad.seq))[["elapsed"]]
round(elapsed, 1)
#> [1] 80.1

res$n.clusters
#> [1] 4
res$sizes
#> [1] 4 4 6 8
table(res$cluster)
#> 
#>   1   2   3   4 
#> 222 162 121 207

Under a minute for 712 sequences of length 70, on one core and with no tuning. The last section says a little more about where that time goes.

How many clusters

The algorithm adds clusters one at a time and stops when the BIC of the partition stops improving, keeping the last partition that did improve.

round(res$bic, 1)
#>        1        2        3        4        5 
#> 318839.3  51490.6  39936.1  16480.9  19853.9

The BIC falls until 4 clusters and rises at 5, so 4 is what is returned. Note the first value is enormous: the one-cluster partition is built from a seed of three sequences and never refitted on the rest, so it is not a fitted model in the way the others are and it is not a useful point of comparison.

This stopping rule is greedy, and on these data the BIC path is not monotone beyond the stopping point. Running with K.max = 8 gives

      1        2        3        4        5        6        7        8
215902.6  64972.3  47082.3  31907.0  16631.1  28127.8  16991.4  15940.1

so the search would stop at 5 clusters while 8 actually scores lower. The rule finds a local optimum, not a global one. If the number of clusters matters to your conclusions, set K.max and look at the whole path rather than trusting the stop.

Be careful what you do with that path, though. The BIC keeps rewarding more clusters on these data well past the point where the extra clusters mean anything, and the last section shows a run that reaches a lower BIC still with a cluster containing a single sequence. Read the path as a guide to the range worth examining, not as an answer to be minimised.

What the clusters are

The clearest summary is how long the members of each cluster spend in each state.

cl <- as.integer(res$cluster)
profile <- t(sapply(seq_len(res$n.clusters), function(k)
  round(colMeans(seqistatd(mvad.seq[cl == k, ])), 1)))
rownames(profile) <- paste0("cluster ", seq_len(res$n.clusters))
profile
#>             EM   FE   HE   JL   SC   TR
#> cluster 1 48.4 15.7  0.0  2.1  0.0  3.8
#> cluster 2 39.3  0.0  0.0 10.4  0.0 20.3
#> cluster 3 15.1  0.0 26.3  4.6 24.0  0.0
#> cluster 4 17.7 22.5 13.5  6.4  5.7  4.3

Together with where the trajectories start and end:

modal <- function(x) names(which.max(table(x)))
share <- function(x) round(100 * max(table(x)) / length(x))
month <- function(k, t) mvad.seq[cl == k, t]

endpoints <- data.frame(
  n         = as.vector(table(cl)),
  states    = res$sizes,
  starts    = sapply(seq_len(res$n.clusters), function(k) modal(month(k, 1))),
  pct.start = sapply(seq_len(res$n.clusters), function(k) share(month(k, 1))),
  ends      = sapply(seq_len(res$n.clusters), function(k) modal(month(k, 70))),
  pct.end   = sapply(seq_len(res$n.clusters), function(k) share(month(k, 70))))
endpoints
#>     n states starts pct.start ends pct.end
#> 1 222      4     FE        64   EM      95
#> 2 162      4     TR        83   EM      78
#> 3 121      6     SC        96   HE      49
#> 4 207      8     FE        64   EM      48

The percentages matter as much as the modal state. Two clusters can share a modal destination and still be quite different: one may send nearly everyone there while the other only sends half, with the rest scattered.

Four routes through the six states, all of which are recognisable:

• one that starts in further education and moves into employment, which by the last month is where nearly all of its members are;
• one that starts in training and also ends in employment, but passes through spells of joblessness on the way;
• one that starts in school and moves into higher education, spending almost no time in employment inside the observation window;
• one mixed route that spreads its time across further education, employment and higher education, and whose members are the least likely to have settled into employment by the end.

Cluster numbering depends on the seed, so refer to clusters by their profile rather than by index when you write results up.

seqdplot(mvad.seq, group = cl, border = NA, cex.legend = 0.8,
         xlab = "month", with.legend = "right")

Reading a cluster model

The point of a mixture of HMMs, as opposed to a distance-based clustering, is that each cluster comes with a generative model you can read. emission.heatmap() shows which symbol each hidden state emits and transition.heatmap() shows how the hidden states follow one another, with the initial probabilities on the bottom row.

k <- 1
emission.heatmap(res$partition[[k]]$emission_probs)

transition.heatmap(res$partition[[k]]$transition_probs,
                   res$partition[[k]]$initial_probs)

On these data the emission matrices come out close to deterministic: each hidden state emits one symbol with probability near 1. That makes the models unusually easy to read. A hidden state is then just a spell in an observed state, and the transition matrix is the route between spells, with the diagonal giving how long spells last.

round(res$partition[[k]]$emission_probs, 2)
#>            symbol_names
#> state_names EM FE HE JL SC TR
#>     State 1  0  0  0  0  0  1
#>     State 2  0  0  0  1  0  0
#>     State 3  1  0  0  0  0  0
#>     State 4  0  1  0  0  0  0
round(res$partition[[k]]$transition_probs, 2)
#>          to
#> from      State 1 State 2 State 3 State 4
#>   State 1    0.94    0.01    0.03    0.02
#>   State 2    0.03    0.83    0.12    0.02
#>   State 3    0.00    0.00    1.00    0.00
#>   State 4    0.01    0.01    0.03    0.95

Look for near-1 entries on the diagonal: a hidden state whose self-transition is 1.00 is absorbing, so the trajectories that reach it never leave within the observation window. In the further-education cluster the employment state is absorbing, which is what the 95% figure in the last month reflects. Where two hidden states emit the same symbol and pass back and forth between each other, they are jointly modelling one long spell in that symbol rather than two distinct phases.

Do the clusters mean anything

None of the background variables in mvad were available to the algorithm, so they are a fair test of whether the partition picked up something real.

covariates <- c("gcse5eq", "Grammar", "funemp", "fmpr", "male", "catholic",
                "livboth")
external <- t(sapply(covariates, function(v) {
  pct <- round(100 * tapply(mvad[[v]] == "yes", cl, mean), 1)
  c(pct, overall = round(100 * mean(mvad[[v]] == "yes"), 1),
    p.value = signif(suppressWarnings(chisq.test(table(cl, mvad[[v]]))$p.value), 2))
}))
colnames(external)[seq_len(res$n.clusters)] <- paste0("cl", seq_len(res$n.clusters))
external
#>           cl1  cl2  cl3  cl4 overall p.value
#> gcse5eq  25.2 13.0 64.5 50.7    36.5 1.2e-23
#> Grammar  12.6 10.5 43.0 15.5    18.1 1.9e-13
#> funemp   14.4 24.1 11.6 15.5    16.4 2.1e-02
#> fmpr     21.2 17.3 34.7 28.0    24.6 2.8e-03
#> male     58.1 62.3 42.1 43.0    52.0 7.6e-05
#> catholic 47.3 46.9 52.9 47.8    48.3 7.4e-01
#> livboth  59.0 62.3 68.6 65.7    63.3 2.8e-01

gcse5eq records whether someone left school with five or more GCSEs at grades A–C, the qualification threshold that gates entry to higher education in this cohort. It separates the clusters sharply, and in the direction you would expect: lowest in the training route, highest in the school-to-higher-education route, with the two employment routes in between. Grammar behaves the same way. Father’s unemployment is highest in the training cluster, which is also the cluster with the most time in joblessness.

catholic and livboth are not associated with the partition at all. That is worth saying out loud: the clusters track educational pathway, not the demography of the sample, which is the substantive finding McVicar and Anyadike-Danes were after.

Practical notes

Restart it, but do not follow the BIC off a cliff. The search is randomised and a single run is a sample from it, not the answer. Over 25 seeds on these data it returned anywhere between 2 and 8 clusters:

Restarting reliably improves the BIC. Drawing runs at random from these 25, the median final BIC is 21710 from a single run, 17734 from the best of five, and 16481 from the best of ten. If you only take one thing from this section, take that: one run is not enough.

best <- NULL
for (s in 1:10) {
  set.seed(s)
  candidate <- hmm.clust(mvad.seq)
  score <- candidate$bic[candidate$n.clusters]
  if (is.null(best) || score < best$bic[best$n.clusters]) best <- candidate
}

But the lowest BIC is not automatically the partition you want. The best of these 25 runs scores 16150 with eight clusters — better than the four-cluster solution used above — and two of its eight clusters hold 1 and 10 sequences. The BIC penalty is not strong enough here to stop the search buying likelihood with clusters that fit a handful of sequences almost exactly. Check the cluster sizes before accepting a partition, whatever its score.

Nor do the runs agree with each other. Comparing each run’s partition to the lowest-BIC one with the adjusted Rand index, most score between 0.1 and 0.3, where 0 is chance agreement and 1 is identity; the four-cluster solution scores 0.18. A better BIC does go with closer agreement, but only loosely (correlation -0.42). These partitions are genuinely different carvings of the same data, not noisy versions of one underlying answer.

So treat the number of clusters as a modelling choice that the BIC informs rather than settles. The four-cluster solution is the one shown here because its clusters are substantial, they correspond to recognisable routes, and they are supported by the background variables. Worth knowing: every partition in the study separates qualifications strongly, at p below 1e-20 regardless of the number of clusters. The external variables confirm there is real structure in these trajectories; they do not tell you how finely to cut it.

Watch the smoothing. Cluster models assign a small probability to symbols they never emitted, controlled by smoothing. Without it a sequence containing a symbol its model rates impossible gets a likelihood of zero and cannot be assigned. With long sequences and a small alphabet the default is usually fine, but if you see warnings about infeasible models, raise it.

Sequences may be weighted. hmm.clust() treats every sequence equally. If your data come from a survey with design weights, as mvad does, the partition is unweighted and any population statement drawn from it needs care.

Performance

A run of this size is comfortable on a laptop. Fixing K.max = 4 so the searches are comparable and taking the median of eight seeds per size, since a single run is a poor estimate of the runtime:

The two halves of the table behave quite differently.

The redistribution sweep — evaluating every sequence against every cluster model, the step the algorithm repeats most often — is clean. It is linear in the number of sequences, and the per-sequence cost settles at about 0.09 ms, so placing all 712 takes 64 milliseconds. Sweeps are not where the time goes.

The total is dominated by the EM fits and is much less predictable. Notice how wide the ranges are: the same 400 sequences take anywhere from 5 to 50 seconds depending only on the seed. That is because the algorithm chooses its own model sizes, and fitting an HMM costs roughly the square of its number of hidden states. Regressing the runtime on both terms, each additional sequence costs about 0.04s and each additional hidden state about 2.4s, and both are significant. So the runtime tracks how much structure the search decides is in your data as well as how much data there is — budget for the spread, not for the median.

Two consequences worth planning around. Long sequences are cheaper per observation than many short ones, since sequence length enters linearly while the sequence count carries the per-sequence overhead. And if a run is taking longer than you want, the lever is K.max, which bounds the most expensive part of the search.

The script that produces this table ships with the package:

source(system.file("benchmark", "scaling.R", package = "DBHC"))

References

McVicar, D. and Anyadike-Danes, M. (2002). Predicting Successful and Unsuccessful Transitions from School to Work by Using Sequence Methods. Journal of the Royal Statistical Society, Series A, 165(2), 317–334.

Budel, G., Frasincar, F. and Boekestijn, D. (2024). DBHC: Discrete Bayesian HMM Clustering. International Journal of Machine Learning and Cybernetics.