Mixture Hidden Markov Models for Sequence Data : The seqHMM Package in R

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: Creative Commons Attribution 3.0 Unported License https://creativecommons.org/licenses/by/3.0/ Mixture Hidden Markov Models for Sequence Data : The seqHMM Package in R © 2019 the Authors Published version Helske, Satu; Helske, Jouni Helske, S., & Helske, J. (2019). Mixture Hidden Markov Models for Sequence Data : The seqHMM Package in R. Journal of Statistical Software, 88, Article 3. https://doi.org/10.18637/jss.v088.i03 2019 JSS Journal of Statistical Software January 2019, Volume 88, Issue 3. doi: 10.18637/jss.v088.i03 Mixture Hidden Markov Models for Sequence Data: The seqHMM Package in R Satu Helske Linköping University University of Oxford University of Jyväskylä Jouni Helske Linköping University University of Jyväskylä Abstract Sequence analysis is being more and more widely used for the analysis of social sequences and other multivariate categorical time series data. However, it is often complex to describe, visualize, and compare large sequence data, especially when there are multiple parallel sequences per subject. Hidden (latent) Markov models (HMMs) are able to detect underlying latent structures and they can be used in various longitudinal settings: to account for measurement error, to detect unobservable states, or to compress information across several types of observations. Extending to mixture hidden Markov models (MHMMs) allows clustering data into homogeneous subsets, with or without external covariates. The seqHMM package in Ris designed for the efficient modeling of sequences and other categorical time series data containing one or multiple subjects with one or multiple interdependent sequences using HMMs and MHMMs. Also other restricted variants of the MHMM can be fitted, e.g., latent class models, Markov models, mixture Markov models, or even ordinary multinomial regression models with suitable parameterization of the HMM. Good graphical presentations of data and models are useful during the whole analysis process from the first glimpse at the data to model fitting and presentation of results. The package provides easy options for plotting parallel sequence data, and proposes visualizing HMMs as directed graphs. Keywords: multi-channel sequences, categorical time series, visualizing sequence data, visualizing models, latent Markov models, latent class models, R. 1. Introduction Social sequence analysis is being more and more widely used for the analysis of longitudinal data consisting of multiple independent subjects with one or multiple interdependent 2seqHMM: Mixture Hidden Markov Models for Sequence Data in R sequences (channels). Sequence analysis is used for computing the (dis)similarities of sequences, and often the goal is to find patterns in data using cluster analysis. However, describing, visualizing, and comparing large sequence data is often complex, especially in the case of multiple channels. Hidden (latent) Markov models (HMMs) can be used to compress and visualize information in such data. These models are able to detect underlying latent structures. Extending to mixture hidden Markov models (MHMMs) allows clustering via latent classes, possibly with additional covariate information. One of the major benefits of using hidden Markov modeling is that all stages of analysis are performed, evaluated, and compared in a probabilistic framework. The seqHMM package (Helske and Helske 2019) for R(RCore Team 2018) is designed for modeling sequence data and other categorical time series with one or multiple subjects and one or multiple channels using HMMs and MHMMs. The package provides functions for the estimation and inference of models, as well as functions for the easy visualization of multi-channel sequences and HMMs. Even though the package was originally developed for researchers familiar with social sequence analysis and the examples are related to life course, knowledge on sequence analysis or social sciences is not necessary for the usage of seqHMM. The package is available from the Comprehensive RArchive Repository (CRAN) at https://CRAN.R-project.org/package=seqHMM and under development on GitHub at https://github.com/helske/seqHMM. There are also other Rpackages on CRAN for HMM analysis of categorical data. The HMM package (Himmelmann 2010) is a compact package designed for fitting an HMM for a single observation sequence. The hmm.discnp package (Turner 2018) can handle multiple observation sequences with possibly varying lengths. For modeling continuous-time processes as hidden Markov models, the msm package (Jackson 2011) is available. Both hmm.discnp and msm support only single-channel observations. The depmixS4 package (Visser and Speekenbrink 2010) is able to fit HMMs for multiple interdependent time series (with continuous or categorical values), but for one subject only. In the msm and depmixS4 packages, covariates can be added for initial and transition probabilities. The mhsmm package (O’Connell and Højsgaard 2011) allows modeling of multiple sequences using hidden Markov and semi-Markov models. There are no ready-made options for modeling categorical data, but users can write their own extensions for arbitrary distributions. The LMest package (Bartolucci, Pandolfi, and Pennoni 2017) is aimed at panel data with a large number of subjects and a small number of time points. It can be used for hidden Markov modeling of multivariate and multi-channel categorical data, using covariates in emission and transition processes. LMest also supports mixed latent Markov models, where the latent process is allowed to vary in different latent subpopulations. This differs from mixture hidden Markov models used in seqHMM, where also the emission probabilities vary between groups. The seqHMM package also supports covariates in explaining group memberships. A drawback in the LMest package is that the user cannot define initial values or zero constraints for model parameters, and thus important special cases such as left-to-right models cannot be used. We start with describing data and methods: a short introduction to sequence data and sequence analysis, then the theory of hidden Markov models for such data, an expansion to mixture hidden Markov models and a glance at some special cases, and then some propositions on visualizing multi-channel sequence data and hidden Markov models. After the theoretic part we take a look at features of the seqHMM package and at the end show an example on using the package for the analysis of life course data. Appendix Ashows the list of notations. Journal of Statistical Software 3 2. Methods 2.1. Sequences and sequence analysis By the term sequence we refer to an ordered set of categorical states. It can be a time series, such as a career trajectory or residential history, or any other series with ordered categorical observations, e.g., a DNA sequence or a structure of a story. Typically, sequence data consist of multiple independent subjects (multivariate data). Sometimes there are also multiple interdependent sequences per subject, often referred to as multi-channel or multidimensional sequence data. As an example we use the biofam data available in the TraMineR package (Gabadinho, Ritschard, Müller, and Studer 2011). It is a sample of 2000 individuals born in 1909–1972, constructed from the Swiss Household Panel survey in 2002 (Müller, Studer, and Ritschard 2007). The data set contains sequences of annual family life statuses from age 15 to 30. Eight observed states are defined from the combination of five basic states: living with parents, left home, married, having children, and divorced. To show a more complex example, we split the original data into three separate channels representing different life domains: marriage, parenthood, and residence. The data for each individual now includes three parallel sequences constituting of two or three states each: single/married/divorced, childless/parent, and living with parents/having left home. Sequence analysis (SA), as defined in the social science framework, is a model-free data-driven approach to the analysis of successions of states. The approach has roots in bioinformatics and computer science (see, e.g., Durbin, Eddy, Krogh, and Mitchison 1998), but during the past few decades SA has also become more common in other disciplines for the analysis of longitudinal data. In social sciences SA has been used increasingly often and is now “central to the life-course perspective” (Blanchard, Bühlmann, and Gauthier 2014). SA is used for computing (dis)similarities of sequences. The most well-known method is optimal matching (McVicar and Anyadike-Danes 2002), but several alternatives exist (see, e.g., Aisenbrey and Fasang 2010;Elzinga and Studer 2014;Gauthier, Widmer, Bucher, and Notredame 2009;Halpin 2010;Hollister 2009;Lesnard 2010). Also a method for analyzing multi-channel data has been developed (Gauthier, Widmer, Bucher, and Notredame 2010). Often the goal in SA is to find typical and atypical patterns in trajectories using cluster analysis, but any approach suitable for compressing information on the dissimilarities can be used. The data are usually presented also graphically in some way. So far the TraMineR package has been the most extensive and frequently used software for social sequence analysis. 2.2. Hidden Markov models In the context of hidden Markov models, sequence data consists of observed states, which are regarded as probabilistic functions of hidden states. Hidden states cannot be observed directly, but only through the sequence(s) of observations, since they emit the observations on varying probabilities. A discrete first order hidden Markov model for a single sequence is characterized by the following: •Observed state sequence y= (y1, y2, . . . , yT)with observed states m∈ {1, . . . , M}. •Hidden state sequence z= (z1, z2, . . . , zT)with hidden states s∈ {1, . . . , S}. 4seqHMM: Mixture Hidden Markov Models for Sequence Data in R •Initial probability vector π={πs}of length S, where πsis the probability of starting from the hidden state s: πs=P(z1=s); s∈ {1, . . . , S}. •Transition matrix A={asr}of size S×S, where asr is the probability of moving from the hidden state sat time t−1to the hidden state rat time t: asr =P(zt=r|zt−1=s); s, r ∈ {1, . . . , S}. We only consider homogeneous HMMs, where the transition probabilities asr are constant over time. •Emission matrix B={bs(m)}of size S×M, where bs(m)is the probability of the hidden state semitting the observed state m: bs(m) = P(yt=m|zt=s); s∈ {1, . . . , S}, m ∈ {1, . . . , M}. The (first order) Markov assumption states that the hidden state transition probability at time tonly depends on the hidden state at the previous time point t−1: P(zt|zt−1, . . . , z1) = P(zt|zt−1).(1) Also, the observation at time tis only dependent on the current hidden state, not on previous hidden states or observations: P(yt|yt−1, . . . , y1, zt, . . . , z1) = P(yt|zt).(2) For a more detailed description of hidden Markov models, see, e.g., Rabiner (1989), MacDonald and Zucchini (1997), and Durbin et al. (1998). HMM for multiple sequences We can also fit the same HMM for multiple subjects; instead of one observed sequence ywe have Nsequences as Y= (y1,...,yN)>, where the observations yi= (yi1, . . . , yiT )of each subject itake values in the observed state space. Observed sequences are assumed to be mutually independent given the hidden states. The observations are assumed to be generated by the same model, but each subject has its own hidden state sequence. HMM for multi-channel sequences In the case of multi-channel sequence data, such as the example described in Section 2.1, for each subject ithere are Cparallel sequences. Observations are now of the form yitc, i= 1, . . . , N;t= 1 ...,T;c= 1 ...,C, so that our complete data is Y={Y1, . . . , Y C}. In seqHMM, multi-channel data are handled as a list of Cdata frames of size N×T. We also define Yias all the observations corresponding to subject i. We apply the same latent structure for all channels. In such a case the model has one transition matrix Abut several emission matrices B1, . . . , BC, one for each channel. We assume that the observed states in different channels at a given time point tare independent of each other given the hidden state at t, i.e., P(yit|zit) = P(yit1|zit)· · · P(yitC|zit). Journal of Statistical Software 5 Sometimes the independence assumption does not seem theoretically plausible. For example, even conditioning on a hidden state representing a general life stage, are marital status and parenthood truly independent? On the other hand, given a person’s religious views, could their opinions on abortion and gay marriage be though as independent? If the goal is to use hidden Markov models for prediction or simulating new sequence data, the analyst should carefully check the validity of independence assumptions. However, if the goal is merely to describe structures and compress information, it can be useful to accept the independence assumption even though it is not completely reasonable in a theoretical sense. When using multi-channel sequences, the number of observed states is smaller, which leads to a more parsimonious representation of the model and easier inference of the phenomenon. Also due to the decreased number of observed states, the number of parameters of the model is decreased leading to the improved computational efficiency of model estimation. The multi-channel approach is particularly useful if some of the channels are only partially observed; combining missing and non-missing information into one observation is usually problematic. One would have to decide whether such observations are coded completely missing, which is simple but loses information, or whether all possible combinations of missing and non-missing states are included, which increases the state space and makes the interpretation of the model more difficult. In the multi-channel approach the data can be used as it is. Missing data Missing observations are handled straightforwardly in the context of HMMs. When observation yitc is missing, we gain no additional information regarding hidden states. In such a case, we set the emission probability bs(yitc)=1for all s∈ {1, . . . , S}. Sequences with varying lengths are handled by setting missing values before and/or after the observed states. Log-likelihood and parameter estimation The unknown transition, emission and initial probabilities are commonly estimated via maximum likelihood. The log-likelihood of the parameters M={π, A, B1, . . . , BC}for the HMM is written as log L= N X i=1 log P(Yi|M),(3) where Yiare the observed sequences in channels c= 1, . . . , C for subject i. The probability of the observation sequence of subject igiven the model parameters is P(Yi|M) = X all z P(Yi|z, M)P(z|M) =X all z P(z1|M)P(yi1|z1,M) T Y t=2 P(zt|zt−1,M)P(yit|zt,M) =X all z πz1bz1(yi11)· · · bz1(yi1C) T Y t=2 azt−1ztbzt(yit1)· · · bzt(yitC), (4) where the hidden state sequences z= (z1, . . . , zT)take all possible combinations of values in the hidden state space {1, . . . , S}and where yit are the observations of subject iat tin channels 1, . . . , C;πz1is the initial probability of the hidden state at time t= 1 in sequence 6seqHMM: Mixture Hidden Markov Models for Sequence Data in R z;azt−1ztis the transition probability from the hidden state at time t−1to the hidden state at t; and bzt(yitc)is the probability that the hidden state of subject iat time temits the observed state at tin channel c. For direct numerical maximization (DNM) of the log-likelihood, any general-purpose optimization routines such as BFGS or Nelder-Mead can be used (with suitable reparameterizations). Another common estimation method is the expectation-maximization (EM) algorithm, also known as the Baum-Welch algorithm in the HMM context. The EM algorithm rapidly converges close to a local optimum, but compared to DNM, the converge speed is often slow near the optimum. The probability (4) is efficiently calculated using the forward part of the forward-backward algorithm (Baum and Petrie 1966;Rabiner 1989). The backward part of the algorithm is needed for the EM algorithm, as well as for the computation of analytic gradients for derivative based optimization routines. For more information on the algorithms, see a supplementary vignette for package seqHMM on CRAN (Helske 2018a). The estimation process starts by giving initial values to the estimates. Good starting values are needed for finding the optimal solution in a reasonable time. In order to reduce the risk of being trapped in a poor local maximum, a large number of initial values should be tested. Inference on hidden states Given our model and observed sequences, we can make several interesting inferences regarding the hidden states. Forward probabilities αit(s)(Rabiner 1989) are defined as the joint probability of hidden state sat time tand the observation sequences yi1,...,yit given the model M, whereas backward probabilities βit(s)are defined as the joint probability of hidden state sat time tand the observation sequences yi(t+1),...,yiT given the model M. From forward and backward probabilities we can compute the posterior probabilities of states, which give the probability of being in each hidden state at each time point, given the observed sequences of subject i. These are defined as P(zit =s|Yi,M) = αitβit P(Yi|M).(5) Posterior probabilities can be used to find the locally most probable hidden state at each time point, but the resulting sequence is not necessarily globally optimal. To find the single best hidden state sequence ˆzi(Yi) = {ˆzi1,ˆzi2,...,ˆziT }for subject i, we maximize P(z|Yi,M)or, equivalently, P(z, Yi|M). A dynamic programming method, the Viterbi algorithm (Rabiner 1989), is used for solving the problem. Model comparison Models with the same number of parameters can be compared with the log-likelihood. For choosing between models with a different number of hidden states, we need to take account of the number of parameters. We define the Bayesian information criterion (BIC) as BIC =−2 log(Ld) + plog N X i=1 T X t=1 1 C C X c=1 I(yitc observed)!,(6) where Ldis computed using Equation 3,pis the number of estimated parameters, I is the indicator function, and the summation in the logarithm is the size of the data. If data are Journal of Statistical Software 7 completely observed, the summation is simplified to N×T. Missing observations in multichannel data may lead to non-integer data size. 2.3. Clustering by mixture hidden Markov models There are many approaches for finding and describing clusters or latent classes when working with HMMs. A simple option is to group sequences beforehand (e.g., using sequence analysis and a clustering method); afterwards one HMM is fitted for each cluster. This approach is simple in terms of HMMs. Models with a different number of hidden states and initial values are explored and compared for one cluster at a time. HMMs are used for compressing information and comparing different clustering solutions, e.g., finding the best number of clusters. The problem with this solution is that it is, of course, very sensitive to the original clustering and the estimated HMMs might not be well suited for borderline cases. Instead of fixing sequences into clusters, it is possible to fit one model for the whole data and determine clustering during modeling. Now sequences are not in fixed clusters but get assigned to clusters with certain probabilities during the modeling process. In this section we expand the idea of HMMs to mixture hidden Markov models (MHMMs). This approach was formulated by van de Pol and Langeheine (1990) as a mixed Markov latent class model and later generalized to include time-constant and time-varying covariates by Vermunt, Tran, and Magidson (2008) (who named the resulting model the mixture latent Markov model; MLMM). The MHMM presented here is a variant of MLMM where only time-constant covariates are allowed. Time-constant covariates deal with unobserved heterogeneity and they are used for predicting cluster memberships of subjects. Mixture hidden Markov model Assume that we have a set of HMMs M={M1,...,MK}, where Mk={πk, Ak, Bk 1, . . . , Bk C} for submodels k= 1, . . . , K. For each subject Yi, denote P(Mk) = wkas the prior probability that the observation sequences of a subject follow the submodel Mk. Now the log-likelihood of the parameters of the MHMM is extended from Equation 3to log L= N X i=1 log P(Yi|M) = N X i=1 log "K X k=1 P(Mk)X all z PYi|z, MkPz|Mk# = N X i=1 log "K X k=1 wkX all z πk z1bk z1(yi11)· · · bk z1(yi1C) T Y t=2 hak zt−1ztbk zt(yit1)· · · bk zt(yitC)i#. (7) Compared to the usual hidden Markov model, there is an additional summation over the clusters in Equation 7, which seems to make the computations less straightforward than in the non-mixture case. Fortunately, redefining MHMM as a special type HMM allows us to use standard HMM algorithms without major modifications. We combine the Ksubmodels into one large hidden Markov model consisting of PK k=1 Skstates, where the initial state vector 8seqHMM: Mixture Hidden Markov Models for Sequence Data in R contains elements of the form wkπk. Now the transition matrix is block diagonal A=      A10· · · 0 0A2· · · 0 . . .. . ..... . . 0 0 · · · AK       ,(8) where the diagonal blocks Ak, k = 1, . . . , K, are square matrices containing the transition probabilities of one cluster. The off-diagonal blocks are zero matrices, so transitions between clusters are not allowed. Similarly, the emission matrices for each channel contain stacked emission matrices Bk. Covariates and cluster probabilities Covariates can be added to MHMMs to explain cluster memberships as in latent class analysis. The prior cluster probabilities now depend on the subject’s covariate values xiand are defined via the multinomial logit model: P(Mk|xi) = wik =ex> iγk 1 + PK j=2 ex> iγj .(9) The first submodel is set as the reference by fixing γ1= (0,...,0)>. As in MHMMs without covariates, we can still use standard HMM algorithms with a slight modification; we now allow initial state probabilities πto vary between subjects, i.e., for subject iwe have πi= (wi1π1, . . . , wiK πK)>. Of course, we also need to estimate the coefficients γ. For direct numerical maximization the modifications are straightforward. In the EM algorithm, regarding the M-step for γ,seqHMM uses the iterative Newton’s method with analytic gradients and Hessian which are straightforward to compute given all other model parameters. This Hessian can also be used for computing the conditional standard errors of coefficients. For unconditional standard errors, which take account of possible correlation between the estimates of γand other model parameters, the Hessian is computed using a finite difference approximation of the Jacobian of the analytic gradients. The posterior cluster probabilities P(Mk|Yi,xi)are obtained as P(Mk|Yi,xi) = P(Yi|Mk,xi)P(Mk|xi) P(Yi|xi) =P(Yi|Mk,xi)P(Mk|xi) PK j=1 P(Yi|Mj,xi)P(Mj|xi)=Li k Li, (10) where Liis the likelihood of the complete MHMM for subject i, and Li kis the likelihood of cluster kfor subject i. These are straightforwardly computed from forward probabilities. Posterior cluster probabilities are used, e.g., for computing classification tables. 2.4. Important special cases The hidden Markov model is not the only important special case of the mixture hidden Markov model. Here we cover some of the most important special cases that are included in the seqHMM package. Journal of Statistical Software 15 0.055 0.033 0.012 0.014 0.084 0.027 0.19 0.99 0.014 0 0 0 single/childless/with parents single/childless/left home divorced/childless/left home married/childless/left home married/children/left home married/childless/with parents States with prob. < 0.05 Figure 4: Illustrating a hidden Markov model as a directed graph. Pies represent five hidden states, with slices showing emission probabilities of combinations of observed states. States with emission probability less than 0.05 are combined into one slice. Edges show the transition probabilities. Initial probabilities of hidden states are given below the pies. same features for different groups. The user has a lot of control over the layout, e.g., dimensions of the grid, widths and heights of the cells, and positions of the legends. We also provide a function mc_to_sc_data for the easy conversion of multi-channel sequence data into a single-channel representation. Plotting combined data is often useful in addition to (or instead of) showing separate channels. 3.3. Visualizing hidden Markov models For the easy visualization of the model structure and parameters, we propose plotting HMMs as directed graphs. Such graphs are easily called with the plot method, with an object of class ‘hmm’ as an argument. Figure 4illustrates a five-state HMM. The code for producing the plot is shown in Section 4.4. Hidden states are presented with pie charts as vertices (or nodes), and transition probabilities are shown as edges (arrows, arcs). By default, the higher the transition probability, the thicker the stroke of the edge. Emitted observed states are shown as slices in the pies. For gaining a simpler view, observations with small emission probabilities (less than 0.05 by default) can be combined into one category. Initial state probabilities are given below or next to the respective vertices. In the case of multi-channel sequences, the data and the model are converted into a single-channel representation with the mc_to_sc function. A simple default plot is easy to call, but the user has a lot of control over the layout. Figure 5 illustrates another possible visualization of the same model. The code is shown in Section 4.4. The ssplot function (see Section 3.2) also accepts an object of class ‘hmm’. The user can easily choose to plot observations, most probable paths of hidden states, or both. The function 16 seqHMM: Mixture Hidden Markov Models for Sequence Data in R 0.05 0.03 0.01 0.01 0.08 0.03 0.2 1 0.01 0 0 0 divorced/childless/left home divorced/childless/with parents divorced/children/left home married/childless/left home married/childless/with parents married/children/left home single/childless/left home single/childless/with parents single/children/left home single/children/with parents Figure 5: Another version of the hidden Markov model of Figure 4with a different layout and modified labels, legends, and colors. All observed states are shown. automatically computes hidden paths if the user does not provide them. Figure 6shows observed sequences with the most probable paths of hidden states given the model. Sequences are sorted according to multidimensional scaling scores computed from hidden paths. The code for creating the plot is shown in Section 4.4. The plot method works for ‘mhmm’ objects as well. The user can choose between an interactive mode, where the model for each (chosen) cluster is plotted separately, and a combined plot with all models in one plot. The equivalent to the ssplot function for MHMMs is mssplot. It plots stacked sequence plots separately for each cluster. If the user asks to plot more than one cluster, the function is interactive by default. 4. Examples with life course data In this section we show examples of using the seqHMM package. We start by constructing and visualizing sequence data, then show how HMMs are built and fitted for single-channel and multi-channel data, then move on to clustering with MHMMs, and finally illustrate how to plot HMMs. Throughout the examples we use the same biofam data described in Section 2.1. We use both the original single-channel data and a three-channel modification named biofam3c, which is included in the seqHMM package. For more information on the conversion, see the documentation of the biofam3c data. 4.1. Sequence data Before getting to the estimation, it is good to get to know the data. We start by loading the original biofam data as well as the three-channel version of the same data, biofam3c. We convert the data into the ‘stslist’ form with the seqdef function. We set the starting age at 15 and set the order of the states with the alphabet argument (for plotting). Colors of the states can be modified and stored as an attribute in the ‘stslist’ object – this way the Journal of Statistical Software 17 15 17 19 21 23 25 27 29 Age MarriageParenthoodResidenceHidden states single married divorced childless children with parents left home State 1 State 2 State 3 State 4 State 5 Observed and hidden state sequences, n = 2000 Figure 6: Using the ssplot function for an ‘hmm’ object makes it easy to plot the observed sequences together with the most probable paths of hidden states given the model. user only needs to define them once. R> library("seqHMM") R> data("biofam", package = "TraMineR") R> biofam_seq <- seqdef(biofam[, 10:25], start = 15, labels = c("parent", + "left", "married", "left+marr", "child", "left+child", "left+marr+ch", + "divorced")) R> data("biofam3c", package = "seqHMM") R> marr_seq <- seqdef(biofam3c$married, start = 15, alphabet = c("single", + "married", "divorced")) R> child_seq <- seqdef(biofam3c$children, start = 15, + alphabet = c("childless", "children")) 18 seqHMM: Mixture Hidden Markov Models for Sequence Data in R R> left_seq <- seqdef(biofam3c$left, start = 15, + alphabet = c("with parents", "left home")) R> attr(marr_seq, "cpal") <- c("violetred2", "darkgoldenrod2", "darkmagenta") R> attr(child_seq, "cpal") <- c("darkseagreen1", "coral3") R> attr(left_seq, "cpal") <- c("lightblue", "red3") Here we show codes for creating Figures 1,2, and 3. Such plots give a good glimpse into multi-channel data. Figure 2: Plotting state distributions We start by showing how to call the simple default plot of Figure 2in Section 3.3. By default the function plots state distributions (type = "d"). Multi-channel data are given as a list where each component is an ‘stslist’ corresponding to one channel. If names are given, those will be used as labels in plotting. R> ssplot(list("Marriage" = marr_seq, "Parenthood" = child_seq, + "Residence" = left_seq)) Figure 1: Plotting sequences Figure 1with the whole sequences requires modifying more arguments. We call for sequence index plots (type = "I") and sort sequences according to the first channel (the original sequences), starting from the beginning. We give labels to yand x-axes and modify the positions of y-labels. We give a title to the plot but omit the number of subjects, which by default is printed. We set the proportion of the plot given to legends and the number of columns in each legend. R> seq_data <- list(biofam_seq[1:10, ], marr_seq[1:10, ], child_seq[1:10, ], + left_seq[1:10, ]) R> ssplot(seq_data, type = "I", sortv = "from.start", sort.channel = 1, + ylab = c("Original", "Marriage", "Parenthood", "Residence"), + xtlab = 15:30, xlab = "Age", ylab.pos = c(1, 1.5), title.n = FALSE, + title = "Ten first sequences", legend.prop = 0.63, + ncol.legend = c(3, 1, 1, 1)) Figure 3: Plotting sequence data in a grid For using the gridplot function, we first need to specify the ‘ssp’ objects of the separate plots. Here we start by defining the first plot for women with the ssp function. It stores the features of the plot, but does not draw anything. We want to sort sequences according to multidimensional scaling scores. These are computed from optimal matching dissimilarities for observed sequences. Any dissimilarity method available in TraMineR can be used instead of the default (see the documentation of the seqdef function for more information). We want to use the same legends for the both plots, so we remove legends from the ‘ssp’ objects. Since we are going to plot to two similar figures, one for women and one for men, we can pass the first ‘ssp’ object to the update function. This way we only need to define the changes and omit everything that is similar. Journal of Statistical Software 19 These two ‘ssp’ objects are then passed on to the gridplot function. Here we make a 2×2 grid, of which the bottom row is for the legends, but the function can also automatically determine the number of rows and columns and the positions of the legends. R> ssp_f <- ssp(list(marr_seq[biofam3c$covariates$sex == "woman", ], + child_seq[biofam3c$covariates$sex == "woman", ], + left_seq[biofam3c$covariates$sex == "woman", ]), + type = "I", sortv = "mds.obs", with.legend = FALSE, title = "Women", + ylab.pos = c(1, 2, 1), xtlab = 15:30, ylab = c("Married", "Children", + "Residence")) R> ssp_m <- update(ssp_f, title = "Men", + x = list(marr_seq[biofam3c$covariates$sex == "man", ], + child_seq[biofam3c$covariates$sex == "man", ], + left_seq[biofam3c$covariates$sex == "man", ])) R> gridplot(list(ssp_f, ssp_m), ncol = 2, nrow = 2, byrow = TRUE, + legend.pos = "bottom", legend.pos2 = "top", row.prop = c(0.65, 0.35)) For more examples on visualization, see a supplementary vignette in package seqHMM on CRAN (Helske 2018c). 4.2. Hidden Markov models We start by showing how to fit an HMM for single-channel biofam data. The model is initialized with the build_hmm function which creates an object of class ‘hmm’. The simplest way is to use automatic starting values by giving the number of hidden states. R> sc_initmod_random <- build_hmm(observations = biofam_seq, n_states = 5) It is, however, often advisable to set starting values for initial, transition, and emission probabilities manually. Here the hidden states are regarded as more general life stages, during which individuals are more likely to meet certain observable life events. We expect that the life stages are somehow related to age, so constructing starting values from the observed state frequencies by age group seems like an option worth a try (these are easily computed using the seqstatf function in TraMineR). We construct a model with four hidden states using age groups 15–18, 19–21, 22–24, 25–27 and 28–30. The fit_model function uses the probabilities given by the initial model as starting values when estimating the parameters. Only positive probabilities are estimated; zero values are fixed to zero. Thus, the amount of 0.1 is added to each value in case of zero-frequencies in some categories (at this point we do not want to fix any parameters to zero). Each row is divided by its sum, so that the row sums equal 1. R> sc_init <- c(0.9, 0.06, 0.02, 0.01, 0.01) R> sc_trans <- matrix(c(0.80, 0.10, 0.05, 0.03, 0.02, 0.02, 0.80, 0.10, + 0.05, 0.03, 0.02, 0.03, 0.80, 0.10, 0.05, 0.02, 0.03, 0.05, 0.80, 0.10, + 0.02, 0.03, 0.05, 0.05, 0.85), nrow = 5, ncol = 5, byrow = TRUE) R> sc_emiss <- matrix(NA, nrow = 5, ncol = 8) R> sc_emiss[1, ] <- seqstatf(biofam_seq[, 1:4])[, 2] + 0.1 R> sc_emiss[2, ] <- seqstatf(biofam_seq[, 5:7])[, 2] + 0.1 20 seqHMM: Mixture Hidden Markov Models for Sequence Data in R R> sc_emiss[3, ] <- seqstatf(biofam_seq[, 8:10])[, 2] + 0.1 R> sc_emiss[4, ] <- seqstatf(biofam_seq[, 11:13])[, 2] + 0.1 R> sc_emiss[5, ] <- seqstatf(biofam_seq[, 14:16])[, 2] + 0.1 R> sc_emiss <- sc_emiss / rowSums(sc_emiss) R> rownames(sc_trans) <- colnames(sc_trans) <- rownames(sc_emiss) <- + paste("State", 1:5) R> colnames(sc_emiss) <- attr(biofam_seq, "labels") R> sc_trans State 1 State 2 State 3 State 4 State 5 State 1 0.80 0.10 0.05 0.03 0.02 State 2 0.02 0.80 0.10 0.05 0.03 State 3 0.02 0.03 0.80 0.10 0.05 State 4 0.02 0.03 0.05 0.80 0.10 State 5 0.02 0.03 0.05 0.05 0.85 R> round(sc_emiss, 3) parent left married left+marr child left+child left+marr+ch State 1 0.928 0.063 0.002 0.002 0.001 0.001 0.002 State 2 0.701 0.218 0.018 0.028 0.001 0.004 0.029 State 3 0.417 0.290 0.050 0.114 0.001 0.006 0.117 State 4 0.204 0.231 0.080 0.201 0.002 0.009 0.256 State 5 0.101 0.157 0.097 0.196 0.002 0.013 0.400 divorced State 1 0.001 State 2 0.001 State 3 0.005 State 4 0.018 State 5 0.034 Now, the build_hmm function checks that the data and matrices are of the right form. R> sc_initmod <- build_hmm(observations = biofam_seq, + initial_probs = sc_init, transition_probs = sc_trans, + emission_probs = sc_emiss) We then use the fit_model function for parameter estimation. Here we estimate the model using the default options of the EM step. R> sc_fit <- fit_model(sc_initmod) The fitting function returns the estimated model, its log-likelihood, and information on the optimization steps. R> sc_fit$logLik [1] -16781.99 Journal of Statistical Software 21 R> sc_fit$model Initial probabilities : State 1 State 2 State 3 State 4 State 5 0.986 0.000 0.014 0.000 0.000 Transition probabilities : to from State 1 State 2 State 3 State 4 State 5 State 1 0.786 0.175 0.0391 0.00000 0.0000 State 2 0.000 0.786 0.0751 0.07568 0.0631 State 3 0.000 0.000 0.8898 0.08342 0.0267 State 4 0.000 0.000 0.0000 0.78738 0.2126 State 5 0.000 0.000 0.0000 0.00136 0.9986 Emission probabilities : symbol_names state_names 0 1 2 3 4 5 6 7 State 1 1 0 0.00000 0.000 0.00000 0.0000 0.000 0.0000 State 2 1 0 0.00000 0.000 0.00000 0.0000 0.000 0.0000 State 3 0 1 0.00000 0.000 0.00000 0.0000 0.000 0.0000 State 4 0 0 0.00195 0.992 0.00581 0.0000 0.000 0.0000 State 5 0 0 0.21508 0.000 0.00000 0.0246 0.713 0.0474 As a multi-channel example we fit a 5-state model for the 3-channel data. Emission probabilities are now given as a list of three emission matrices, one for each channel. The alphabet function from the TraMineR package can be used to check the order of the observed states – the same order is used in the build functions. Here we construct a left-to-right model where transitions to earlier states are not allowed, so the transition matrix is upper-triangular. This seems like a valid option from a life-course perspective. Also, in the previous single-channel model of the same data the transition matrix was estimated almost upper triangular. We also give names for channels – these are used when printing and plotting the model. We estimate model parameters using the local step with the default L-BFGS algorithm using parallel computation with 4 threads. R> mc_init <- c(0.9, 0.05, 0.02, 0.02, 0.01) R> mc_trans <- matrix(c(0.80, 0.10, 0.05, 0.03, 0.02, 0, 0.90, 0.05, 0.03, + 0.02, 0, 0, 0.90, 0.07, 0.03, 0, 0, 0, 0.90, 0.10, 0, 0, 0, 0, 1), + nrow = 5, ncol = 5, byrow = TRUE) R> mc_emiss_marr <- matrix(c(0.90, 0.05, 0.05, 0.90, 0.05, 0.05, 0.05, 0.90, + 0.05, 0.05, 0.90, 0.05, 0.30, 0.30, 0.40), nrow = 5, ncol = 3, + byrow = TRUE) R> mc_emiss_child <- matrix(c(0.9, 0.1, 0.9, 0.1, 0.1, 0.9, 0.1, 0.9, 0.5, + 0.5), nrow = 5, ncol = 2, byrow = TRUE) R> mc_emiss_left <- matrix(c(0.9, 0.1, 0.1, 0.9, 0.1, 0.9, 0.1, 0.9, 0.5, + 0.5), nrow = 5, ncol = 2, byrow = TRUE) R> mc_obs <- list(marr_seq, child_seq, left_seq) 22 seqHMM: Mixture Hidden Markov Models for Sequence Data in R R> mc_emiss <- list(mc_emiss_marr, mc_emiss_child, mc_emiss_left) R> mc_initmod <- build_hmm(observations = mc_obs, initial_probs = mc_init, + transition_probs = mc_trans, emission_probs = mc_emiss, + channel_names = c("Marriage", "Parenthood", "Residence")) R> mc_initmod Initial probabilities : State 1 State 2 State 3 State 4 State 5 0.90 0.05 0.02 0.02 0.01 Transition probabilities : to from State 1 State 2 State 3 State 4 State 5 State 1 0.8 0.1 0.05 0.03 0.02 State 2 0.0 0.9 0.05 0.03 0.02 State 3 0.0 0.0 0.90 0.07 0.03 State 4 0.0 0.0 0.00 0.90 0.10 State 5 0.0 0.0 0.00 0.00 1.00 Emission probabilities : Marriage : symbol_names state_names single married divorced State 1 0.90 0.05 0.05 State 2 0.90 0.05 0.05 State 3 0.05 0.90 0.05 State 4 0.05 0.90 0.05 State 5 0.30 0.30 0.40 Parenthood : symbol_names state_names childless children State 1 0.9 0.1 State 2 0.9 0.1 State 3 0.1 0.9 State 4 0.1 0.9 State 5 0.5 0.5 Residence : symbol_names state_names with parents left home State 1 0.9 0.1 State 2 0.1 0.9 State 3 0.1 0.9 State 4 0.1 0.9 State 5 0.5 0.5 Journal of Statistical Software 23 R> mc_fit <- fit_model(mc_initmod, em_step = FALSE, local_step = TRUE, + threads = 4) We store the model as a separate object for the ease of use and then compute the BIC. R> hmm_biofam <- mc_fit$model R> BIC(hmm_biofam) [1] 28842.7 4.3. Clustering and mixture hidden Markov models When fitting mixture hidden Markov models, the starting values are given as lists, with one component per cluster. For multi-channel data, emission probabilities are given as a list of lists. Here we fit a model for two clusters with 5 and 4 hidden states. For the cluster with five states we use the same starting values as for the multi-channel HMM described earlier. Covariates are defined with the usual formula and data arguments. Here we use sex and birth cohort to explain cluster memberships. We fit a model using 100 random restarts of the EM algorithm followed by the local L-BFGS method. Again we use parallel computation. R> mc_init2 <- c(0.9, 0.05, 0.03, 0.02) R> mc_trans2 <- matrix(c(0.85, 0.05, 0.05, 0.05, 0, 0.90, 0.05, 0.05, 0, 0, + 0.95, 0.05, 0, 0, 0, 1), nrow = 4, ncol = 4, byrow = TRUE) R> mc_emiss_marr2 <- matrix(c(0.90, 0.05, 0.05, 0.90, 0.05, 0.05, 0.05, + 0.85, 0.10, 0.05, 0.80, 0.15), nrow = 4, ncol = 3, byrow = TRUE) R> mc_emiss_child2 <- matrix(c(0.9, 0.1, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5), + nrow = 4, ncol = 2, byrow = TRUE) R> mc_emiss_left2 <- matrix(c(0.9, 0.1, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5), + nrow = 4, ncol = 2, byrow = TRUE) R> mhmm_init <- list(mc_init, mc_init2) R> mhmm_trans <- list(mc_trans, mc_trans2) R> mhmm_emiss <- list(list(mc_emiss_marr, mc_emiss_child, mc_emiss_left), + list(mc_emiss_marr2, mc_emiss_child2, mc_emiss_left2)) R> biofam3c$covariates$cohort <- cut(biofam3c$covariates$birthyr, + c(1908, 1935, 1945, 1957)) R> biofam3c$covariates$cohort <- factor(biofam3c$covariates$cohort, + labels = c("1909-1935", "1936-1945", "1946-1957")) R> init_mhmm <- build_mhmm(observations = mc_obs, initial_probs = mhmm_init, + transition_probs = mhmm_trans, emission_probs = mhmm_emiss, + formula = ~ sex + cohort, data = biofam3c$covariates, + channel_names = c("Marriage", "Parenthood", "Residence"), + cluster_names = c("Cluster 1", "Cluster 2")) R> set.seed(1011) R> mhmm_fit <- fit_model(init_mhmm, local_step = TRUE, threads = 4, + control_em = list(restart = list(times = 100))) R> mhmm <- mhmm_fit$model 24 seqHMM: Mixture Hidden Markov Models for Sequence Data in R The summary method automatically computes some features for an MHMM, e.g., standard errors for covariates and prior and posterior cluster probabilities for subjects. A print method shows some summaries of these: estimates and standard errors for covariates (see Section 2.3), log-likelihood and BIC, and information on most probable clusters and prior probabilities. Parameter estimates for transitions, emissions, and initial probabilities are omitted by default. The classification table shows mean probabilities of belonging to each cluster by the most probable cluster (defined from the posterior cluster probabilities). A good model should have values close to 1 on the diagonal. R> summary(mhmm, conditional_se = FALSE) Covariate effects : Cluster 1 is the reference. Cluster 2 : Estimate Std. error (Intercept) -1.209 0.138 sexwoman 0.213 0.141 cohort1936-1945 -0.785 0.172 cohort1946-1957 -1.238 0.165 Log-likelihood: -12969.57 BIC: 26592.66 Means of prior cluster probabilities : Cluster 1 Cluster 2 0.857 0.143 Most probable clusters : Cluster 1 Cluster 2 count 1753 247 proportion 0.876 0.124 Classification table : Mean cluster probabilities (in columns) by the most probable cluster (rows) Cluster 1 Cluster 2 Cluster 1 0.9775 0.0225 Cluster 2 0.0013 0.9987 4.4. Visualizing hidden Markov models The figures in Section 3.3 illustrate the five-state multi-channel HMM fitted in Section 4.2. A basic HMM graph is easily called with the plot method. Figure 7illustrates the default plot. R> plot(hmm_biofam) Journal of Statistical Software 31 Müller NS, Studer M, Ritschard G (2007). “Classification de Parcours de Vie à l’Aide de l’Optimal Matching.” XIVe Rencontre de la Société Francophone de Classification (SFC 2007), pp. 157–160. Nocedal J (1980). “Updating Quasi-Newton Matrices with Limited Storage.” Mathematics of Computation,35(151), 773–782. doi:10.1090/s0025-5718-1980-0572855-7. O’Connell J, Højsgaard S (2011). “Hidden Semi Markov Models for Multiple Observation Sequences: The mhsmm Package for R.” Journal of Statistical Software,39(4), 1–22. doi: 10.18637/jss.v039.i04. Rabiner L (1989). “A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition.” Proceedings of the IEEE,77(2), 257–286. doi:10.1109/5.18626. RCore Team (2018). R: A Language and Environment for Statistical Computing.RFoundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/. Rinnooy Kan AHG, Timmer GT (1987a). “Stochastic Global Optimization Methods Part I: Clustering Methods.” Mathematical Programming,39(1), 27–56. doi:10.1007/ bf02592070. Rinnooy Kan AHG, Timmer GT (1987b). “Stochastic Global Optimization Methods Part II: Multi-Level Methods.” Mathematical Programming,39(1), 57–78. doi:10.1007/ bf02592071. Turner R (2018). hmm.discnp: Hidden Markov Models with Discrete Non-Parametric Observation Distributions.Rpackage version 2.1-3, URL https://CRAN.R-project.org/ package=hmm.discnp. van de Pol F, Langeheine R (1990). “Mixed Markov Latent Class Models.” Sociological Methodology,20, 213–247. doi:10.2307/271087. Vermunt JK, Tran B, Magidson J (2008). “Latent Class Models in Longitudinal Research.” In S Menard (ed.), Handbook of Longitudinal Research: Design, Measurement, and Analysis, pp. 373–385. Elsevier, Burlington, MA. Visser I, Speekenbrink M (2010). “depmixS4: An RPackage for Hidden Markov Models.” Journal of Statistical Software,36(7), 1–21. doi:10.18637/jss.v036.i07. Ypma J (2018). nloptr:RInterface to NLopt.Rpackage version 1.2.1, URL https://CRAN. R-project.org/package=nloptr. 32 seqHMM: Mixture Hidden Markov Models for Sequence Data in R A. Notations Symbol Meaning YiObservation sequences for subject i, i = 1, . . . , N yit Observations for subject iat time t, t = 1, . . . , T yitc Observation for subject iat time tin channel c, c= 1, . . . , C mc∈ {1, . . . , Mc}Observed state space for channel c zit Hidden state at time tfor subject i s∈ {1, . . . , S}Hidden state space A={asr}Transition matrix of size S×S asr =P(zt=r|zt−1=s)Transition probability between hidden states sand r Bc={bs(mc)}Emission matrix of size S×Mcfor channel c bs(mc) = P(yitc =mc|zit =s)Emission probability for observed state mcin channel c given hidden state s bs(yit) = bs(yit1)· · · bs(yitC)Joint emission probability of observations at time tin channels 1, . . . , C given hidden state s π= (π1, . . . , πS)Vector of initial probabilities πs=P(z1=s)Initial probability for hidden state s ˆzi(Yi)The most probable hidden state sequence for subject i xiCovariates for subject i Mk, k = 1, . . . , K Submodel for cluster k(latent class/cluster) wik Probability of cluster kfor subject i γkRegression coefficients for cluster k {πk, Ak, Bk 1, . . . , Bk C, γk}Model parameters for cluster k Table 2: Notations. Affiliation: Satu Helske Institute for Analytical Sociology Linköping University 601 74 Norrköping, Sweden E-mail: [email protected] Jouni Helske Department of Science and Technology Linköping University 601 74 Norrköping, Sweden E-mail: [email protected] Journal of Statistical Software http://www.jstatsoft.org/ published by the Foundation for Open Access Statistics http://www.foastat.org/ January 2019, Volume 88, Issue 3 Submitted: 2016-01-14 doi:10.18637/jss.v088.i03 Accepted: 2017-11-22