Longitudinal engagement data have a structural problem: eight indicators per student per course are too many dimensions to compare students on, and averaging them across time destroys exactly the thing that matters: how each student changes. The VaSSTra method solves this in three compressions. Multivariate observations become a small set of interpretable states; each student’s states, ordered in time, become a sequence; and students whose sequences resemble each other become trajectories, a person-centered summary of change that a single regression coefficient cannot express.
This tutorial runs the complete workflow on the longitudinal engagement data from the VaSSTra chapter of Learning Analytics Methods and Tutorials, in this order:
| Task | Verb |
|---|---|
| Complete analysis | vasstra() |
| Compare cluster counts | evaluate(), plot() |
| Fit statistics | fit_indices() |
| Rename groups | set_labels() |
| Tidy results | as.data.frame(), summary() |
| Individual steps | step1_states() … step4_describe() |
| Explicit grids | state_choices(), trajectory_choices() |
The engagement data ship analysis-ready: 142 students at
eight course positions, with eight course-standardized indicators and
attached role metadata, so no column plumbing is needed. The chapter’s
substantive decision (three latent profiles named Disengaged, Average,
and Active) is expressed entirely by the labels: three labels mean three
states, estimated by latent profile analysis (mclust "EEI",
the chapter’s tidyLPA Model 1) because that is the default state method.
The LCS distance with Ward clustering follows the chapter’s trajectory
analysis.
library(VaSSTra)
data("engagement", package = "VaSSTra")
fit <- vasstra(
engagement,
state_labels = c("Disengaged", "Average", "Active"),
state_colors = c(Disengaged = "#e01a4f", Average = "#FFDF3C", Active = "#28a41f"),
dissimilarity = "lcs",
cluster_method = "ward.D2",
positive_states = "Active",
negative_states = "Disengaged"
)
#> Using n_states = 3 to match the supplied labels.
#> Selected n_trajectories = 3 (lcs + ward.D2, silhouette = 0.514); see `diagnostics$selection`.
fit
#> VaSSTra Analysis
#> 142 subjects | 8 times | 3 states | 3 trajectories
#> lcs + ward.D2 | silhouette 0.514state_colors pins the state palette on the fit, so every
plot below is drawn on the chapter’s colours. A companion argument,
state_order, fixes the order states read in every plot;
here the discovery order is already low-to-high, so it is left at its
default. Both are optional.
The one remaining automated decision, the number of trajectories, is
reported in the message and stored with its full comparison table in the
fit, so it can be defended later rather than taken on faith. Leaving the
labels out too (vasstra(engagement)) automates the state
count as well; automation never selects a solution whose smallest group
holds under five percent of observations, the conventional threshold
below which latent classes are treated as spurious.
A fitted clustering is a claim, and the claim needs two kinds of
evidence: that the chosen count beats its neighbors, and that the chosen
solution is internally sound. evaluate() provides both in
one object: a candidate table across counts and a per-group quality
table for the fitted solution.
evaluation <- evaluate(fit)
evaluation
#> VaSSTra clustering evaluation: states (lpa)
#> Fitted: 3 states | silhouette 0.260
#> Candidates:
#> n_states method lpa_model silhouette bic aic classification_entropy
#> 2 lpa EEI 0.35 21679 21553 0.90
#> 3 lpa EEI 0.26 20120 19949 0.89
#> 4 lpa EEI 0.20 19362 19146 0.89
#> 5 lpa EEI 0.20 19009 18747 0.87
#> 6 lpa EEI 0.18 18601 18294 0.89
#> min_size max_size size_ratio eligible best fitted
#> 509 627 1.2 TRUE FALSE FALSE
#> 266 551 2.1 TRUE FALSE TRUE
#> 152 414 2.7 TRUE FALSE FALSE
#> 155 367 2.4 TRUE TRUE FALSE
#> 33 350 10.6 FALSE FALSE FALSE
#> Fitted states:
#> state n proportion silhouette
#> Disengaged 266 0.23 0.33
#> Average 551 0.48 0.23
#> Active 319 0.28 0.26
#>
#> VaSSTra clustering evaluation: trajectories (lcs + ward.D2)
#> Fitted: 3 trajectories | silhouette 0.514
#> Candidates:
#> n_trajectories dissimilarity method silhouette mean_within_distance min_size
#> 2 lcs ward.D2 0.44 6.4 41
#> 3 lcs ward.D2 0.51 4.5 28
#> 4 lcs ward.D2 0.40 4.2 28
#> 5 lcs ward.D2 0.37 3.9 8
#> 6 lcs ward.D2 0.35 3.6 8
#> max_size size_ratio eligible best fitted
#> 101 2.5 TRUE FALSE FALSE
#> 73 2.6 TRUE TRUE TRUE
#> 43 1.5 TRUE FALSE FALSE
#> 43 5.4 TRUE FALSE FALSE
#> 43 5.4 TRUE FALSE FALSE
#> Fitted trajectories:
#> trajectory n proportion mean_within_distance silhouette
#> Trajectory 1 41 0.29 4.3 0.58
#> Trajectory 2 73 0.51 5.0 0.44
#> Trajectory 3 28 0.20 3.9 0.61The two rows of panels answer different questions. The trajectory selection curve peaks at the fitted three groups (silhouette 0.51 against 0.44 and 0.40 for two and four), so the automated choice is also the local optimum, the strongest position a selection can be in. The state curve instead declines from two: with LPA the count was fixed by the three labels, and the panel makes the cost of that substantive choice visible (0.26 at three against 0.35 at two) instead of hiding it. The silhouette-by-group bars then locate where cohesion comes from: the two flanking groups are the most compact (0.58 and 0.61), while the large middle trajectory (0.44) is the first place to look when judging whether three groups over-merge distinct journeys.
For the model-based evidence reviewers expect,
fit_indices() returns the complete information-criterion
family as one tidy row: the same content as tidyLPA’s fit table, plus
the classification-quality columns that usually have to be assembled by
hand.
fit_indices(fit)
#> VaSSTra fit indices: states (lpa)
#> n_states method lpa_model log_likelihood n_parameters aic bic sabic caic
#> 3 lpa EEI -9940 34 19949 20120 20012 20154
#> awe clc kic icl entropy prob_min prob_max silhouette min_size max_size
#> 20461 20145 19986 20384 0.89 0.94 0.96 0.26 266 551
#> size_ratio
#> 2.1Three of these numbers carry the interpretation. Entropy (0.89) says
the posterior classification is sharp overall. prob_min
(0.94) says even the worst class assigns its members with high
confidence. Entropy summarizes the whole matrix and
prob_min guards its weakest class, so reporting both
prevents a good average from hiding a bad class. The criterion family
(AIC through ICL) is there for the comparison table:
fit_indices(fit, compare = TRUE)
#> VaSSTra fit indices: states (lpa)
#> n_states method lpa_model log_likelihood n_parameters aic bic sabic caic
#> 2 lpa EEI -10752 25 21553 21679 21600 21704
#> 3 lpa EEI -9940 34 19949 20120 20012 20154
#> 4 lpa EEI -9530 43 19146 19362 19226 19405
#> 5 lpa EEI -9322 52 18747 19009 18844 19061
#> 6 lpa EEI -9086 61 18294 18601 18407 18662
#> awe clc kic icl entropy prob_min prob_max silhouette min_size max_size
#> 21930 21655 21581 21831 0.90 0.97 0.97 0.35 509 627
#> 20461 20145 19986 20384 0.89 0.94 0.96 0.26 266 551
#> 19794 19398 19192 19700 0.89 0.93 0.97 0.20 152 414
#> 19531 19114 18802 19480 0.87 0.87 0.96 0.20 155 367
#> 19213 18633 18358 19062 0.89 0.88 0.97 0.18 33 350
#> size_ratio eligible best fitted
#> 1.2 TRUE FALSE FALSE
#> 2.1 TRUE FALSE TRUE
#> 2.7 TRUE FALSE FALSE
#> 2.4 TRUE TRUE FALSE
#> 10.6 FALSE FALSE FALSEEvery criterion keeps falling through six classes, a nearly linear
decline with no elbow. This is precisely why the chapter’s three
profiles rest on interpretability and the classification-quality columns
rather than on a BIC minimum. The table makes that argument honest: the
best flag marks what pure criteria would pick,
fitted marks what was chosen, and the two need not
agree.
The states are only useful if they can be named, and naming requires knowing how the profiles differ. The overview draws the four state views together; each answers one question.
The profile lines answer what makes each state itself: Active sits around one standard deviation above course average on every indicator, Disengaged about one below, Average on the line. Engagement here is a level, not a trade-off between, say, forum and lecture activity, which is exactly what justifies ordered labels. The bars carry the same means but grouped by indicator, which is the view for spotting whether any single indicator breaks the ordering (none does). The heatmap adds the exact magnitudes the lines only sketch: Disengaged is most extreme on session count (−1.27) and least on lecture views (−0.81). The sizes panel grounds prevalence: Average is the plurality (48%), and no state is rare enough to question.
The default sequence view is the heatmap: every aligned sequence as one row, rows ordered by similarity, so persistence shows as long unbroken bands and instability as speckle. It is the one view that keeps all 142 individual journeys visible at once.
Two further views split the temporal question the heatmap compresses. The distribution answers the aggregate question (does the composition of the cohort drift over the eight positions), while the index plot answers the individual one on a fixed subject order. The distinction matters because a flat distribution can hide heavy individual churn.
Grouping the sequences (LCS distance, Ward linkage) yields three trajectories of 41, 73, and 28 students. The grouped index plot is the core evidence that the groups mean something: one group is dominated by Active rows throughout, one by Disengaged, and the largest mixes Average with excursions. These are visibly different journeys, not arbitrary partitions of similar ones.
The three groups line up with the three states, so we can name them
now by the state each mostly holds. set_labels() renames in
place: every later table, sequence, and plot follows, and no fitted
value changes.
Passing several views in type composes them into one
figure with one row per trajectory and one column per
view, so each group’s transition network sits beside the
sequences it summarises and their changing composition. This is the
standard per-cluster VaSSTra display in one call. A shared legend labels
every panel at once, and the transition column requires the suggested
cograph package.
The rows now read as whole stories. The Active group’s network is dominated by the Active self-loop and its distribution stays green, and the Disengaged group mirrors it in red. The mixed middle group shows the busiest network: the many off-diagonal transitions that its higher entropy and mobility counts quantify in the next section.
The per-trajectory summary turns those pictures into numbers, one row per trajectory.
as.data.frame(fit, unit = "trajectory")
#> trajectory n proportion mean_within_distance n_observed unique_states
#> 1 Mostly active 41 0.29 4.3 8 1.9
#> 2 Mostly average 73 0.51 5.0 8 2.0
#> 3 Mostly disengaged 28 0.20 3.9 8 1.8
#> transitions entropy complexity volatility integrative_potential
#> 1 1.8 0.42 0.33 0.44 0.77
#> 2 2.3 0.47 0.39 0.50 0.11
#> 3 1.9 0.40 0.32 0.44 0.00
#> negative_exposure
#> 1 0.0095
#> 2 0.1317
#> 3 0.8204The two time-weighted columns exist because
positive_states and negative_states were
declared in the fit: integrative potential is the late-weighted share of
time in Active, negative exposure the same for Disengaged. Here they
separate the groups almost completely (0.77 versus 0.11 versus 0.00 for
integrative potential). Later time points count more, so these indices
read as where the student is heading, which raw state
proportions cannot say. Entropy and transition counts add the stability
dimension: the large middle group is the most mixed in state and also
the most mobile.
Per-student indices for downstream modeling (one row per student with trajectory membership, entropy, complexity, and both exposure indices) come from the same fit:
head(as.data.frame(fit))
#> user_id trajectory n_observed unique_states transitions entropy
#> 1 D2C5F64E Mostly active 8 3 3 0.99
#> 2 D7E3D0DC Mostly active 8 3 4 0.82
#> 3 2926CF64 Mostly disengaged 8 1 0 0.00
#> 4 985AFF35 Mostly disengaged 8 2 2 0.60
#> 5 F89100C7 Mostly active 8 3 4 0.89
#> 6 AB944042 Mostly disengaged 8 2 2 0.34
#> complexity volatility integrative_potential negative_exposure
#> 1 0.65 0.71 0.42 0.083
#> 2 0.68 0.79 0.58 0.083
#> 3 0.00 0.17 0.00 1.000
#> 4 0.41 0.48 0.00 0.556
#> 5 0.71 0.79 0.39 0.194
#> 6 0.31 0.48 0.00 0.889Interpretation happens after fitting, so labels stay
assignable at any point without refitting, which is exactly what Section
5 did once the groups were visible. set_labels() renames in
place, propagating the names through every table, sequence, and plot,
including the recorded positive and negative states, and changing no
fitted value. A full vector renames every group; a named vector renames
only what it names, which is handy for refining a single label after a
closer look:
# rename only one group, by name
fit <- set_labels(fit, trajectories = c("Mostly average" = "Fluctuating"))The named trajectories now flow through every summary of the fit:
summary(fit)
#> trajectory n proportion mean_within_distance n_observed unique_states
#> 1 Mostly active 41 0.29 4.3 8 1.9
#> 2 Mostly average 73 0.51 5.0 8 2.0
#> 3 Mostly disengaged 28 0.20 3.9 8 1.8
#> transitions entropy complexity volatility integrative_potential
#> 1 1.8 0.42 0.33 0.44 0.77
#> 2 2.3 0.47 0.39 0.50 0.11
#> 3 1.9 0.40 0.32 0.44 0.00
#> negative_exposure
#> 1 0.0095
#> 2 0.1317
#> 3 0.8204Everything above is also four explicit steps, each with the same
automated defaults. This is the route to take when intermediate objects
are needed, when states come from another method (enter at
step2_sequences() with a state column), or when the method
itself is under comparison:
states <- step1_states(engagement, n_states = 3)
sequences <- step2_sequences(states)
trajectories <- step3_trajectories(sequences, dissimilarity = "lcs",
method = "ward.D2")
description <- step4_describe(trajectories, positive_states = "Active",
negative_states = "Disengaged")When the clustering method is itself in question,
state_choices() and trajectory_choices() fit
the full grid and return one tidy row per candidate. Nothing is selected
silently. Recommendations are marked within each method (silhouette for
hard clustering, BIC for LPA, since model-based criteria only compare
within LPA), and the inspected candidate is fitted with
fit_state_choice():
state_options <- state_choices(engagement, n_states = 2:4,
method = c("lpa", "kmeans"))
state_options
#> VaSSTra state choices
#> 6 candidates | 6 successful | 2 recommended
#> candidate_id n_states method lpa_model recommendation_criterion silhouette
#> 3 4 lpa EEI bic 0.20
#> 4 2 kmeans <NA> silhouette 0.35
#> bic min_size eligible
#> 19362 152 TRUE
#> NA 519 TRUE
plot(state_options, metric = "silhouette")states <- fit_state_choice(state_options, n_states = 3, method = "lpa",
labels = c("Disengaged", "Average", "Active"))The choice plot and evaluate() show the same silhouette
on the same scale, so the method comparison here and the count
comparison there compose into one argument. Method overview and formal
background: the VaSSTra
chapter.