Analyzing and Visualizing State Sequences in R with TraMineR

This article describes the many capabilities oﬀered by the TraMineR toolbox for categorical sequence data. It focuses more speciﬁcally on the analysis and rendering of state sequences. Addressed features include the description of sets of sequences by means of transversal aggregated views, the computation of longitudinal characteristics of individual sequences and the measure of pairwise dissimilarities. Special emphasis is put on the multiple ways of visualizing sequences. The core element of the package is the state sequence object in which we store the set of sequences together with attributes such as the alphabet, state labels and the color palette. The functions can then easily retrieve this information to ensure presentation homogeneity across all printed and graphical displays. The article also demonstrates how TraMineR ’s outcomes give access to advanced analyses such as clustering and statistical modeling of sequence data.


Introduction
This article is concerned with categorical sequence data and more specifically with state sequences, where the position of each successive state receives a meaningful interpretation in terms of age, date, or more generally of elapsed time or distance from the beginning of the sequence.Its aim is to examine a series of questions about such state sequences and to present the various solutions that we implemented in the R (R Development Core Team 2011) package TraMineR for answering them.
The addressed methods are for sets of sequences and most of them are holistic (Billari 2001b) in that they consider each sequence as a whole; i.e., as a conceptual unit.The discussion is mainly oriented towards the analysis of sequences describing individual life courses.Nevertheless, most of the discussed concepts and tools should be applicable in other domains such as text, biology, quality control or web logs analysis, to cite just a few.
Sequences are complex objects, and we need special tools for describing and displaying them.We consider, therefore, questions regarding the exploration and description of sets of sequences such as: • Which characteristics of sequences are we interested in?
• What kind of indicators can we compute for a sequence set?
• What are suited plots for rendering sequences?
• How can we measure similarity between sequences?With a more analytical or explanatory concern, we also consider issues such as: • How can we identify groups with similar patterns and build typologies of sequences?
• How can we analyze the relationship of sequences with covariates?
In the social sciences, state sequences are of interest for studying life trajectories such as occupational histories, professional careers or cohabitational life courses.Some of the typical questions arising in this area are: • Do life courses obey some social norm?Which are the standard trajectories?What kind of departures do we observe from these standards ?How do life course patterns evolve over time ?
• Why are some people more at risk to follow a chaotic trajectory or to stay stuck in a state?How does the trajectory complexity evolve across birth cohorts?
• How is the life trajectory related to sex, social origin and other cultural factors?
Empirical answers to such questions require us to consider collections of life sequences, to examine them from both a transversal and a longitudinal perspective, and to study their relationships with covariates.
The primary objective of sequence methods is then to extract simplified workable information from sequential data sets; that is, to efficiently summarize and render these sets and to categorize the sequential patterns into a limited number of groups.This is essentially an exploratory task that consists of computing summary indicators, as well as sorting, grouping and comparing sequences.The resulting groups and real-value indicators may then be submitted to classical inferential methods and serve, for instance, as response variables or explanatory factors for regression-like models.
A common approach for categorizing patterns consists of computing pairwise distances between them by means of sequence alignment algorithms (such as optimal matching) or other suitable metrics and using this information for clustering the sequences.This method has been applied to various data since the pioneering work of Abbott and Forrest (1986).A review can be found in Abbott and Tsay (2000).The expected outcome of such a strategy is a typology, with each cluster grouping cases with similar trajectories.Through binary logistic regression or classification trees, for example, we can then study how each cluster membership is related to covariates.
A more recent complementary approach considered in the literature (Elzinga and Liefbroer 2007;Widmer and Ritschard 2009) is to focus on sequence indicators measuring for instance the longitudinal diversity and complexity of the sequences and to analyze them by means of conventional statistical tools for real-value variables.
With a somewhat more aggregated point of view, an approach considered, for instance, by Billari (2001a) consists of looking at the sequence of transversal characteristics measured at each position, such as the diversity of states observed at each given age.Comparing the evolution of such transversal characteristics for different groups defined by birth cohorts or sex, for instance, provides instructive insights (Widmer and Ritschard 2009).However, when working with transversal indicators we lose the specific information on individual follow-ups.
We may indeed imagine many other ways of looking at categorical sequences such as correspondence analysis of the states (Deville and Saporta 1983) or advanced Markov modeling (Berchtold and Raftery 2002); i.e., the study of how the probability of a given state depends on the previously observed states.Transforming state sequences into event sequences and resorting to tools for mining frequent subsequences permits us to gain interesting knowledge about the typical sequencing of states or events (Billari, Frnkranz, and Prskawetz 2006;Ritschard, Gabadinho, Müller, and Studer 2008).A very common approach in the life course literature is event history or survival analysis (Mayer and Tuma 1990;Yamaguchi 1991;Hosmer and Lemeshow 1999;Blossfeld, Golsch, and Rohwer 2007) which focuses on the occurrence of a specific event or somewhat equivalently on the duration-time to event-until a given state transition.Though not addressed here, all these techniques may usefully complement the considered state sequence techniques.
The TraMineR R package is available from the Comprehensive R Archive Network at http: //CRAN.R-project.org/package=TraMineR and offers many analysis and visualization tools for either state or event sequences.These tools include already known methods, as well as new developments.We focus in this article on the functions intended for state sequence analysis.The paper is organized as follows.In Section 2, we introduce the TraMineR library, describe the mvad data set used for illustration and give a first example of analysis that can be run with TraMineR.Section 3 defines the different forms of sequential data that are supported by the package.In Section 4, we introduce the central concept of state sequence object.Section 5 introduces two basic visualization tools and describes the general plotting principles used by the package.Section 6 is devoted to the summarization and visual rendering of sets of sequences, while Section 7 is concerned with individual sequence indicators.In Section 8, we present the metrics that were implemented for measuring pairwise dissimilarities between sequences.In Section 9, we illustrate how dissimilarities measures can be used for further statistical analysis.Finally, we make some concluding remarks in Section 10.

The TraMineR R package
TraMineR (Gabadinho, Ritschard, Studer, and Müller 2009) is a package for mining and visualizing sequences of categorical data describing life courses in R (R Development Core Team 2011), the name TraMineR being a contraction of Life Trajectory Miner for R. It puts together most of the features proposed separately by other software for sequential data and offers many original tools for managing, analyzing and rendering categorical sequences.
They all compute the optimal-matching edit distance between pairs of sequences and each of them offers specific useful facilities for describing sets of sequences.TraMineR is, to our knowledge, the first such toolbox for the free R statistical and graphical environment.Its salient characteristics are: • R and TraMineR are free and open source; • Since TraMineR is developed in R, it takes advantage of many already optimized procedures of R as well as of its powerful graphical capabilities; • R runs under several OS including Linux, MacOS X, Unix and Windows; any R script with TraMineR functions runs unmodified under all operating systems; • Specific TraMineR functions can be combined in the same script with any of the numerous basic statistical procedures of R as well as with those of any other R-package.
TraMineR is readily installed from within R via install.packages("TraMineR").It features a unique set of procedures for analyzing and visualizing state sequence data, such as: • Handling a large number of state sequence representations with simple functions for transforming to and from different formats; • A whole series of easy to use plot functions for rendering sets of sequences (density plot, frequency plot, index plot, representative sequence plot and more); • Individual longitudinal characteristics of sequences (length, time in each state, longitudinal entropy, complexity, turbulence and more); • Sequence of transversal characteristics by position (transversal state distribution, transversal entropy, modal state); • Other aggregated characteristics (transition rates, average duration in each state, sequence frequency); • A choice of metrics for evaluating distances between sequences.
Table 1 gives an overview of the key functions for analyzing state sequences that will be described in the remainder of the article.It is worth mentioning that the package provides also tools for event sequences such as finding the most frequent and most discriminating subsequences and extracting association rules between subsequences, as well as tools for ANOVAlike analyses of sequences (Studer, Ritschard, Gabadinho, and Müller 2011) that will not be discussed here.See the User's guide (Gabadinho et al. 2009) for a detailed description of the package usage.
To continue, we examine how the diversity of states within each sequence is related to sex, to whether the father is unemployed and to whether the qualification grade at end of compulsory school was good.We compute the longitudinal entropy and regress it on the covariates: Results show that males experience less diverse states and that youngsters with good grades at the end of compulsory school experience more diverse states.Whether the father is unemployed does not have a significant effect.

Sequence representations
State sequences can be represented in many different ways, depending on the data source and on how the information is organized.Data organization and conversion between formats is discussed in detail in Ritschard, Gabadinho, Studer, and Müller (2009), where an ontology of longitudinal data presentations is given that may help identify the kind of data at hand.Here, we limit the discussion to the sequence data representations that TraMineR can handle and import.Those formats are listed in Table 4 together with the conversions that can currently be done with the provided seqformat() function.

State sequences
We consider sequences of discrete or categorical data.Formally, we define a state sequence of length ℓ as an ordered list of ℓ elements successively chosen from a finite set A of size a = |A| that is called the alphabet.A natural way of representing a sequence x is by listing the successive elements that form the sequence x = (x 1 , x 2 , . . ., x ℓ ), with x j ∈ A. With reference to this expanded form of sequences, state sequences are characterized by two properties.Firstly, they are formed by elements that are states; i.e., something that can last as opposed, for instance, to events that occur at given time points.Secondly, the position of each element conveys meaningful information in terms of age, date or, more generally, elapsed time or distance from the beginning of the sequence.Position indexes providing time information may be either absolute calendar values (day, year, month, ...) or relative process time (age, process duration, ...).
In TraMineR, the expanded form is called STate-Sequence (STS) format.In this format, the successive states (statuses) of an individual are given either in consecutive columns, or as a character string with states separated by a given symbol such as '-' or '/', the former being the default separator.Each position (column) is supposed to correspond to a predetermined time unit.

Other sequence representations
Sequence data can be represented in more compact ways than STS essentially by giving only one of several same successive states.In that case, we have to explicitly stamp the successive distinct states with their starting position or duration.  of two sequences in the different formats.The two considered sequences describe family formation histories of two individuals, the states being single (S), married (M), married with children (MC) and divorced (D).
A first efficient way of representing a state sequence is by listing the distinct successive states with their associated durations.We get thus a sequence of couples (x j , t j ) where x j is a state and t j its duration; i.e., the number of times it is repeated.This is the State-Permanence-Sequence (SPS) format (Aassve, Billari, and Piccarreta 2007).By considering only the distinct successive states without their associated durations, we get the Distinct-Successive-States (DSS) sequence format.This DSS form holds the basic state sequencing information, but loses all time (t j ) and, more generally, alignment data.
In the SPELL format there is one line for each spell.Each spell is characterized by the state (supposed constant during the spell) and the spell start and end times.STS sequences can easily be derived from such representation.

State sequence objects
The general philosophy of the library is to ensure that the various results and plots outputted for a same set of sequence data use the same state labels and colors.Likewise, any information on case weights and possible missing information about some positions in sequences should be treated the same way throughout the analysis.To achieve this goal, the TraMineR functions for state sequence analysis require as an argument a state sequence object that includes both the sequential data and its attributes.Thus, the first step when using TraMineR for state sequence analysis is to create a state sequence object.This is done with the seqdef() function from data organized in either of the STS, SPS or SPELL forms described in the previous section.set introduced in Section 2. The main attributes are listed in Table 5, together with their default values and the dedicated functions to retrieve or set them.

Creating state sequence objects
In the mvad data set, the retained activity status variables are stored in columns 17 to 86.We display these statuses for the first six considered months (September 1993 to February 1994) of the first two records: R> The default input format for the seqdef() function is STS, which is appropriate for the mvad data set.If the input data is in another format it must be specified with the informat argument and seqdef() will automatically make the required conversion.

Alphabet and state labels
The alphabet is the list of states allowed in the sequences.Both short and long labels of the states forming the alphabet are attached to the object.Long labels serve mainly for color legends in plots, while short state names are primarily used in printed outputs.Shorter names produce cleaner and shorter output when printing the sequences.
By default, the sorted list of the distinct states found in the input data (as returned by the seqstatl() function) defines the alphabet and is used as state names and labels.This can be changed with optional arguments, which is necessary, for example, when the alphabet contains states that do not appear in the retained sequences.
Below we specify short state names with the states argument and long state labels with labels.These arguments expect vectors of names or labels that are ordered conformably with the alphabet.The alphabet argument can be used to change the order of the states, in which case the vectors passed with the states and labels arguments should conform to the newly defined order.

Other important attributes and properties
We briefly comment here upon some other important attributes that will be used in conjunction with the alphabet and state labels by TraMineR's functions.

State colors and position names
The sequence plot functions provided by the library need a distinct color for each state.A color palette is therefore attached to the sequence object.A default color palette from RColorBrewer (Neuwirth 2007) is automatically selected as long as the alphabet size does not exceed 12. Position names serving mainly for labeling the ticks of the x-axis but that are also useful for increasing readability of tabulated output are also an attribute of the object.If left unspecified, position names are taken from the corresponding column names of the original data frame.The interval between the x-axis tick-marks is an additional attribute that can be set (xtstep argument) for optimizing rendering.

Case weights
Survey data often come with case weights that account for the sampling scheme and unit nonresponses.Using such case weights is important to compensate for sampling bias and thus get results that are more realistic.When weights are attached to the state sequence object, the TraMineR functions that can handle weights automatically produce weighted results.To disable the use of weights, add option weighted=FALSE to the function.
The weight variable in mvad contains case weights that account for the selective attrition during the survey and we attach them to the sequence object as shown below.Unless otherwise specified, we will use this weighted sequence object from here on.

Missing values
Missing values in the expanded (STS) form of a sequence occur, for example, when: • Sequences do not start on the same date while using a calendar time axis; • The follow-up time is shorter for some individuals than for others yielding sequences that do not end up at the same position; • The observation at some positions is missing due to nonresponse, yielding internal gaps in the sequences.
The way missing values should be handled may be different for each of these situations.In the first case, we may want to maintain explicitly the starting missing values to preserve alignment across sequences or possibly left-align sequences by switching to a process time axis.In the second case, the ending missing terms could just be ignored.
To allow such differentiated treatments, TraMineR distinguishes left, in-between and right missing values.We can specify how each of the missing types should be encoded with the left, gaps and right arguments.By default, gaps and left-missing states are coded as explicit missing values while all missing values encountered after the last valid (rightmost) state in a sequence are considered void elements; i.e., the sequence is considered to end after the last valid state.
The specific treatment of each type of missing value will depend upon whether the analysis method envisaged supports missing values; and, if yes, which kind it supports.Most of the proposed functions, such as seqdist() for computing distances between sequences, have optional arguments for dealing with missing states.

Subsets and attributes inheritance
Subsets of sequence objects can be defined by specifying row and column indexes (or names) as for R matrices and data frames.Every subset of a state sequence object inherits its 'parent' attributes.The alphabet and color palette, for instance, remain the same for all subsets.This is of particular importance when comparing graphics that render different subsets of a same sequence object.

Visualizing individual state sequences
State sequence visualization is one of the most important features of the package.This section introduces two basic plotting functions, namely the index plot intended to render a set or subset of individual state sequences and the frequency plot that visualizes them according to their frequencies.We explain also the common design of most of TraMineR's plotting functions.

Sequence index plots
A sequence index plot (Figure 2) individually renders the selected state sequences.Each of them is represented by horizontally stacked boxes that are colored according to the state at the successive positions.The resulting bars are put above each other to vertically align the  positions.We thus visualize, for each case, the individual longitudinal succession of states as well as, through the length of each color segment, the duration spent in each successive state.
The alignment also permits easy transversal comparisons at each position.The sequence index plot shown in Figure 2 was obtained with the command below: 1 R> seqiplot(mvad.seq,border = NA, with.legend= "right") Since we have attached case weights to the mvad.seqsequence object, the width of the bar representing each sequence is proportional to its weight.This default behavior could be changed with the weighted=FALSE argument.The plotted sequences are selected with the idxs argument by providing either a vector of indexes, or 0 for requesting all the sequences.The default value is 1:10 and Figure 2 displays therefore only the first 10 sequences of mvad.seq.The seqIplot() alias produces full index plots that display all the sequences in the set without spaces between sequences and without borders around unit states.The usefulness of such plots has, for instance, been stressed by Scherer (2001) and Brzinsky-Fay et al. (2006).However, when the number of displayed sequences is large, they may produce burden pictures that are often hard to interpret. 2We can partially overcome this drawback by sorting the sequences according to the values of a suitably chosen covariate-passed with the sortv argument.Good choices are, for instance, the distance to the most frequent sequence or the scores of a multidimensional scaling analysis 3 of the dissimilarities between sequences 1 seqiplot(), as most other plotting functions described in this paper, is just an alias for calling a generic seqplot() state sequence plot function with the appropriate type argument and suitable default option values.The border=NA option suppresses the border that surrounds, by default, each unit state in the sequence.
2 When plotting several hundred of sequences, saving index plots may also produce heavy files in vectorial formats such as PostScript and PDF; generating plots in bitmap formats such as PNG or JPEG is recommended in such cases.
3 The scores are obtained from the dissimilarity matrix with the cmdscale() function.(Figure 3).Both solutions suppose that we can compute such dissimilarities; this will be addressed in Section 8.

Sequence frequencies
The seqtab() function returns a table with the counts and percent frequencies of the sequences sorted in decreasing order of their frequencies.In the next example, we request the four most frequent sequences of mvad.seq with idxs=1:4.In the printed outcome, sequences are displayed in the shorter and more readable SPS format: R> seqtab(mvad.seq,idxs = 1:4) The most frequent sequence in the mvad.seqobject is a spell of two years of school followed by 46 months of higher education.It accounts, however, for only 4.7% of the total weights of the 712 cases considered.The second most frequent sequence, which concerns 3.5% of the weighted individuals, is indeed very similar to the previous one.

Sequence frequency plots
A graphical view of the sequence frequency table where bar widths are proportional to the frequencies is obtained with the seqfplot() function.Figure 4 shows the plot of the weighted and unweighted frequencies4 obtained with:  If we look at the unweighted results, the most frequent sequence is to stay employed during the entire follow-up period (be in state EM during 70 months).This sequence, which was the fourth most frequent in the weighted frequency table with 2.5% of the total weight, accounts for 7% of the 712 cases considered.
The probability for two individuals to follow exactly the same 70-month trajectory is small, yielding a large number of different patterns.The 10 most frequent sequences account for only about 20% of all the trajectories, which reflects this high diversity.

Reading and controlling state sequence plots
The way index plots render individual state sequences with horizontally stacked boxes is common to other functions of the library that visualize specific state sequences.The position in the sequence is read on the x-axis.The first value on this axis is the selected origin.The sequence is read from left to right in the same way as printed outputs.Tick labels for the x-axis are retrieved, by default, from the plotted sequence object.
The values on the y-axis are the indexes of the plotted sequences.The index refers to the considered ranking of the sequences.For instance, in sequence index plots, the default order is that in the state sequence object unless a specific sort variable is provided with the sortv argument.In sequence frequency plots, sequences are sorted according to their frequency in the data set, while in representative sequence plots (Section 9.1), sequences are sorted according to their representativeness score.

Analyzing and Visualizing State Sequences in R with TraMineR
The indexes on the y-axis-and hence the sequences-are displayed bottom-up.Thus, when sequences are sorted, the first ranked one is at the bottom of the plot.5This respects the usual standard for y-axes.It may, however, be confusing when compared with the corresponding printed outputs where sequences are displayed top-down.
Other aspects of the graphic (title, font size, axes display, axis label, state legend, ...) can be controlled with dedicated options described in detail in the reference manual.There is also an option to produce separate plots by levels of a covariate.

Computing and plotting overall and transversal statistics
We now turn to the facilities offered by TraMineR for visualizing and computing overall and transversal descriptive statistics of a set of sequences.The functions discussed here all require a state sequence object as main argument and admit a series of optional parameters.
We illustrate with the mvad.seqweighted sequence object created on page 11.

Overall statistical characteristics
We consider, first, global synthesized information that is based neither on individual longitudinal characteristics nor on transversal characteristics by position.More specifically, we focus on the overall state distribution and transition rates between states.

Mean time spent in each state
A first synthetic information is given by the mean-not necessarily consecutive-time spent in the different states; that is, the mean number of times each state is observed in a sequence.This characterizes the overall state distribution.As an example, we plot the mean times for two subsets defined by the funemp covariate that indicates whether the respondent's father was unemployed at the time of the survey (Figure 5).The graphic with the distinct plots by levels of the funemp covariate is obtained by passing funemp as group argument to the plotting function.This option is common to all the plotting functions presented in this article.
R> seqmtplot(mvad.seq,group = mvad$funemp, ylim = c(0, 30)) We can see that the mean time spent in joblessness and training is higher for interviewees with unemployed fathers, while the time they spent in 'school', 'further education' and 'higher education', is lower.
Mean time values are obtained with the seqmeant() function.

Transition rates
Another interesting information about a set of sequences is the transition rate between each couple of states (s i , s j ); i.e., the probability to switch at a given position from state s i to state s j .Let n t (s i ) be the number of sequences that do not end in t with state s i at position t and let n t,t+1 (s i , s j ) be the number of sequences with state s i at position t and state s j at position t + 1.The transition rate p(s j | s i ) between states s i and s j is obtained as .
with L the maximal observed sequence length.
The seqtrate() function returns the matrix of transition rates for the provided sequence object.By default, the rates are assumed position-independent; i.e., the same whatever t.
The outcome is a single matrix where each row i gives a transition distribution from the originating state s i in t to the states in t + 1; that is, each row total equals one.Hence, transition rates provide information about the most frequent state changes observed in the data together with, on the diagonal, an assessment of the stability of each state.
In the following example we compute the transition rate matrix for the mvad.seqsequence object: R> mvad.trate<-seqtrate(mvad.seq)R> round(mvad.trate, 2) Time-varying transition rates can be obtained with option time.varying=TRUE, in which case a 3-dimensional array with a distinct transition rate matrix for each of the positions t = 1, 2, . . ., L − 1 is returned.The matrix for position t is computed by considering only the states at t and t + 1.The third dimension of the array corresponds to the position t index.

Transversal state distributions
Time varying transition rates are transversal characteristics computed at the successive considered positions.In the same vein, it is of interest to look at the transversal distribution of the states at each position of the considered sequences.A state distribution plot, as produced by seqdplot(), displays the general pattern of the whole set of trajectories.When interpreting such graphics, one must remember that, unlike sequence index plots and sequence frequency plots, they do not render individual sequences or individual follow-ups.They provide aggregated views made of successive slices, each of which represents transversal characteristics.

Sequence of modal states
An interesting summary that can be derived from the state distributions is the sequence made of the most frequent state at each position.It is obtained with the seqmodst() function and plotted with seqmsplot().Figure 7 shows how such modal state sequences are displayed.The height of the bar at each position is proportional to the frequency of the displayed state at that position.The number of occurrences of the modal state sequence is also displayed.Since the shown sequences of modal states do not belong to the sequence dataset, the number of occurrences is 0 for both considered groups.

Transversal entropy of state distributions
In addition to the state distribution, the seqstatd() function provides for each position in the sequence the number of valid states and the Shannon entropy of the transversal state distribution.Shannon's entropy, also known as the entropy index, has been applied to social science data by, for instance, Billari (2001a) and Fussell (2005).Letting p i denote the proportion of cases in state i at the considered position, the entropy is where a is the size of the alphabet.The entropy is 0 when all cases are in the same state and is maximal when we have the same proportion of cases in each state.The entropy can be seen as a measure of the diversity of states observed at the considered position.
Plotting the transversal entropies can be useful to find out how the diversity of states evolves along the time axis.We plot transversal entropies with seqHtplot().Figure 8 shows the curves by end of compulsory school qualification group.For the first group, the entropy of the state distributions noticeably decreases at the end of the follow-up period.This is a consequence of the increasing proportion of youngsters entering into full employment (Figure 6).For the second group, the entropy index slightly increases at the end of the considered period, which may be explained by the emergence of two balanced subgroups, namely those who continue higher education and those who enter into employment.We focus now on the characterization and summarization of longitudinal characteristics of individual sequences.Essentially, the aim is to define measures that inform on how each sequence is constituted; i.e., on whether it takes a simple or more complex form.

Individual sequence characteristics
The interpretation of complexity indexes will depend indeed on the context.Consider for instance the number of transitions-changes of state-in a sequence.When looking at work trajectories, for example, sequences with numerous transitions may correspond to unusual disrupted trajectories.In other contexts such as family formation, sequences with fewer transitions may indicate that an individual failed to pass through the usual stages of the family formation (leaving the parental home, cohabitation with a partner, birth of one or more children, etc...).
In the SPS form (see Section 3) a state sequence is represented as an ordered list of successive distinct states with their associated durations; i.e., as a sequence of couples (x j , t j ) where x j is a state and t j its duration.This suggests that we can distinguish characteristics of the state sequencing-the distinct successive states (DSS)-from those of the durations. 8We first examine two indicators of the state sequencing and one based on the durations.More synthetic measures are addressed in Section 7.2.

Number of transitions
Perhaps the simplest indicator is the number of transitions in the sequence; i.e., the number of state changes.The number of transitions in a sequence x is readily obtained from the length ℓ d (x) of its DSS sequence.It is ℓ d (x) − 1.We get the number of transitions for each sequence of state sequence object with seqtransn().

Number of subsequences
8 The two pieces of information can be extracted separately with seqdss() and seqdur().The number ϕ(x) of subsequences that can be extracted from the DSS sequence provides also useful information on the sequencing of the states.This measure is returned by the seqsubsn() function.It is used in the turbulence measure presented below.A subsequence y of x is composed of elements of x occurring in the same order than in x.
The maximal number of subsequences is reached only for a sequence made of the repetition of the alphabet.In Figure 9, for example, sequences 5 and 9 have the maximal number of transitions, while the number of subsequences is maximal for sequence 9 only.

Within sequence entropy
Regarding the durations, we consider the total time spent in each state; i.e., in case of multiple spells in a same state, the sum of the lengths of these spells.For example, in (EM,4)-(TR,2)-(EM,64), the first sequence in the object mvad.seq,there are two spells in state EM with respective durations 4 and 64.Hence, the time spent in state EM is 68 months, as shown by the output of the seqistatd() function:

R> seqistatd(mvad.seq[1:4, ])
EM FE HE JL SC TR 1 68 0 0 0 0 2 2 0 36 34 0 0 0 3 10 34 0 2 0 24 4 14 0 0 9 0 47 The total time spent in each state characterizes the state distribution within a sequence.The entropy of this distribution can be seen as a measure of the diversity of its states.We call it within or longitudinal entropy to distinguish from the transversal entropy considered in Section 6.2 on page 20.
The seqient() function returns the longitudinal Shannon entropies; i.e., for each sequence the value of where a is the size of the alphabet and Ã i the proportion of occurrences of the ith state in the considered sequence.When the state remains the same during the whole sequence, the entropy equals 0, while the maximum entropy is reached when the same time is spent inside the sequence in each possible element of the alphabet.By default the entropy is normalized by dividing the value of h(Ã 1 , . . ., Ã s ) by its theoretical maximum, log a.9 Figure 9 helps to get a more concrete idea of what the entropy measures.We see that the within-sequence entropy does not account for the state order in the sequence.For instance, sequences 7 and 9 have the same maximal normalized entropy of 1.

Composite complexity measures
The previous measures are based either on the sequencing or on the durations.We look now at composite measures that account simultaneously for those two aspects.

Turbulence
The turbulence T (x) of a sequence x is a composite measure proposed by Elzinga (Elzinga and Liefbroer 2007) that accounts for the number ϕ(x) of distinct subsequences of the DSS sequence and the variance s 2 t (x) of the consecutive times t j spent in the ℓ d (x) distinct states.The formula is is the maximum value that s 2 t (x) can take given the total duration ℓ(x) = j t j of that sequence.This maximum is 2 where t(x) is the mean consecutive time spent in the distinct states.
From a prediction point of view, the higher the differences in state durations and hence the higher their variance, the less uncertain the sequence.In that sense, small duration variance indicates high complexity.
The vector containing the turbulences of the sequences in a sequence object is obtained with the seqST() function.

Complexity index
The complexity index, introduced in Gabadinho, Ritschard, Studer, and Müller (2010), is a composite measure that combines the number of transitions in the sequence with the longi-tudinal entropy.It reads where h max is the theoretical maximum value of the entropy given the alphabet; i.e., h max = log a.We get the vector of complexity indexes with the seqici() function.
The minimum value of 0 can only be reached by a sequence with a single distinct state; i.e., with no transition and an entropy of 0. C(x) reaches its maximum 1 if and only if the sequence x is such that i) x contains each of the states in the alphabet, ii) the same time ℓ(x)/a is spent in each state, and iii) the number of transitions is ℓ(x) − 1.

Complexity index versus turbulence
It is instructive to look at how the turbulence and complexity indexes behave for the examples in Figure 9.The turbulence produces significantly higher values for sequences 3 and 4, which have a rather low 'sequencing' complexity but a null variance of their state durations.Indeed, this variance does not account for states that are not visited, which tends to give high turbulence values to seemingly simple sequences such as sequence 3 with two spells of same length and hence a null variance of their durations.Similarly, the turbulence exceeds the complexity index for sequences 3, 4, 5, and 7, which all have a zero variance in duration and, hence, a relatively high turbulence value.The longitudinal entropy that intervenes in the complexity index is another way of looking at the time spent in the states.It accounts, on its side, for all states, including the nonvisited ones, and, therefore, discriminates clearly between the sequences with zero duration variance.

Measuring sequence (dis)similarity
We examine now how we can measure the dissimilarity between two state sequences.As we will see in Section 9, once we have pairwise dissimilarities we will be able to run many types of powerful classical and specific statistical analysis methods on sequence data.
Many sequence dissimilarity measures have been proposed in the literature, of which the most popular in social sciences is the optimal matching (OM) edit distance.TraMineR offers a general seqdist() function that can compute the OM dissimilarity as well as a set of other dissimilarity measures.Table 6 lists the available methods and their required parameters.
The seqdist() function can output the matrix of pairwise dissimilarities or the vector of distances to a provided reference sequence.We can also compute multichannel dissimilarities (Pollock 2007) with the seqdistmc() function.
Dissimilarity measures can be classified into measures based on the count of matching attributes and those defined as the (minimal) cost of transforming one sequence into the other.Another interesting distinction is between those that make position-wise comparisons; i.e., that do not allow shifting a sequence or part of it, and those accounting for similar shifted patterns (see Table 6).Without shift, x = ABAB and y = BABA are very distant, while they are quite similar if we shift y by just one position.

Dissimilarities based on counts of common attributes
Let A(x, y) be a count of common attributes between sequences x and y.It is a proximity measure since the higher the counts, the closer the sequences.We derive a dissimilarity measure from it through the following general formula where d(x, y) is the distance between objects x and y.The dissimilarity is maximal when A(x, y) = 0; i.e., when the two sequences have no common attribute.It is zero when the sequences are identical, in which case we have A(x, y) = A(x, x) = A(y, y).Let us briefly describe the implemented count-based dissimilarity measures.
The simple Hamming distance (Hamming 1950) is the number of positions at which two sequences of equal length differ.It can equivalently be defined as ℓ − A H (x, y), with ℓ = |x| = |y| the common sequence length and A H (x, y) the number of matching positions. 10We get the Hamming distance with Equation (1) by using A H (x, y)/2 as proximity measure.
We obtain another simple distance measure by using the length A P (x, y) of the longest common prefix (LCP) between two sequences; i.e., by counting the number of successive common positions starting from the beginning of the sequences11 (see for instance Elzinga 2007b).The reversed longest common prefix (RLCP) or longest common suffix is similar to the LCP but looks for the common elements from the end rather than from the beginning of the sequences.
Another implemented metric is based on the length A S (x, y) of the longest common subsequence (LCS).12Notice that consecutive states in the common subsequence are not necessarily consecutive in the compared sequences.For example, the length of the LCS between sequences 1 and 3 of mvad.seq(see page 11 and Figure 2 on page 13) is 12 and we get 59 between sequences 2 and 5.
Quite obviously, we can only have A S (x, y) ≥ A P (x, y); i.e., the length of the LCS cannot be smaller than the length of the LCP, and hence the LCS distance cannot be greater than the LCP distance.We have also A S (x, y) ≥ A H (x, y).When compared with metrics based on position-wise counts such as the simple Hamming and the LCP distances, the LCS metric reduces distances by accounting for non-aligned matches; i.e., position-shifted similarities.

Edit distances
An edit distance is defined as the minimal cost of transforming one sequence into the other.This cost depends, indeed, on the allowed transforming operations and their individual costs.
Basically, two types of operations are considered: i) the substitution of one element by an other one, and ii) the indel; i.e., the insertion or deletion of an element, which generates a one-position shift of all the elements on its right.The generalized Hamming (HAM) and dynamic Hamming distances (DHD) (Lesnard 2006) accept only substitutions and hence no shift.The former assumes position-independent substitution costs while the second allows for position-dependent costs.The Optimal Matching (OM) distance, first considered by Levenshtein (1966) and popularized in the social sciences by Abbott (Abbott and Forrest 1986), accounts for both operations.

Setting indels and substitution costs
Usually the indel cost is set as a constant independent of the concerned position and state.Setting a high indel cost relatively to substitution costs favors substitutions while low values favor indels.We can prohibit shifts by setting the indel cost sufficiently high. 13ubstitution costs are generally organized in matrix form.A three-dimensional matrix is necessary in the case of position varying costs as used, for instance, by the DHD metric.
In the time invariant case, the substitution-cost matrix is a square symmetrical matrix of dimension a × a, where a is the number of distinct states in the alphabet.The element (i, j) in the matrix is the cost of substituting state s i with state s j .The user can either specify its own substitution-cost matrix,14 or compute one by means of the seqsubm() function with option method = "CONSTANT" or method = "TRATE".With "CONSTANT", all costs are set as the user provided cval constant (2 by default).With "TRATE", the costs are determined from the estimated transition rates as where p(s i | s j ) is the probability of observing state s i at time t + 1 given that state s j has been observed at time t (see page 17).The idea is to set a high cost when changes between s i and s j are seldom observed and lower cost when they are frequent.
Here is how we get the time-invariant transition-rate-based substitution cost matrix for the mvad data: R> scost <-seqsubm(mvad.seq,method = "TRATE") R> round(scost, 3) EM FE HE JL SC TR EM 0.000 1.971 1.987 1.957 1.988 1.961 FE 1.971 0.000 1.993 1.977 1.991 1.993 HE 1.987 1.993 0.000 1.997 1.981 1.999 JL 1.957 1.977 1.997 0.000 1.992 1.976 SC 1.988 1.991 1.981 1.992 0.000 1.995 TR 1.961 1.993 1.999 1.976 1.995 0.000 The minimum cost is 0 for the substitution of each state by itself, and the maximum is less than 2; i.e., the value that we would get for a transition not observed in the data.In accordance with what we observed in the transition rate matrix (page 18), we get the lowest costs for substituting EM (employment) to JL (joblessness) or TR (training).Remember, however, that-unlike transition rates-the substitution costs are symmetric and hence we have the lowest cost for changing JL or TR into EM.

Implemented edit distances
Generalized Hamming (HAM) and Dynamic Hamming (DHD) dissimilarities are intended for sequences of equal lengths only.The former generalizes the basic Hamming distance by allowing for state dependent substitution costs.Indeed, the count of nonmatching positions is the cost of substituting a state at each position when all costs are set to 1. DHD is the extension proposed by Lesnard (2006) to account for time-varying costs.For the mvad data set, the flexibility in substitution costs allowed by the DHD metric has only a limited impact as can be seen in Figure 10.
In addition to substitutions, Optimal Matching (OM) allows also insertions/deletions.It thus applies to sequences of unequal lengths.Since the cost of the transformation may vary with the order of the successive indels and substitutions, OM is defined as the minimal costin terms of insertions, deletions and substitutions-of transforming one sequence into the other one.The cost minimization is achieved through dynamic programming, the algorithm implemented in TraMineR being essentially that of Needleman and Wunsch (1970) with standard optimizations.
The 712 × 712 pairwise distance matrix for our mvad data computed by using the transitionrate-based costs and an indel cost of 1 is obtained with the command: R> mvad.om<-seqdist(mvad.seq,method = "OM", indel = 1, sm = scost) The mvad.om distance matrix requires only 3.96 megabytes memory space.However the number n of sequences in the data can be an important issue when computing dissimilarity matrices since both the computing time and the size of the resulting matrix increase exponentially with n.If necessary, we can divide the size by 2 by requesting only the upper triangle of the matrix with the full.matrix=FALSEargument.Most R functions accept the resulting upper-triangle objects as dissimilarity argument.

Comparing dissimilarity measures
Choosing a dissimilarity measure and setting substitution and indel costs is an important step in sequence analysis.Though popular in social sciences, distances based on such costs have raised questions in the literature (see for instance Dijkstra and Taris 1995;Wu 2000;Elzinga 2007b).The meaning of the substitution costs, their required symmetry and the sensitivity of the results to the chosen values have been pointed out as important issues.More recently, the meaning of indels was also addressed (Hollister 2009;Lesnard 2010).Favoring insertions and deletions reduces the importance of time shifts in the comparison, while favoring substitutions gives more importance to position-wise similarities.
Comparing the results obtained with various settings can also be useful for selecting the appropriate measure.Figure 10 compares the discussed dissimilarity measures for the distance to the most frequent sequence on the mvad data.We observe that, apart from the LCP metric, the measures yield very similar results.The few significant differences between HAM (or DHD) and LCS (or OM) illustrate how LCS and OM reduce dissimilarity by allowing for shifts in the comparison of the sequences.The mean difference between OM, obtained with costs derived from substitution rates, and LCS is only 0.4% of the maximal distance.The largest difference is 0.63%.These small differences are a consequence of low transition rates which lead to substitution costs comprised between 1.96 and 2; i.e., close to 2. With a constant substitution cost of 2 and an indel cost equal to 1, OM is just LCS (Elzinga 2007b).

Normalized distances
When dealing with sequences of different lengths, we may want to normalize the distances to account for these differences.More specifically, the aim of normalization is to relativize distances such that a dissimilarity of say 10 between sequences of length 100 becomes less important than a dissimilarity of 10 between sequences of length 5.While the maximal distance between a pair of sequences depends on their length, normalization aims at setting it to 1 or, at least, to a value that does not depend on the lengths.15 With seqdist() we control normalization by means of the norm argument.When setting it to TRUE, the normalization applied is determined by the selected metric.For LCP, RLCP and LCS, we apply Elzinga (2007b)'s normalization.It works as follows.Letting A(x, y) be the (non normalized) proximity measure, we first normalize this similarity The normalized distance is, then, just the complement to 1 of the normalized similarity.
which gives values comprised between 0 and 1.
For the OM distance, as well as for HAM and DHD, we apply Abbott's normalization, which consists of dividing the distance by the length of the longest of the two sequences It results that for OM with an indel cost of 1 and a constant substitution cost of 2, the maximal normalized OM distance is 2. Though OM is in this latter case equivalent to LCS, their normalized values differ.
We can also force the normalization method by specifying either "gmean" for Elzinga's normalization or "maxlength" for Abbott's solution.Alternatively, we can use "maxdist", which consists of dividing each distance by its maximal theoretical value.For LCP and LCS distances, the maximal possible value is the sum ℓ x +ℓ y of the lengths of the two sequences x and y, for HAM it is the length ℓ of the sequences, while the maximum theoretical OM distance is where c I > 0 is the indel cost, max(S) the greatest substitution cost and |ℓ x − ℓ y | the absolute value of the difference in lengths of the two sequences.With the unit indel cost and the scost transition-rate-based substitution cost matrix, this yields 139.94 for the mvad data.This is very close from twice the sequence length; i.e., 2 • 70 = 140.
It is worth mentioning that the triangle inequality property of the original distance may in some cases be lost through the "maxlength" and "maxdist" normalizations.

Dissimilarity based sequence analysis
Beside information on the similarity between any pair of sequences, a distance matrix opens access to many classical statistical and data analysis tools.It permits for instance to extract representative sequences such as medoids, to run any clustering technique based on pairwise dissimilarities and to apply multidimensional scaling.It even permits to compute pseudovariances and run ANOVA-like analyses as explained in Studer et al. (2011).We demonstrate in this section how the mvad.omdissimilarity (distance) matrix obtained with the command shown page 27 can be exploited for further statistical analysis.

Representative sequences
A major concern when analyzing sets of categorical sequences is to find useful ways of summarizing them.Possible solutions could be to determine some central or typical sequence such as the modal-most frequent-sequence or the medoid-most central-sequence.However, such solutions are of limited interest since they provide usually only a too rough idea of the main patterns in the set.A more general approach consists in finding sets of representatives and TraMineR provides the versatile generic seqrep() function for extracting such sets from the dissimilarity matrix.The function allows control over the amount of information that the representative set should convey.The sets returned by seqrep() exhibit, thus, the key features of the whole set they are extracted from, which proves useful, for example, when labeling clusters of sequences.
The principle of the search algorithm (Gabadinho, Ritschard, Studer, and Müller 2011) is to sort the sequences according to a representativeness criterion16 and to remove the redundancy by browsing the sorted sequences.The redundancy threshold is set as a percentage (10% by default) of the maximum theoretical dissimilarity D max between two sequences and the representative set will thus not contain any pair of sequences that are nearer each other than this threshold.The size of the representative set can be controlled by fixing either the minimal expected coverage of the representative set or the number nrep of representatives.
The coverage of a representative sequence is the percentage of sequences that are in its neighborhood; i.e., the number of sequences with a distance to the representative less than a selected threshold.17The total coverage of the representative set corresponds to the percentage of the n original sequences that have a representative in their neighborhood.A series of other individual and global measures to evaluate the quality of the obtained representatives is also computed.
The list of representative sequences is obtained by printing the outcome of seqrep() and we get the quality measures with the summary() method.The seqrplot() function generates representative sequence plots.

Example 1: Medoid and the centrality criterion
A first simple example18 of a representative sequence is the medoid of a set of sequences.
The medoid is the most central object; i.e., the one with minimal sum of distances to all other objects in the set (Kaufman and Rousseeuw 2005).It is a special case of representative se-quence obtained by selecting the centrality sorting criterion (criterion="dist") and setting the size of the representative set to 1 (nrep=1).
R> medoid <-seqrep(mvad.seq,diss = mvad.om,criterion = "dist", + nrep = 1) R> print(medoid, format = "SPS") The medoid of a set of sequences usually yields poor coverage.To increase coverage we should allow for more than one representative.When seeking for more than one representative, an initial sort of the sequences according to the density of their neighborhood yields better results.Neighborhood density is, therefore, the default criterion used by seqrep().
The command below finds and plots the representative set that, with a neighborhood radius of 10% (default pradius value), covers at least 25% (default coverage value) of the sequences in each of the two gcse5eq groups: R> seqrplot (mvad.seq,group = mvad$gcse5eq,diss = mvad.om,border = NA) In the resulting plot (Figure 11) the selected representative sequences are plotted bottomup according to their representativeness score with bar width proportional to the number of sequences assigned to them.At the top of the plot, two parallel series of symbols standing each for a representative are displayed horizontally on a scale ranging from 0 to the maximal theoretical distance D max .The location of the symbol associated with the representative, r i , indicates on axis A the discrepancy within the subset R i of sequences assigned to r i and on axis B the mean distance to the representative.
We learn from the plots that respectively five and one representatives are necessary for each of the two groups to achieve the 25% coverage and that the actual coverage is 29% in both cases.

Clustering sequences
Clustering is an exploratory data analysis method aimed at finding automatically homogeneous groups or clusters in the data (Kaufman and Rousseeuw 2005).In life course studies (e.g., McVicar and Anyadike-Danes 2002;Widmer and Ritschard 2009), the method has typically been used in combination with OM distances to identify distinct groups of sequences with similar patterns; that is, to define a typology of sequences.We already showed, in Section 2, how we can make a cluster analysis of sequences using the cluster library (Maechler, Rousseeuw, Struyf, and Hubert 2005).We used agnes() to make a hierarchical clustering with the Ward method, but pam() (partitioning around medoids) or diana() (divisive analysis), for example, could also be used.The four clusters solution was retained after examining the dendrogram (Figure 12) of the clustering tree obtained with: R> plot(clusterward, which.plots= 2, labels = FALSE) Figure 13 on page 34 obtained with the command below shows the representative sequences by cluster, complementing the plots of the transversal state distributions shown in Figure 1 page 7. The threshold for the coverage of the representative set is set to 35% using the coverage=.35argument.
R> seqrplot (mvad.seq,group = cl4.lab,diss = mvad.om,coverage = 0.35,+ border = NA) Looking at these two figures helps interpreting and labeling the clusters.They show that clustering from the OM distances identifies four distinct patterns of school to work transitions.In the first cluster the trajectories are clearly oriented toward early transition to employment, with, in some cases, a spell of training.The second cluster is dominated by trajectories containing a spell of school or further education followed by higher education.Cluster 3 corresponds to slow transition to employment with first an important spell of further education.In the last cluster, the transitions from school to work are more chaotic with frequent spells of training and joblessness.

Figure 3 :
Figure 3: Unsorted and sorted full-sequence index plots.

Figure 5 :
Figure 5: Mean time spent in each state by father's unemployment status.

Figure 6 :
Figure 6: Transversal state distributions by end of compulsory school qualification group.

Figure 7 :
Figure 7: Modal state sequence by end of compulsory school qualification group.

Figure 10 :
Figure 10: Distances to the most frequent sequence obtained with various metrics, mvad data.

Figure 11 :
Figure 11: Representative sequences by end of compulsory school qualification group.

Table 1 :
TraMineR's key functions.theirmonthlyfollow-up over the course of 6 years starting in the month where they were first eligible to leave compulsory education (July 1993).Each individual is characterized by a unique identifier, 13 covariates and 72 monthly activity state variables from July 1993 to June 1999.Since the first two months of the follow-up are summer holidays, we look hereafter at trajectories from September 1993 yielding sequences of 70 monthly statuses.The states are school, FE (further education), employment, training, joblessness, and HE (higher education).See Table2for a description of the variables in mvad.
It contains the data used by McVicar and Anyadike-Danes (2002) for studying the school-towork transition in Northern Ireland.The figures cover 712 individuals, the sequences being

Table 2 :
List of variables in the mvad data set.
Table 4 displays the same example

Table 4 :
Sequence data representations; some formats handled by the seqformat() function.

Table 5 :
Main sequence object attributes.
We show below how to create a state sequence object from the mvad data

6
However, unlike for graphical displays, functions returning statistics and sequence characteristics do not have a group argument.We can retrieve the values by levels of a covariate with the row indexing mechanism 7 or with the by() function: Figure 8: Transversal entropy by end of compulsory school qualification group.

Table 6 :
List of available metrics for computing distances with the seqdist() function.