A Modern Approach to Transition Analysis and Process Mining with Markov Models in Education

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ A Modern Approach to Transition Analysis and Process Mining with Markov Models in Education © The Author(s) 2024 Published version Helske, Jouni; Helske, Satu; Saqr, Mohammed; López-Pernas, Sonsoles; Murphy, Keefe Helske, J., Helske, S., Saqr, M., López-Pernas, S., & Murphy, K. (2024). A Modern Approach to Transition Analysis and Process Mining with Markov Models in Education. In M. Saqr, & S. López-Pernas (Eds.), Learning Analytics Methods and Tutorials : A Practical Guide Using R (pp. 381-427). Springer. https://doi.org/10.1007/978-3-031-54464-4_12 2024 A Modern Approach to Transition Analysis and Process Mining with Markov Models in Education Jouni Helske, Satu Helske, Mohammed Saqr, Sonsoles López-Pernas, and Keefe Murphy 1 Introduction In the previous two chapters, we have learned about sequence analysis [1, 2] and its relevance to educational research. This chapter presents a closely-related method: Markovian models. Specifically, we focus on a particular type of Markovian model, where the data are assumed to be categorical and observed at discrete time intervals, as per the previous chapters about sequence analysis, although in general Markovian models are not restricted to categorical data. One of the main differences between sequence analysis and Markovian modelling is that the former relies on deterministic data mining, whereas the latter uses probabilistic models [3]. Moreover, sequence analysis takes a more holistic approach by analysing sequences as a whole, whereas Markovian modelling focuses on the transitions between states, their probability, and the reasons (covariates) which explain why these transitions happen. J. Helske () Department of Mathematics and Statistics, University of Jyväskylä, Jyväskylän yliopisto, Finland INVEST Research Flagship Centre, University of Turku, Turku, Finland e-mail: jouni.helske@iki.fi S. Helske INVEST Research Flagship Centre, University of Turku, Turku, Finland Department of Social Research, University of Turku, Turku, Finland M. Saqr · S. López-Pernas School of Computing, University of Eastern Finland, Joensuu, Finland K. Murphy Department of Mathematics and Statistics, Hamilton Institute, Maynooth University, Maynooth, Ireland © The Author(s) 2024 M. Saqr, S. López-Pernas (eds.), Learning Analytics Methods and Tutorials, https://doi.org/10.1007/978-3-031-54464-4_12 381 382 J. Helske et al. We provide an introduction and practical guide to the topic of Markovian models for the analysis of sequence data. While we try to avoid advanced mathematical notations, to allow the reader to continue to other, more advanced sources when necessary, we do introduce the basic mathematical concepts of Markovian models. When doing so, we use the same notation as in the R package ‘seqHMM’ [4], which we also use in the examples. In particular, we illustrate first-order Markov models, hidden Markov models, mixture Markov models, and mixture hidden Markov models with applications to synthetic data on students’ collaboration roles throughout a complete study program. The chapter proceeds to describe the theoretical underpinnings on each method in turn, then showcases each method with code, before presenting some conclusions and further readings. In addition to the aforementioned applications to collaboration roles and achievement sequences, we also provide a demonstration of the utility of Markovian models in another context, namely process mining. In the process mining application, we leverage Markov models and mixture Markov models to explore learning management system logs. Finally, we conclude with a brief discussion of Markovian models in general and provide some recommendations for further reading of advanced topics in this area as a whole. 2 Methodological Background 2.1 Markov Model The simple first-order Markov chain or Markov model (MM) can be used to model transitions between successive states. In the first-order MM, given the current observation, the next observation in the sequence is independent of the past—this is called the Markov property (the order of MM determines on how many previous observations the next observation depends on). For example, when predicting a student’s school success in the fourth year under a first-order model, we only need to consider their success in the third year, while their success in the first and second year give no additional information for the prediction (see Fig. 1 for an illustration). As such, the model is said to be memoryless. As an example, consider the data described in Table 1 which includes four sequences of length ten. The alphabet—that is, the list of all possible states appearing in the data—consists of two types of observed state; low achievement Fig. 1 Illustration of the Markov Model. The nodes . Y1to . Y4refer to states at time points 1 to 4. The arrows indicate dependencies between states A Modern Approach to Transition Analysis and Process Mining with Markov... 383 Table 1 Four example sequences of school achievement with individuals A–D across the rows and years 1–10 across the columns 1 2 3 4 5 6 7 8 9 10 A L L L H L H L H H H B L H H L H L H L L H C H H L H L L H L H H D H H L L L H L L L H Table 2 Transition matrix showing the probabilities of transitioning from one state to another (low or high achievement). The rows and columns describe the origin state and the destination state, respectively . →Low . →High Low .→8/20 = 0.4 12/20 = 0.6 High .→10/16 = 0.625 6/16 = 0.375 success (L) and high achievement success (H). Here, the individuals are assumed to be independent from one another: Say t describes the position in the sequence, or in this example, the year (in other words, here t runs from 1 to 10). If we assume that the probability of observing L or H at any given point t depends on the previous observation only, we can estimate the transition probabilities .aLL (from state L to state L), .aLH (L to H), .aHL (H to L), and .aHH (H to H) by calculating the number of observed transitions from each state to all states and scaling these with the total number of transitions from that state. Mathematically, we can write the transition probability . ars from state r to state s as . ars =P(z t=s|zt−1=r), s,r ∈{L, H}, which simply states that the observed state . ztin year t being L or H depends on which of the two states were observed in the previous year .t−1. For example, to compute .aLH =P(z t=H|zt−1=L), the probability of transitioning from the origin state L to the destination state H, we divide the twelve observed transitions to state H from state L by 20, which is the total number of transitions from L to any state. The basic MM assumes that the transition probabilities remain constant in time (this property is called time-homogeneity). This means, for example, that the probabilities of transitioning from the low-achievement state to the highachievement state is the same in the ninth year as it was in the second year. We can collect the transition probabilities in a transition matrix (which we call A) which shows all of the possible transition probabilities between each pair of origin and destination states, as illustrated in Table 2. For example, when a student has low achievement in year t, they have a 40% probability to have low achievement in year .t+1 and a higher 60% probability to transition to high achievement instead, regardless of the year t. Notice that the probabilities in each row must add up to 1 (or 100%). 384 J. Helske et al. Lastly, we need to define probabilities for the starting states of the sequences, i.e., the initial probabilities . πs=P(z 1=s), s ∈{L, H}. In the example, half of the students have low achievement and the other half have high achievement in the first year, so .πL=πH=0.5. This basic MM is very simple and is often not realistic in the context of educational sciences. We can, however, extend the basic MM in several ways. First of all, we can include covariates to explain the transition and/or initial probabilities. For example, if we think that transitioning from low to high achievement becomes more challenging as the students get older we may add time as an explanatory variable to the model, allowing the probability of transitioning from low to high achievement to decrease in time. We could also increase the order of the Markov chain, accounting for longer histories. This may be more realistic, but at the same time increasing the order makes the model considerably more complex, the more so the longer the history considered. Secondly, one of the most useful extensions is the inclusion of hidden (or latent) states that cannot be observed directly but can be estimated from the sequence of observed states. An MM with time-constant hidden states is typically called the mixture Markov model (MMM). It can be used to find latent subpopulations, or in other words, to cluster sequence data. A model with time-varying hidden states is called the hidden Markov model (HMM), which allows the individuals to transition between the hidden states. Allowing for both time-constant and timevarying hidden states leads to a mixture hidden Markov model (MHMM). Unless otherwise specified, from now on when talking about hidden states we refer always to time-varying hidden states, while time-constant hidden states are referred to as clusters. 2.2 Mixture Markov Model Consider a common case in sequence analysis where individual sequences are assumed to be clustered into subpopulations such as those with typically high and low achievement. In the introductory sequence analysis chapter, the clustering of sequences was performed based on a matrix of pairwise dissimilarities between sequences. Alternatively, we can use the MMM to group the sequences based on their initial and transition probabilities, for example, into those who tend to stay in and transition to high achievement states and those that tend to stay in and transition to low achievement states, as illustrated in Table 3. In MMMs, we have a separate transition matrix .Akfor each cluster k (for .k=1,...,K clusters/subpopulations), and the initial state distribution defines the probabilities to start (and stay) in the hidden states corresponding to a particular cluster. This probabilistic clustering provides group membership probabilities for A Modern Approach to Transition Analysis and Process Mining with Markov... 385 Table 3 Two transition matrices showing the probabilities of transitioning from one state of achievement to another in two clusters of Low achievement and High achievement. The rows and columns describe the origin state and the destination state, respectively (a) Low achievement (b) High achievement Cluster: Low achievement . →Low . →High Cluster: High achievement . →Low . →High Low .→0.8 0.2 Low .→0.6 0.4 High .→0.4 0.6 High .→0.1 0.9 each sequence; these define how likely it is that each individual is a member of each cluster. We can easily add (time-constant) covariates to the model to explain the probabilities of belonging to each cluster. By incorporating covariates in this way we could, for example, find that being in a high-achievement cluster is predicted by gender or family background. However, we note that this is distinct from the aforementioned potential inclusion of covariates to explain the transition and/or initial probabilities. An advantage of this kind of probabilistic modelling approach is that we can use traditional model selection methods such as likelihood-based information criteria or cross-validation for choosing the best model. For example, if the number of subpopulations is not known in advance—as is typically the case—we can compare models with different clustering solutions (e.g., those obtained with different numbers of clusters, different subsets of covariates, or different sets of initial probabilities, for example) and choose the best-fitting model with, for example, the Bayesian information criterion (BIC) [5]. 2.3 Hidden Markov Model The HMM can be useful in a number of cases when the state of interest cannot be directly measured or when there is measurement error in the observations. In HMMs, the Markov chain operates at the level of hidden states, which subsequently generate or emit observed states with different probabilities. For example, think about a progression of a student’s ability as a hidden state and school success as the observed state. We cannot measure true ability directly, but we can estimate the student’s progress by their test scores that are emissions of their ability. There is, however, some uncertainty in how well the test scores represent students’ true ability. For example, observing low test scores at some point in time does not necessarily mean the student has low ability; they might have scored lower than expected in the test due to other reasons such as being sick at that particular time. Such uncertainty can be reflected in the emission probabilities; for example, in the high-ability state students get high test scores eight times out of ten and low test scores with a 20% probability, while in the low-ability state the students get low test scores nine times out of ten and high test scores with a 10% probability. These probabilities are collected in an emission matrix as illustrated in Table 4. 386 J. Helske et al. Table 4 Emission matrix showing the probabilities of each hidden state (low or high ability) emitting each observed state (low or high test scores) Low scores High scores Low ability 0.9 0.1 High ability 0.2 0.8 Fig. 2 Illustration of the HMM. The nodes . Z1to . Z4refer to hidden states at time points 1 to 4, while the nodes . Y1to . Y4refer to observed states. The arrows indicate dependencies between hidden and/or observed states Again, the full HMM is defined by a set of parameters: the initial state probabilities . π s , the hidden state transition probabilities . a rs , and the emission probabilities of observed states .b s (m). What is different to the MM is that in the HMM, the initial state probabilities . π s define the probabilities of starting from each hidden state. Similarly, the transition probabilities .a rs define the probabilities of transitioning from one hidden state to another hidden state. The emission probabilities . b s (m) (collected in an emission matrix B) define the probability of observing a particular state m (e.g., low or high test scores) given the current hidden state s (e.g., low or high ability). When being in a certain hidden state, observed states occur randomly, following the emission probabilities. Mathematically speaking, instead of assuming the Markov property directly on our observations, we assume that the observations are conditionally independent given the underlying hidden state. We can visualise the HMM as a directed acyclic graph (DAG) illustrated in Fig. 2.Here Z are the unobserved states (such as ability) which affect the distribution of the observed states Y (test scores). At each time point t, the state . z t can obtain one of S possible values (there are two hidden states in the example of low and high ability, so .S=2), which in turn defines how . Y t is distributed. 2.4 Mixture Hidden Markov Models Combining the ideas of both time-constant clusters and time-varying hidden states leads to the concept of mixture hidden Markov model (MHMM). Here we assume A Modern Approach to Transition Analysis and Process Mining with Markov... 387 Table 5 Two transition matrices showing the probabilities of transitioning from one state of ability to another in two clusters, the Stayers and the Movers. The rows and columns describe the origin state and the destination state, respectively (a) Stayers (b) Movers Cluster: Stayers . →Low . →High Cluster: Movers . →Low . →High Low .→1 0 Low .→0.6 0.4 High .→0 1 High .→0.3 0.7 Table 6 Two emission matrices showing the probabilities of each hidden state (low or high ability) emitting each observed state (low or high test scores) (a) Stayers (b) Movers Cluster: Stayers Low scores High scores Cluster: Movers Low scores High scores Low ability 0.9 0.1 Low ability 0.7 0.3 High ability 0.1 0.9 High ability 0.2 0.8 that the population of interest consists of a finite number of subpopulations, each with their own HMM with varying transition and emission probabilities. For example, we could expect to find underlying groups which behave differently when estimating the progression of ability through the sequence of test scores, such as those that consistently stay on a low-ability or high-ability track (stayers) and those that move between low and high ability (movers). In this case, we need two transition matrices: the stayers’ transition matrix allows for no transitions while the movers’ transition matrix allows for transitioning between low and high ability, as illustrated in Table 5. Similarly, we need two emission matrices that describe how the observed states are related to hidden states, as illustrated in Table 6. In this example, there is a closer match between low/high ability and low/high test scores in the Stayers cluster in comparison to the Movers cluster. Mathematically, when estimating a MHMM we first fix the number of clusters K, and create a joint HMM consisting of K submodels (HMMs). The number of hidden states does not have to be fixed but can vary by submodel, so that the HMMs have more hidden states for some clusters and fewer for others (in our example, because the transition matrix of the Stayers cluster is diagonal, we could also split the cluster into two single state clusters, one corresponding to low and another to high ability). This can increase the burden of model selection, so often a common number of hidden states is assumed for each cluster for simplicity. In any case, the initial state probabilities of this joint model define how sequences are assigned to different clusters. We estimate this joint model using the whole data and calculate cluster membership probabilities for each individual. The idea of using mixtures of HMMs has appeared in literature under various names with slight variations, e.g., [6], [7], and [4]. Notably, MHMMs inherit from MMMs the ability to incorporate covariates to predict cluster memberships. 388 J. Helske et al. 2.5 Multi-Channel Sequences There are two options to analyse multi-channel (or multi-domain or multidimensional) sequence data with Markovian models. The first option is to combine observed states in different channels into one set of single-channel sequences with an expanded alphabet. This option is simple, and works for MMs, HMMs, MMMs, and MHMMs, but can easily lead to complex models as the number of states and channels increases considerably. The second option, which can only be used when working with HMMs and MHMMs, is to treat the observed states in each channel independently given the current hidden state. This can be easily performed by defining multiple emission probability matrices, one for each channel. The assumption of conditional independence simplifies the model, but is sometimes unrealistic, in which case it is better to resort to the first option and convert the data into single-channel sequences. Both options are discussed further in Chapter 13 [8], a dedicated chapter on multi-channel sequences, where applications of distance-based and Markovian clustering approaches are presented. In this chapter, we henceforth focus on single-channel sequences. 2.6 Estimating Model Parameters The model parameters, i.e. the elements of the initial probability vectors . π, transition probability matrices A, and emission probability matrices B, can be estimated from data using various methods. Typical choices are the Baum-Welch algorithm (an instance of the expectation-maximisation, i.e., the EM algorithm) and direct (numerical) maximum likelihood estimation. It is possible to restrict models, for example, by setting some parameters to fixed values (typically zeros), for example, to make certain starting states, transitions, or emissions impossible. After the parameter estimation, in addition to studying the estimated model parameters upon convergence, we can, for example, compute cluster-membership probabilities for each individual and find the most probable paths of hidden state sequences using the Viterbi algorithm [9]. These can be further analysed and visualised for interpretation. 3 Review of the Literature Markovian methods have been used across several domains in education and have gained renewed interest with the surge in learning analytics and educational data mining. Furthermore, the introduction of specialised R packages (e.g., seqHMM [10]) A Modern Approach to Transition Analysis and Process Mining with Markov... 395 4.2.2 Hidden Markov Models The structure of an HMM is set with the build_hmm() function. In contrast to build_mm(), other build_*() functions such as build_hmm() do not directly estimate the model parameters. For build_hmm(), in addition to observations (an stslist), we need to provide the n_states argument which tells the model how many hidden states to construct. Using again the collaboration roles sequences, if we want to estimate an HMM with two hidden states, we can write: set.seed(1) hidden_markov_model <- build_hmm( observations = roles_seq, n_states = 2 ) The set.seed() call ensures that we will always end up with the same exact initial model with hidden states in the same exact order even though we use random values for the initial parameters of the model (which is practical for reproducibility). We are now ready to estimate the model with the fit_model() function. The HMM we want to estimate is simple, so we rely on the default values and again use the print() method to provide information about the estimated model: fit_hmm <- fit_model(hidden_markov_model) fit_hmm$model Initial probabilities : State 1 State 2 0.657 0.343 Transition probabilities : to from State 1 State 2 State 1 0.9089 0.0911 State 2 0.0391 0.9609 Emission probabilities : symbol_names state_names Isolate Mediator Leader State 1 0.4418 0.525 0.0336 State 2 0.0242 0.478 0.4980 The estimated initial state probabilities show that it is more probable to start from hidden state 1 than from hidden state 2 (66% vs. 34%). The high transition 396 J. Helske et al. Fig. 4 HMM with two hidden states (pie charts), with transitions between hidden states shown as labelled edges 0.091 0.039 0.66 0.34 Isolate Mediator Leader others probabilities on the diagonal of the transition matrix indicate that the students typically tend to stay in the hidden state they currently are in. Transition probabilities between the hidden states are relatively low and also asymmetric: it is more likely that students move from state 1 to state 2 than from state 2 to state 1. Looking at the emission matrices, we see that the role of the students in state 2 is mostly Leader or Mediator (emission probabilities are 50% and 48%). On the other hand, state 1 captures more of those occasions where students are isolated or exhibit at most a moderate level of participation (mediators). We can also visualise this with the plot() method of seqHMM (see Fig. 4): plot(fit_hmm$model, ncol.legend = 4, legend.prop = 0.2, edge.label.color = "black", vertex.label.color = "black" ) The plot values mainly shows the same information. By default, to simplify the graph, the plotting method combines all states with less than 5% emission probabilities into one category. This threshold can be changed with the combine.slices argument (setting combine.slices = 0 plots all states). For simple models, using n_states is sufficient. It automatically draws random starting values that are then used for the estimation of model parameters. However, as parameter estimation of HMMs and mixture models can be sensitive to starting values of parameters, it may be beneficial to provide starting values manually using the initial_probs,transition_probs, and emission_probs arguments. This is also necessary in case we want to define structural zeros for some of these components, e.g., if we want to restrict the initial probabilities so that each sequence starts from the same hidden state, or if we want to set the lower diagonal part of the transition matrix to zero, which means that the model does not allow transitioning back to previous states (this is called a left-to-right model) [4]. It is also possible to mix random and user-defined starting values by using simulate_*() functions (e.g. simulate_transition_probs()) for some of the model components and user-defined values for others. A Modern Approach to Transition Analysis and Process Mining with Markov... 397 In the following example we demonstrate estimating a three-state HMM with user-defined starting values for the initial state probabilities and the transition matrix but simulate starting values for the emission matrices. For simulating starting values with simulate_emission_probs(), we need to define the number of hidden states, and the number of observed symbols, i.e., the length of the alphabet of the sequences. # Set seed for randomisation set.seed(1) # Initial state probability vector, must sum to one init_probs <- c(0.3, 0.4, 0.3) # a 3x3 transition matrix, each row should sum to one trans_probs <- rbind( c(0.8, 0.15, 0.05), c(0.2, 0.6, 0.2), c(0.05, 0.15, 0.8) ) # Simulate emission probabilities emission_probs <- simulate_emission_probs( n_states = 3, n_symbols = length(alphabet(roles_seq)) ) # Build the HMM hidden_markov_model_2 <- build_hmm( roles_seq, initial_probs = init_probs, transition_probs = trans_probs, emission_probs = emission_probs ) Our initial probabilities suggest that it is slightly more likely to start from the second hidden state than the first and the third. Furthermore, the starting values for the transition matrices suggest that staying in hidden states 1 and 3 is more likely than staying in hidden state 2. All non-zero probabilities are, however, mere suggestions and will be estimated with the fit_model() function. We now estimate this model 50 times with the EM algorithm using randomised starting values: set.seed(1) fit_hmm_2 <- fit_model(hidden_markov_model_2, control_em = list(restart = list(times = 50)) ) 398 J. Helske et al. We can get the information on the EM estimation as follows: fit_hmm_2$em_results $logLik [1] -3546.155 $iterations [1] 488 $change [1] 9.947132e-11 $best_opt_restart [1] -3546.155 -3546.155 -3546.155 -3546.155 -3546.155 -3546.155 -3546.155 [8] -3546.155 -3546.155 -3546.155 -3546.155 -3546.155 -3546.155 -3546.155 [15] -3546.155 -3546.155 -3546.155 -3546.155 -3546.155 -3546.155 -3546.155 [22] -3546.155 -3546.155 -3546.155 -3546.155 The loglik element gives the log-likelihood of the final model. This value has no meaning on its own, but it can be used to compare HMMs with the same data and model structure (e.g., when estimating the same model from different starting values). The iterations and change arguments give information on the last EM estimation round: how many iterations were used until the (local) optimum was found and what was the change in the log-likelihood at the final step. The most interesting element is the last one: best_opt_restart shows the likelihood for 25 (by default) of the best estimation rounds. We advise to always check these to make sure that the best model was found several times from different starting values: this way we can be fairly certain that we have found the actual maximum likelihood estimates of the model parameters (global optimum). In this case all of the 25 log-likelihood values are identical, meaning that it is likely that we have found the best possible model among all HMMs with three hidden states. plot(fit_hmm_2$model, legend.prop = 0.15, ncol.legend = 3, edge.label.color = "black", vertex.label.color = "black", combine.slices = 0, trim = 0.0001 ) Interpreting the results in Fig. 5 we see that the first hidden state represents about equal amounts of isolate and mediator roles, the second hidden state represents mainly Leaders and some Mediator roles, and the third hidden state represents mainly Mediator roles and partly Leader roles. Interestingly, none of the students start as Mediator/Leader, while of the other two the Isolate/Mediator state is more typical (two thirds). There are no transitions from the first to the second state nor vice versa, and transition probabilities to the third state are considerably higher than away from it. In other words, it seems that the model has two different origin states and one destination state. A Modern Approach to Transition Analysis and Process Mining with Markov... 399 0.12 0.1 0.034 0.0077 0.69 0.31 0 Isolate Mediator Leader Fig. 5 HMM with three hidden states (pie charts), with transitions between hidden states shown as labelled edges 1 3 5 7 9 11 13 15 17 19 Observations Hidden states Isolate Mediator Leader State 1 State 2 State 3 n = 200 Fig. 6 Observed and hidden state sequences from the HMM with three hidden states 400 J. Helske et al. We can visualise the observed and/or hidden state sequences with the ssplot() function. The ssplot() function can take an stslist object or a model object of class mm or hmm (see Fig. 6). Here we want to plot full sequence index plots (type = "I") of both observed and hidden states (plots = "both") and sort the sequences using multidimensional scaling of hidden states (sortv = "mds.hidden"). See the seqHMM manual and visualisation vignette for more information on the different plotting options. ssplot(fit_hmm_2$model, # Plot sequence index plot (full sequences) type = "I", # Plot observed and hidden state sequences plots = "both", # Sort sequences by the scores of multidimensional scaling sortv = "mds.hidden", # X axis tick labels xtlab = 1:20 ) By looking at the sequences, we can see that even though none of the students start in hidden state 3, the majority of them transition there. In the end, most students end up alternating between mediating and leadership roles. Is the three-state model better than the two-state model? As already mentioned, we can use model selection criteria to test that. To make sure that the three-state model is the best, we also estimate a HMM with four hidden states and then use the Bayesian information criterion for comparing between the three models. Because the four-state model is more complex, we increase the number of re-estimation rounds for the EM algorithm to 100. # Set seed for randomisation set.seed(1) # Build and estimate a HMM with four states hidden_markov_model_3 <- build_hmm(roles_seq, n_states = 4) fit_hmm_3 <- fit_model(hidden_markov_model_3, control_em = list(restart = list(times = 100)) ) fit_hmm_3$em_results$best_opt_restart [1] -3534.304 -3534.304 -3534.304 -3534.304 -3534.304 -3534.304 -3534.304 [8] -3534.304 -3534.304 -3534.304 -3534.304 -3534.304 -3534.304 -3534.305 [15] -3534.305 -3534.306 -3534.308 -3534.310 -3534.332 -3534.335 -3534.335 [22] -3534.335 -3534.336 -3534.337 -3534.337 A Modern Approach to Transition Analysis and Process Mining with Markov... 401 The best model was found only 13 times out of 101 estimation rounds from randomised starting values. A cautious researcher might be wise to opt for a higher number of estimation rounds for increased certainty, but here we will proceed to calculating the BIC values. BIC(fit_hmm$model) [1] 7430.028 BIC(fit_hmm_2$model) [1] 7208.427 BIC(fit_hmm_3$model) [1] 7259.37 Generally speaking, the lower the BIC, the better the model. We can see that the three-state model (fit_hmm_2) has the lowest BIC value, so three clusters is the best choice (at least among HMMs with 2–4 hidden states). 4.2.3 Mixture Markov Models The MMM can be defined with the build_mmm() function. Similarly to HMMs, we need to either give the number of clusters with the n_clusters argument, which generates random starting values for the parameter estimates, or give starting values manually as initial_probs and transition_probs. Here we use random starting values: # Set seed for randomisation set.seed(123) # Define model structure (3 clusters) mmm <- build_mmm(roles_seq, n_clusters = 3) Again, the model is estimated with the fit_model() function: fit_mmm <- fit_model(mmm) The results for each cluster can be plotted one at a time (interactively, the default), or in one joint figure. Here we opt for the latter (see Fig. 7). At the same time we also illustrate some other plotting options: 402 J. Helske et al. Cluster 1 Isolate Mediator Leader Cluster 2 Isolate Mediator Leader Cluster 3 Isolate Mediator Leader Fig. 7 MMM with three clusters plot(fit_mmm$model, # Plot all clusters at the same time interactive = FALSE, # Set the number of rows (1) and columns (3) for cluster plots nrow = 1, ncol = 3, # Omit legends with.legend = FALSE, # Choose another layout for the vertices (see plot.igraph) layout = layout_in_circle, # Omit pie graphs from vertices pie = FALSE, # Set state colours vertex.label.color = c("black", "black", "white"), # Set state label colours vertex.color = cpal(roles_seq), # Increase the size of the circle vertex.size = 80, # Plot state labels instead of initial probabilities vertex.label = "names", # Set state label in the centre of the circle vertex.label.dist = 0, # Omit labels for transition probabilities edge.label = NA ) The following code plots the sequence distribution plot of each cluster (Fig. 8). In Cluster 1, we see low probabilities to downward mobility and high probabilities for upward mobility, so this cluster describes leadership trajectories. In Cluster 2, we can see that the thickest arrows lead to mediator and isolates roles, so this cluster describes trajectories with less central roles in collaboration. In Cluster 3, we see the highest transition probabilities for entering the mediator role but also A Modern Approach to Transition Analysis and Process Mining with Markov... 403 1 2 3 4 5 6 7 8 9 11 13 15 17 19 0.0 0.4 0.8 Proportion Isolate Mediator Leader Cluster 1, n = 48 1 2 3 4 5 6 7 8 9 11 13 15 17 19 0.0 0.4 0.8 Proportion Isolate Mediator Leader Cluster 2, n = 87 1 2 3 4 5 6 7 8 9 11 13 15 17 19 0.0 0.4 0.8 Proportion Isolate Mediator Leader Cluster 3, n = 65 Fig. 8 State distribution plots by most probable clusters estimated with the mixture Markov model. (a) Cluster 1. (b) Cluster 2. (c) Cluster 3 some transitions from mediator to leader, so this cluster describes trajectories with more moderate levels of participation in comparison to cluster 1. This behavior is easier to see when visualising the sequences in their most probable clusters. The plot is interactive, so we need to hit ‘Enter’ on the console to generate each plot. Alternatively, we can specify which cluster we want to plot using the which.plots argument. cl1 <- mssplot(fit_mmm$model, # Plot Y axis yaxis = TRUE, # Legend position with.legend = "bottom", # Legend columns ncol.legend = 3, # Label for Y axis ylab = "Proportion" ) We can add covariates to the model to explain cluster membership probabilities. For this, we need to provide a data frame (argument data) and the corresponding formula (argument formula). In the example data we use the data frame called cov_data that we created at the beginning of the tutorial with columns ID and GPA, 404 J. Helske et al. where the order of the ID variable matches to that of the sequence data roles_seq (note that the ID variable is not used in the model building, so the user needs to make sure that both matrices are sorted by ID). We can now use the information about students’ GPA level as a predictor of the cluster memberships. Numerical estimation of complex models from random starting values may lead to convergence issues and other problems in the estimation (you may, for example, get warnings about the EM algorithm failing). To avoid such issues, giving informative starting values is often helpful. This model is more complex than the model without covariates and estimation from random starting values leads to convergence issues (not shown here). To facilitate model estimation, we use the results from the previous MMM as informative starting values. Here we also remove the common intercept by adding 0 to the formula, which simplifies the interpretation of the covariate effects later (instead of comparing to a reference category, we get separate coefficients for each of the three GPA categories). set.seed(98765) mmm_2 <- build_mmm( roles_seq, # Starting values for initial probabilities initial_probs = fit_mmm$model$initial_probs, # Starting values for transition probabilities transition_probs = fit_mmm$model$transition_probs, # Data frame for covariates data = cov_data, # Formula for covariates (one-sided) formula = ~ 0 + GPA ) Again, the model is estimated with the fit_model() function. Here we use the EM algorithm with 50 restarts from random starting values: set.seed(12345) fit_mmm_2 <- fit_model( mmm_2, # EM with randomised restarts control_em = list( restart = list( # 50 restarts times = 50, # Store loglik values from all 50 + 1 estimation rounds n_optimum = 51 ) ) ) A Modern Approach to Transition Analysis and Process Mining with Markov... 411 legend.prop = 0.4, ncol.legend = 1, ncol = 3, interactive = FALSE, combine.slices = 0 ) Based on the two plots, we can determine that Cluster 1 describes students who start as leaders but then transition to alternating between mediator and leader. Cluster 2 describes students who start by alternating between isolate and mediator roles and then mainly transition to alternating between mediator and leader roles. Cluster 3 describes students who start as alternating between isolate and mediator roles, after which they transition between isolate/mediator and mediator/leader. cluster_names(fit_mhmm_2$model) <- c( "Downward transition", "Upward transition", "Alternating" ) With summary(fit_mhmm_2$model) we get the parameter estimates and standard errors for the covariates and information about clustering: summary(fit_mhmm_2$model) Covariate effects : Downward transition is the reference. Upward transition : Estimate Std. error GPALow -0.455 0.464 GPAMiddle 0.440 0.310 GPAHigh -2.743 0.727 Alternating : Estimate Std. error GPALow 1.3560 0.324 GPAMiddle 0.3461 0.316 GPAHigh 0.0468 0.250 Log-likelihood: -3519.243 BIC: 7237.543 Means of prior cluster probabilities : Downward transition Upward transition Alternating 0.302 0.181 0.517 Most probable clusters : Downward transition Upward transition Alternating count 61 30 109 proportion 0.305 0.15 0.545 412 J. Helske et al. Classification table : Mean cluster probabilities (in columns) by the most probable cluster (rows) Downward transition Upward transition Alternating Downward transition 0.95727 0.0267 0.0161 Upward transition 0.03007 0.8037 0.1662 Alternating 0.00975 0.0962 0.8940 We can see, that the prior probabilities of belonging to each cluster are very different: half of the students can be described as alternating, while of the rest, a downward transition is more typical (31%). Based on the classification table, the Downward transition cluster is rather crisp, while the other two are partly overlapping (see the MMM example for more information on interpreting the classification table). The Covariate effects tables show that, in comparison to Alternating cluster, students with low GPA are less likely to end up in the Upward or Downward transition clusters and students with high GPA are less likely to end up in Upward transition cluster. Again, we can calculate the probabilities of belonging to each cluster by GPA levels: exp(fit_mhmm_2$model$coefficients)/rowSums(exp(fit_mhmm_2$model$coefficients)) Downward transition Upward transition Alternating GPALow 0.1813217 0.11502283 0.7036555 GPAMiddle 0.2521406 0.39144399 0.3564154 GPAHigh 0.4734128 0.03048189 0.4961054 The table shows that students with low GPA typically belong to the Alternating cluster (70% probability) while students with high GPA mainly end up in the Downward transition cluster (47%) or the Alternating cluster (50%). Most students with middle GPA end up in the Upward transition cluster (39%), but the probabilities are almost as high for the Alternating cluster (36%) and also fairly high for the Downward transition cluster (25%). In light of this, it is worth noting that the covariates do not merely explain the uncovered clusters; as part of the model, they drive the formation of the clusters. In other words, an otherwise identical model without the dependence on the GPA covariate may uncover different groupings with different probabilities. If we are not sure how many clusters or hidden states we expect, or if we wish to investigate different combinations of covariates, we can estimate several models and compare the results with information criteria or cross-validation. Estimating a large number of complex models is, however, very time-consuming. Using prior information for restricting the pool of potential models is useful, and sequence analysis can also be used as a helpful first step [10, 37]. A Modern Approach to Transition Analysis and Process Mining with Markov... 413 4.3 Stochastic Process Mining with Markovian Models Process mining is a relatively recent method for the analysis of event-log data (timestamped logs) which aims to understand the flow and dynamics of the process under study. In education, process mining has been used extensively to analyse learners’ online logs collected from Learning Management Systems (LMS), to understand how they utilize learning resources and transitions between learning activities to mention a few [38, 39]. In this book, we have devoted a full chapter for process mining where we explained how process mining can be performed in R [40]. Yet, in this chapter we will present a novel method that we propose to perform stochastic process mining using MMs. While process mining can be performed using different software, techniques and algorithms, MMs offer a powerful framework for process mining with several advantages over the commonly used methods. First, it is more theoretically aligned with the idea of a transition from an action to an action and that actions are temporally dependent on each other. Second, MMs allow for data to be clustered into similar transition patterns, a possibility not offered by other process mining methods (see the process mining chapter of this book [40]). Third, contrary to other process mining methods, MMs do not require researchers to arbitrarily exclude—or trim—a large part of the data to “simplify” the model. For instance, most of the process mining analyses require an arbitrary cutoff to trim the data so that the process model is readable. This trimming significantly affects the resulting model and makes it hard to replicate. Most importantly, MMs have several fit statistics that we can use to compare and judge the model fit as we have seen before. Several R packages can perform stochastic process mining; in this tutorial we will rely on the same package we discussed earlier and combine it with a powerful visualisation that allows us to effectively visualise complex processes. In the next example, we will analyse data extracted from the learning management system logs and offer a detailed guide to process mining. We will also use MMMs to cluster the data into latent patterns of transitions. Given that the traditional plotting function in seqHMM works well with a relatively short alphabet, we will use a new R package called qgraph for plotting. The package qgraph offers powerful visualisations which makes plotting easier, and more interpretable especially for larger models. Furthermore, qgraph allows researchers to use a fixed layout for all the plotted networks so the nodes can be compared to each other more easily. Let us now go through the analysis. The next chunk of code imports the prepared sequence data from the sequence analysis chapter. The data belong to a learning analytics course and the events are coded trace logs of students’ actions such as Course view, Instructions, Practicals, Social, etc. Then, we build a sequence object using the function seqdef() from TraMineR. 414 J. Helske et al. library(qgraph) library(rio) library(seqHMM) library(tidyverse) library(TraMineR) seq_data <- import(paste0(URL, "1_moodleLAcourse/LMS_data_wide.xlsx")) seq_data_all <- seqdef(seq_data, var = 7:54 ) Before proceeding further, it is advisable to visualise the sequences. Figure 12 shows the sequence index plot, sorted according to the first states. The data are much larger than the collaboration roles and achievement sequences analysed previously; there are 9478 observations with an alphabet of 12 states. Unlike in the previous example, the sequence lengths vary considerably. Due to this, shorter sequences contain missing values to fill the empty cells in the data frame. However, there Fig. 12 Sequence index plot for the learning management system logs A Modern Approach to Transition Analysis and Process Mining with Markov... 415 are no internal gaps. When creating the sequence object with the seqdef function, TraMineR allows for distinguishing between real missing values (NA, where the true state is unknown) and technical missing values (void) used to pad the sequences to equal lengths. The seqHMM package is able to account for both types of missing values and treats them slightly differently, for example when calculating the most probable paths of hidden states. seqplot(seq_data_all, type = "I", ncol = 4, sortv = "from.start", legend.prop = 0.2, cex.legend = 0.7, border = NA, ylab = "Sequence (sorted)", xlab = "Time" ) A simple transition analysis can be performed by estimating and plotting the transition probabilities. This can be performed using the TraMineR package. Yet, this simple approach has drawbacks and it is advisable to estimate the MM and use their full power. The next code estimates the transition probabilities of the full dataset using the function seqtrate() from TraMineR package (shown in Table 7). overalltransitions <- seqtrate(seq_data_all) As we mentioned earlier, we will use a novel plotting technique that is more suitable for large process models. Below, we plot the transition probabilities with the qgraph() function from the qgraph package (Fig. 13). We use some arguments to improve the process model visualisation. First, we use the argument cut = 0.15 to show the edges with probabilities below 0.15 in lower thickness and colour intensity. This cut makes the graph easier to read and less crowded, and gives emphasis to the edges which matter. The argument minimum = 0.05 hides small edges below the probability threshold of 0.05. We use edge.labels = TRUE to show the transition probabilities as edge labels. The argument color gets the colour palette from the sequence with the function cpal() and the argument curveAll = TRUE ensures the graph shows curved edges. The “colorblind” theme makes sure that the colours can be seen by everyone regardless of colour vision abilities. Lastly, the mar within the figure area. # get the labeles to use them as node names. Labelx <- alphabet(seq_data_all) transitionsplot <- qgraph( overalltransitions, cut = 0.15, minimum = 0.05, labels = Labelx, edge.labels = TRUE, edge.label.cex = 0.65, color = cpal(seq_data_all), curveAll = TRUE, theme = "colorblind", mar = c(4, 3, 4, 3) ) 416 J. Helske et al. Table 7 Transition probabilities From\To Applications Assignment Course_view Ethics Feedback General Group_work Instructions La_types Practicals Social Theory Applications .0.46 .0.07 .0.13 .0.01 .0.01 .0.19 .0.05 .0.01 .0.01 .0.05 .0.00 . 0.00 Assignment .0.00 .0.70 .0.19 .0.00 .0.01 .0.02 .0.03 .0.02 .0.02 .0.02 .0.00 . 0.00 Course_view .0.01 .0.07 .0.35 .0.01 .0.03 .0.03 .0.28 .0.10 .0.02 .0.08 .0.02 . 0.01 Ethics .0.01 .0.00 .0.12 .0.61 .0.01 .0.04 .0.10 .0.01 .0.03 .0.04 .0.01 . 0.02 Feedback .0.00 .0.02 .0.23 .0.00 .0.56 .0.00 .0.11 .0.04 .0.01 .0.02 .0.00 . 0.00 General .0.04 .0.05 .0.18 .0.01 .0.00 .0.49 .0.06 .0.06 .0.05 .0.03 .0.01 . 0.02 Group_work .0.00 .0.01 .0.19 .0.00 .0.01 .0.01 .0.73 .0.02 .0.00 .0.01 .0.01 . 0.00 Instructions .0.00 .0.02 .0.33 .0.00 .0.03 .0.04 .0.12 .0.37 .0.02 .0.03 .0.04 . 0.00 La_types .0.01 .0.06 .0.24 .0.01 .0.00 .0.10 .0.07 .0.05 .0.38 .0.03 .0.01 . 0.03 Practicals .0.00 .0.02 .0.17 .0.00 .0.01 .0.01 .0.03 .0.02 .0.01 .0.73 .0.00 . 0.01 Social .0.00 .0.01 .0.25 .0.00 .0.00 .0.01 .0.12 .0.11 .0.01 .0.02 .0.48 . 0.00 Theory .0.00 .0.02 .0.15 .0.03 .0.00 .0.02 .0.06 .0.01 .0.05 .0.05 .0.00 .0.60 A Modern Approach to Transition Analysis and Process Mining with Markov... 417 0.05 0.05 0.05 0.05 0.06 0.06 0.06 0.06 0.07 0.07 0.07 0.08 0.1 0.1 0.1 0.11 0.11 0.12 0.12 0.12 0.13 0.15 0.17 0.18 0.19 0.19 0.19 0.23 0.24 0.25 0.28 0.33 0.35 0.37 0.38 0.46 0.48 0.49 0.56 0.6 0.61 0.7 0.73 0.73 Applications Assignment Course_view Ethics Feedback General Group_work Instructions La_types Practicals Social Theory Fig. 13 Process map for the overall process The seqtrate() function only computes the transition probabilities but does not compute the initial probabilities. While it is not difficult to calculate the proportions of starting in each state, we can also estimate a simple Markov model which does the same with a short command. We do so using the build_mm() function as per Sect. 4.2, recalling that the build_mm() function is distinct from build_hmm(), build_mmm(), and build_mhmm() in that it is the only build function that automatically estimates the parameters of the model. The plotting now includes an extra option called pie = overallmodel$initi al_probs which tells qgraph to use the initial probabilities from the fitted MM as the sizes of the pie charts in the borders of the nodes in Fig. 14. For instance, the pie around Course view is around half of the circle corresponding to 0.48 initial probability to start from Course view. Please also note that the graph is otherwise equal to the one generated via seqtrate() apart from these initial probabilities. 418 J. Helske et al. 0.05 0.05 0.05 0.05 0.06 0.06 0.06 0.06 0.07 0.07 0.07 0.08 0.1 0.1 0.1 0.11 0.11 0.12 0.12 0.12 0.13 0.15 0.17 0.18 0.19 0.19 0.19 0.23 0.24 0.25 0.28 0.33 0.35 0.37 0.38 0.46 0.48 0.49 0.56 0.6 0.61 0.7 0.73 0.73 Applications Assignment Course_view Ethics Feedback General Group_work Instructions La_types Practicals Social Theory Fig. 14 Process map for the overall process with initial probabilities overallmodel <- build_mm(seq_data_all) overallplot <- qgraph( overallmodel$transition_probs cut = 0.15, minimum = 0.05, labels = Labelx, mar = c(4, 3, 4, 3), edge.labels = TRUE, edge.label.cex = 0.65, color = cpal(seq_data_all), curveAll = TRUE, theme = "colorblind", pie = overallmodel$initial_probs ) A Modern Approach to Transition Analysis and Process Mining with Markov... 419 Having plotted the transitions of the full dataset, we can now look for transition patterns, that is typical transition patterns (i.e., clusters) that are repeated within the data. The procedure is the same as before. In the next example, we use the function build_mmm() to build the model with four clusters as a demonstration. Ideally, researchers need to estimate several models and choose the best model based on model selection criteria (such as BIC) values as well as interpretability. The steps involved in fitting the model are as before; we make use of the function fit_model() to estimate the model. The results of the running the code will be an MM for each cluster (with distinct initial and transition probabilities). Given the number of sequences in the dataset, their length, and the number of states, the computational burden is larger than for previous applications in this chapter. For illustrative purposes, instead of repeated EM runs with random starting values, we use single EM run followed by global optimisation, using the argument global_step = TRUE. One benefit of this global (and local) step in fit_model over the EM algorithm is the flexibility to define a maximum runtime (in seconds) for the optimization process (argument maxtime in control_global). This can be valuable for larger problems with predefined runtime (e.g., in a shared computer cluster). Note, however, that relying on the runtime can lead to non-reproducible results even with fixed seed if the optimisation terminates due to the time limit. Finally, we run additional local optimisation step using the results of the global optimisation, for more accurate results. The last argument threads = 16 instructs to use parallel computing to enable faster fitting (please, customise according to the number of cores in your computer). As for the starting values, we use the transition probabilities computed from the full data for all clusters, and random values for the initial probabilities. While in theory many of the global optimisation algorithms should eventually find the global optimum, in practice there are no guarantees that it is found in limited time. Thus, as earlier, in practice it is advisable to try different global/local optimisation algorithms and/or EM algorithm with different initial values to make it more likely that the global optimum is found (see [4] for further discussion). set.seed(1) trans_probs <- simulate_transition_probs(12, 4, diag_c = 5) init_probs <- as.numeric(prop.table(table(seq_data_all[,1])[1:12])) init_probs <- replicate(4, init_probs, simplify = FALSE) builtseqLMS <- build_mmm( seq_data_all, transition_probs = trans_probs, initial_probs = init_probs ) fitLMS <- fit_model( builtseqLMS, global_step = TRUE, control_global = list( maxtime = 3600, 420 J. Helske et al. maxeval = 1e5, algorithm = "NLOPT_GD_STOGO_RAND"), local_step = TRUE, threads = 16 ) fitLMS$global_results$message fitLMS$logLik [1] "NLOPT_SUCCESS: Generic success return value." [1] -114491.2 Before plotting the clusters, let us do some cleanups. First, we get the transition probabilities of each cluster and assign them to a variable. In that way, they are easier to manipulate and work with. In the same way, we can extract the initial probabilities for each cluster. #extract transition probabilities of each cluster Clustertp1 <- fitLMS$model$transition_probs$`Cluster 1` Clustertp2 <- fitLMS$model$transition_probs$`Cluster 2` Clustertp3 <- fitLMS$model$transition_probs$`Cluster 3` Clustertp4 <- fitLMS$model$transition_probs$`Cluster 4` #extract initial probabilities of each cluster Clusterinitp1 <- fitLMS$model$initial_probs$`Cluster 1` Clusterinitp2 <- fitLMS$model$initial_probs$`Cluster 2` Clusterinitp3 <- fitLMS$model$initial_probs$`Cluster 3` Clusterinitp4 <- fitLMS$model$initial_probs$`Cluster 4` Plotting the process maps can be performed in the same way we did before. However, if we need to compare clusters, it is best if we use a unified layout. An average layout can be computed with the function averageLayout() which takes the transition probabilities of the four clusters as input and creates—as the name implies—an averaged layout. Another option is to use the same layout of the overallplot in the previous example. This can be obtained from the plot object overallplot$layout. This can be helpful if you would like to plot the four plots corresponding to each cluster with the same layout as the overall plot (see Fig. 15). Labelx <- colnames(Clustertp1) # we need to get the labels Averagelayout <- averageLayout( list(Clustertp1, Clustertp2, Clustertp3, Clustertp4) ) A Modern Approach to Transition Analysis and Process Mining with Markov... 427 33. Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019) Welcome to the tidyverse. J Open Source Softw 4:1686. https://doi.org/10.21105/joss.01686 34. Gabadinho A, Ritschard G, Müller NS, Studer M (2011) Analyzing and visualizing state sequences in R with TraMineR. J Stat Softw 40. https://doi.org/10.18637/jss.v040.i04 35. Saqr M, López-Pernas-Pernas S, Helske S, Hrastinski S (2023) The longitudinal association between engagement and achievement varies by time, students’ profiles, and achievement state: a full program study. Comput Educ 199:104787 36. Saqr M, López-Pernas-Pernas S (2022) How CSCL roles emerge, persist, transition, and evolve over time: a four-year longitudinal study. Comput Educ 189:104581. https://doi.org/10.1016/j. compedu.2022.104581 37. Helske S, Keski-Säntti M, Kivelä J, Juutinen A, Käriälä A, Gissler M, Merikukka M, Lallukka T (2023) Predicting the stability of early employment with its timing and childhood social and health-related predictors: a mixture markov model approach. Longitud Life Course Stud 14:73–104 38. Peeters W, Saqr M, Viberg O (2020) Applying learning analytics to map students’ selfregulated learning tactics in an academic writing course. In: Proceedings of the 28th international conference on computers in education. Asia-Pacific Society for Computers in Education, pp 245–254 39. Saqr M, Matcha W, Jovanovic J, Gaševi´ c D, López-Pernas-Pernas S, et al (2022) Transferring effective learning strategies across learning contexts matters: a study in problem-based learning. Australas J Educ Technol 39(3)35–57 40. López-Pernas-Pernas S, Saqr M (2024) The why, the how, and the when of educational process mining in R. In: Saqr M, López-Pernas-Pernas S (eds) Learning analytics methods and tutorials: a practical guide using R, Chap. 14. Springer, Cham 41. Tikka S, Helske J (2023) dynamite: an R package for dynamic multivariate panel models. https://doi.org/10.48550/ARXIV.2302.01607 42. Bartolucci F, Pandolfi S, Pennoni F (2017) LMest: an R package for latent Markov models for longitudinal categorical data. J Stat Softw 81:1–38. https://doi.org/10.18637/jss.v081.i04 43. Vermunt JK, Magidson J (2016) Guide for latent GOLD 5.1: basic, advanced, and syntax. Statistical Innovations Inc., Belmont 44. Berchtold A (1999) The double chain Markov model. Commun Stat Theory Methods 28:2569– 2589. https://doi.org/10.1080/03610929908832439 45. Maitre O, Emery K, Oliver Buschor with contributions from, Berchtold A (2020). march: Markov chains. https://CRAN.R-project.org/package=march 46. Gabadinho A, Ritschard G (2016) Analyzing state sequences with probabilistic suffix trees: the PST R package. J Stat Softw 72:1–39. https://doi.org/10.18637/jss.v072.i03 Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made. The images or other third party material in this chapter are included in the chapter’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.