**Table 3** Prime implicants for the presence of the outcome, by context and time point

Source: authors' elaboration on NEPS data. Conservative solutions. Consistency threshold: 0.75. Upper-case letters represent the presence and lower-case letters the absence of the condition *m* = months

7Applying consistency thresholds between 0.76 and 0.85 produces results identical to those displayed in Table 3, with the exception of the second prime implicant, which in this case consists of "child\*high-paredu\*PARTNERSHIP."

8In Boolean language, upper-case letters represent the presence and lower-case letters the absence of the condition; the star sign represents the logical AND (set intersection) and the plus sign represents the logical OR (set union).

combinations of factors sufficient for remaining in employment or in education diverged across the two contexts: In East Germany, the only prime implicant from 24 years on is constituted by women in a partnership *and* with children. In West Germany, where the transition to fertility was generally delayed compared to the East (Kreyenfeld 2004), for 24-year-old women we still find the alternative combinations of not having children *and* having highly educated parents *or* or not having children *and* being in a partnership. However, from age 31, the only prime implicant left is that of women without children and with highly educated parents.

For both contexts, we observe a collapse over time into a single prime implicant. However, the implications are very different for East and West Germany: The prime implicant for East Germany (being in a partnership *and* with children) corresponds to the prevailing configuration of attributes for women aged 24 and older. Moreover, for these Eastern German women, the outcome is skewed to the right or, in other words, the overwhelming majority were in employment or in education (consistently with what shown by Huinink et al. 1995). Hence, the prime implicant and the outcome are in a relation of almost perfect set coincidence, which should not be interpreted as causal but rather as constitutive (Borgna 2013). In substantive terms, as noted earlier, this is driven by a strong standardization of individual life courses (Fasang and Raab 2014; Struffolino et al. 2016). In West Germany, in contrast, the lasting prime implicant corresponds to a limited portion of the population and can therefore properly be read as "outcome-enabling": Among women older than 24, only those without children *and* with highly educated parents could systematically still be found in employment or in education. The prime implicant of women in a partnership *and* without children disappears, because for ages 31 and 35 it corresponds to a group that is too small to be considered an existing configuration, because the great majority of women in a partnership eventually had children, and those who did not have any originate from families with highly educated parents (see also Table 2).

Overall, while we detected a greater heterogeneity in the life courses of women in West Germany from the mere observation of the truth tables, the truth-table analyses suggest that this greater diversity does not translate into a greater variety of pathways towards employment or education. While further research would be needed for any definitive, substantial conclusions, the results of this empirical exercise are not only in line with previous evidence on women's labor-market participation in East and West Germany but also provide some new insights into the different degree of complexity of their life courses.

#### **4 Concluding Remarks**

In this contribution, we have put forward a sequential mixed-methods design to unpack the temporal and configurational complexity of social phenomena. By applying SA, researchers can analyze the unfolding of phenomena over time in terms of sequencing of states and duration of events. More specifically, through discrepancy analysis of state sequences, researchers can identify the crucial points in time when trajectories (start to) diverge from each other. In a second research stage, QCA can be applied to investigate whether, at such crucial time points, some particular combinations of factors explain why the cases display a given state (or outcome). By analyzing these cross-sectional turning points sequentially, we preserve the longitudinal and holistic perspective on trajectories.

We have illustrated the usefulness of our framework with an empirical analysis of employment trajectories in divided Germany. We have shown that, especially for West Germany, gender significantly explains the divergence of trajectories at what we characterized as turning points in women's life course. Moreover, while a plurality of factors initially explain women's presence in employment or education in both contexts, only one configuration of factors for each context is systematically linked to the outcome at the end of the period considered. Hence, by combining SA and QCA, we were able to unveil the dynamic of conditions sufficient for women to be in employment or education.

Notwithstanding the specificities of our empirical illustration, we believe that our proposed framework holds the promise of being useful to a wide range of researchers from a variety of disciplinary fields and who are interested in various kinds of micro-, meso-, and macro-level processes that involve a temporal and configurational dimension. We limit ourselves to three fields for which this SA-QCA design clearly appears to be valuable. First, recent developments in social stratification research have emphasized the role of cumulative advantage (DiPrete and Eirich 2006) and intersectionality (Platt 2011) but struggle to combine the two aspects in a single analytical framework. Second, a holistic perspective on both the context of implementation and the process of change is central in policy analysis and evaluation, especially for the approach of theory-based evaluation (Befani and Sager 2006). Third, the definition of "critical junctures" in the historical institutionalist tradition (Thelen 1999) appears to bear more than one similarity to our treatment of outcome-enabling configurations at turning points along the sequences. Similarly, Schneider and Rohlfing (2013) and more recently Williams and Gemperle (2017) have called for an integration of QCA and process-tracing methodologies to identify the "situational" causal mechanisms. Future research could explore the potential of combining our SA-QCA framework with more qualitative methods, like process tracing to the study of social processes and critical junctures in particular.

**Acknowledgements** We thank Anette E. Fasang for generously sharing the routines for data preparation and providing helpful comments and suggestions on an earlier version of this chapter. For insightful comments, we are also grateful to David Brady, Markus Siewert, Eva Thomann, as well as to the participants in the Colloquia *Quantico* at the ISW-Humboldt University of Berlin and *CO:STA* at the WZB Berlin Social Science Center. We thank Ana Santiago-Vela for excellent research assistance.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Combining Sequence Analysis and Hidden Markov Models in the Analysis of Complex Life Sequence Data**

**Satu Helske, Jouni Helske, and Mervi Eerola**

#### **1 Introduction**

Longitudinal data often consists of multiple parallel sequences that ought to be analyzed jointly. For example, life course data may contain sequences of employment, family formation, and residence. Such data is often referred to as multichannel or multidimensional sequence data. A multichannel approach often gives a simpler representation of the data as opposed to combining states across life domains (the extended alphabet approach); the latter approach rapidly grows the state space as the number of channels and/or states grows. If some data is only partially observed, the multichannel approach also allows for handling data as it is instead of having to make difficult decisions on how to combine observed and unobserved states (Helske and Helske 2018).

Joint analysis of complex multidimensional data poses several challenges. Multichannel sequence analysis (Gauthier et al. 2010) has been the standard tool for the analysis of multichannel sequence data (for empirical applications see, e.g., Eerola and Helske 2016; Müller et al. 2012; Spallek et al. 2014). This approach is

S. Helske (-)

Institute for Analytical Sociology, Linköping University, Linköping, Sweden

Department of Sociology, University of Oxford, Oxford, UK

J. Helske

Department of Science and Technology, Linköping University, Linköping, Sweden

Department of Mathematics and Statistics, University of Jyvaskyla, Jyvaskyla, Finland

M. Eerola Centre of Statistics, University of Turku, Turku, Finland

Department of Mathematics and Statistics, University of Jyvaskyla, Jyvaskyla, Finland e-mail<EMAIL_ADDRESS>
simple and fast in computing dissimilarities between sequences, and cluster analysis is often used for grouping similar sequences. Describing and visualizing results is, however, often challenging.

We propose an approach for compressing the information within multichannel sequences and for facilitating the interpretation of such data by finding (1) groups of similar trajectories and (2) similar phases within trajectories belonging to the same group. For the first task we use the standard multichannel sequence analysis approach and for the second task we propose using hidden Markov models (HMMs). With the help of HMMs the data can then be illustrated with a graph showing typical phases within trajectories and the transitions between them and/or shown as simplified (single-channel) trajectories consisting of these typical phases. We illustrate this approach with an empirical application to complex longitudinal life course data but such an approach, and HMMs in general, are useful in various longitudinal problems across disciplines.

Hidden Markov models have been widely used in economics, bioinformatics, and engineering (see, e.g., MacDonald and Zucchini 1997; Durbin et al. 1998; Rabiner 1989), often to study single long sequences such as time series. In social sciences, such models are commonly referred to as latent Markov (chain) models (Wiggins 1955, 1973; Van de Pol and De Leeuw 1986); typically they have been used for analysing panel data with a few measurement points. In the social science framework, Vermunt et al. (1999) extended the HMM to include individual covariates and Bartolucci et al. (2007) further developed it for multichannel observations. See also Taushanov and Berchtold (2018) in this bundle.

Hidden Markov modelling have been applied in various longitudinal settings; for accounting for measurement error and unobserved heterogeneity (e.g., Van de Pol and Langeheine 1990; Poulsen 1990; Breen and Moisio 2004; Vermunt et al. 2008; Pavlopoulos and Vermunt 2015), for finding latent sub-populations (e.g., Van de Pol and Langeheine 1990; McDonough et al. 2010; Bassi 2014), and for detecting true unobservable states (e.g., various periods of the bipolar disorder in Lopez 2008).

To the best of our knowledge, few papers apply HMMs to multichannel social sequence data and they all consider binary observations. Bartolucci et al. (2007) studied criminal trajectories of 11,400 offenders, applying HMMs to ten-channel data with six time points. Ip et al. (2015) analysed and classified 18-item profiles of food security among 248 Latino farm worker households in the USA for eight time points. Rijmen et al. (2008) studied 12 parallel trajectories of emotions at 63 time points among 32 anorectic patients. Our analysis extends this framework into multichannel data with much longer and multinomial sequences.

The rest of the paper is structured as follows. We start by giving an introduction to HMMs (we assume that the reader is familiar with sequence analysis and refer to the introduction chapter in this book for the less experienced). We then proceed to framing our goals in the context of complex life course data. We continue by describing the data and the empirical analysis and show the results. We conclude with discussing the usefulness of the method, the challenges it poses, and mention some future directions.

#### **2 Hidden Markov Model**

In hidden Markov models, observations are related to a hidden process following a Markov chain. Hidden states can only be detected through the observed sequence(s), as they generate or "emit" observations on varying probabilities.

Let us assume we have multichannel sequence data with *N* individuals, *T* timepoints, and *C* channels and a hidden Markov model with *S* hidden states. Now **z***i* = *(zi*1*, zi*2*,...,ziT )* represents the hidden state sequence for individual *i* = 1*,...,N* from time 1 to time *t* and *yitc* denotes the observation of individual *i* at time *t* = 1*,...,T* in channel *c* = 1*,...,C*.

Figure 1 illustrates the structure of an HMM for two-channel data. The first order Markov assumption states that the probability of transitioning to the hidden state at time *t* only depends on the hidden state at the previous time point *t* −1. Here we also assume the same latent structure applies to all channels, i.e., hidden state *zit* emits observed states *yitc* in all channels *c* and observations *yit*1*,...,yitC* are assumed conditionally independent given the hidden state *zit* .

The following probabilities characterize a discrete first-order hidden Markov model for multichannel data:

• *Initial probability* vector *π* = {*πs*} of length *S*, where *πs* is the probability of starting from the hidden state *s*:

$$
\pi\_s = P(z\_{l1} = s); \quad s \in \{1, \ldots, S\}.
$$

• *Transition probability* matrix *A* = {*asr*} of size *S*×*S*, where *asr* is the probability of moving from the hidden state *s* at time *t* − 1 to the hidden state *r* at time *t*:

$$a\_{sr} = P(z\_{lt} = r | z\_{l(t-1)} = s); \quad s, \ r \in \{1, \ldots, S\}.$$

• *C emission probability* matrices *Bc* = {*bs(mc)*} of size *S* × *Mc*, where *bs(mc)* is the probability of the hidden state *s* emitting the observed state *mc* in channel *c* and *Mc* is the number of observed states in channel *c*:

$$b\_s(m\_c) = P(\mathbf{y}\_{ltc} = m\_c | z\_{lt} = s); \quad s \in \{1, \dots, S\}, \ m\_c \in \{1, \dots, M\_c\}.$$

**Fig. 1** Illustration of hidden and observed state sequences in a hidden Markov model for twochannel data of individual *i*. The hidden state at time *t* is illustrated with *zit* inside a circle and the observed state at time *t* in channel *c* with *yitc* inside a rectangle. Arrows indicate dependencies between states

Typically, the maximum likelihood estimates of these probabilities are calculated with the Baum–Welch algorithm, i.e., the expectation–maximization (EM) algorithm for HMMs (Baum and Petrie 1966; Rabiner 1989). The most probable path of hidden states for each subject given their observations and the model can be computed using the Viterbi algorithm (Viterbi 1967; Rabiner 1989). Missing observations are handled straightforwardly. When observation *yitc* is missing, it does not contribute to the estimation of model parameters nor hidden states. See Helske and Helske (2018) for a more extensive presentation on HMMs for multichannel data.

#### **3 Combining Sequence Analysis and Hidden Markov Models for Complex Life Sequences**

For analysing complex life sequence data, we aim to compress the information into two types of components:

The first component corresponds to finding clusters or latent classes of individuals who have experienced similar life events in similar order and timing. The other, time-varying components should correspond to life stages during which individuals are more likely to have similar experiences, e.g., observed states within the sequences. These life stages could be either stable episodes between two transitions (e.g., employed and married without children) or characterized by transitions in some of the life domains (e.g., moving between unemployment and short-term jobs). Individuals may, and typically do, go through several different life stages during their life course.

SA followed by cluster analysis is a typical strategy for grouping life trajectories. Hidden Markov models, in turn, may be used for finding time-varying latent structures and transitions between them. At first, we use multichannel SA to compute pairwise dissimilarities and then group individuals into clusters. Separate HMMs are then fitted for each cluster. The number and nature of the hidden states are determined independently for each group.

We estimate left-to-right HMMs where transitions to previous hidden states are not allowed. We had several reasons to do this. First, left-to-right models are simpler to estimate since some of the parameters are restricted to zero. Second, due to the nature of the life trajectories, also the observed states tend to show a left-to-right behaviour and many of the HMMs would end up being estimated close to left-toright models anyway. Third, we find that left-to-right models are often easier to interpret in the context of life course: individuals go through different life stages but even if they return to have a similar life stage compared to a previous one – say re-marriage after a divorce – this second life stage comes with a different history compared to the first time.

#### **4 Data**

We illustrate the analysis of complex life sequence data using a subsample of the German National Educational Panel Survey (NEPS) (Blossfeld et al. 2011). We restricted the analysis to the life courses of an age cohort born in 1955–1959. Only individuals who were born in Germany or moved there before the age of 14 are included.

The data consisted of monthly life statuses of 1731 individuals in three life domains (labour market participation, partnerships, and parenthood) from age 15 to age 50. For each individual, there were three parallel sequences of length 434, which made altogether 2,253,762 data points (of which 2,232,730 were observed and 21,032 were missing). Using the monthly time scale also allowed for the detection of smaller fluctuations in life courses, e.g. recurrent transitions between short-term unemployment and employment.

#### *4.1 Sequences*

The sequences in three life domains were constructed as follows:

**Labour market participation** with 4 states:

#### **Partnerships** with 4 states:

#### **Parenthood** with 2 states:

The coding for parenthood was very simple. A practical reason was that this record was available for most individuals, whereas more detailed information was often missing. On the other hand, we can argue that specifically the experience of becoming a parent is relevant as one step in the developmental process into adulthood.

For the latter two life domains, the status of each month was usually determined from the latest event. An exception was made for the rare partnerships that lasted for less than a month; there separation was coded from the following month onward. In a case of multiple records per month in the career domain, the final status was given according to assumed importance: school and vocational training came before employment, which in turn dominated over vocational preparation, unemployment, and other non-employment statuses.

Altogether 306 individuals (17.7%) had some missing information in one or two life domains. Thus, at each time point we had at least some information from each individual.

#### **5 Analysis**

We have little prior knowledge on the structure of the model; hence, how many clusters to choose and how many hidden states to include in each cluster? Since the complexity of these types of life course trajectories varies a lot (e.g., some individuals have no family-related transitions while others have many), we expected the groups to have varying numbers of hidden states.

#### *5.1 Sequence Analysis and Clustering*

We started by applying multichannel sequence analysis and computed the dissimilarities between the sequences. These were then used in cluster analysis.

The dissimilarities between sequences were determined according to the generalized Hamming distance with user-defined substitution costs (see Table 1). We set the highest cost to be the same in all life domains to give them equal weight. We gave no cost for substituting missing states since we wanted to determine dissimilarity based on the observed trajectories. Regarding the costs within different life domains, our choices were mainly based on how far the states are regarded on the pathway to adulthood and, in terms of labour market participation, also on how close the other states can be regarded to employment which is often the favourable state. The metric compares observed states time point by time point and gives a cost for mismatches. It generally works well in a multichannel problem where timing is important (Studer and Ritschard 2016) and resulted in meaningful clusters with high goodness-of-fit.

We used Ward's clustering method for the Hamming dissimilarities and chose six clustering solutions with 7–12 clusters for further examination. The choice was based on goodness-of-fit statistics, the dendrogram, and the interpretability of the clusters. Ward's method was chosen because it typically produces usable and relatively even-sized clusters compared to most of the other clustering methods (Aassve et al. 2007; Helske et al. 2015). Also, the method is hierarchical (agglomerative), so when two smaller clusters are merged, all other clusters remain the same. This

**Table 1** Substitution costs for Hamming distances in three life domains: labour market participation, partnerships, and parenthood

means that among the 7 + 8 + 9 + 10 + 11 + 12 = 57 clusters in the six sets of clustering results, only 7 + 2 + 2 + 2 + 2 + 2 = 17 were unique, resulting in significant decrease in the number of models to be estimated compared to nonhierarchical clustering.

#### *5.2 Hidden Markov Models for Clusters*

At the next step, we estimated five HMMs with 4–8 or 5–9 hidden states separately for each of the 17 unique clusters—fewer hidden states for clusters with simpler observed trajectories, more for the more complex ones. Since the goal was to find general life stages between adolescence and middle age in a given group, having too few or too many hidden states was not plausible nor interpretable. When increasing the number of hidden states, at some point they lost their distinctive nature (consecutive states had very similar emission probabilities) and/or they were rarely "visited" in the most probable paths of hidden states.

A well-known problem with the HMM estimation is that most of the optimization methods are sensitive to initial estimates of the parameters. In order to reduce the risk of being trapped in a poor local optimum, we estimated the models numerous times with random starting values. We continued re-estimation until we had found the same optimum for at least 100 times (which turned out to be much more than necessary).

For each cluster, we compared the HMMs with a different number of hidden states to find the best model. Bayesian information criterion (BIC) and other information criteria are common choices for comparison of HMMs with different numbers of hidden states. Another common option for model selection is crossvalidation.

We chose to use BIC as it generally selects parsimonious models. Unfortunately, here BIC kept suggesting models with more and more states. We did, however, use BIC as one source of information for choosing the number of hidden states by looking for turning points in BIC after which additional hidden states offered little improvement. In addition to BIC, the choice of the number of hidden states was based on the interpretability of the model and the prevalence of the hidden states in the individual trajectories.

#### *5.3 Software*

Analyses were conducted with the R software (R Core Team 2015) by using the packages TraMineR (Gabadinho et al. 2011) for sequence analysis, cluster (Maechler et al. 2015) for cluster analysis, and seqHMM (Helske and Helske 2018) for hidden Markov modelling. For the estimation of HMMs we used the automatic re-estimation routine for the EM algorithm provided in the model estimation function.

#### **6 Results**

The number of hidden states per cluster varied between six and eight. The model of eight clusters resulted in the smallest BIC (even the highest likelihood) and was chosen as the best model. We present a few different ways to describe the results: a table showing the most typical transitions in each cluster, a figure illustrating the structure of the HMMs, and a figure of the most probable hidden states, i.e., the trajectories of general life stages for each individual.

Table 2 describes each cluster in terms of some important transitions and states: typical labour market participation (showing the timing of completing education and the type of employment after that), partnership histories (age at first partnership, the type and number of partnerships), and parenthood (the timing of the first child). It also shows the number of individuals in each cluster and the proportion of the sample, as well as the hidden states described with the most important transitions.

Figure 2 illustrates the HMM structure for each of the eight clusters. It shows the HMMs as directed graphs where the pies represent hidden states and the slices show the emission probabilities of observed states within each hidden state (to draw attention to the most prevalent observations, we only show probabilities that are greater than 0.05). The arrows indicate transition probabilities between the hidden states—the thicker the arrow, the higher the probability.

Figure 3 illustrates the most probable hidden state paths. We have assigned similar colours to similar hidden states across clusters.

As an example of how to interpret these figures, let us look at the smallest of the clusters titled Single parents (cluster H). All individuals start from the first hidden

**Fig.**OUunemployed;single,cohabiting,divorced/separated;child(ren).Hidden states are described by the most probable observations

state (State 1, indicated with light blue in the hidden state paths), a life stage where they are childless singles and mostly studying. For almost all, the next transition is to State 2 (dark blue), moving to employment. A few make a straight transition to State 3 (light pink), a life stage of becoming parents and being out of the workforce. State 4 (darker purple) describes a life stage during which individuals are singles, have children, and are employed. This is the most prevalent life stage for the members of this cluster and many stay there until the end of the follow-up. A few move out of employment, mostly to unemployment (State 5, light purple). During the last life stage, experienced by almost half of the members, individuals form their first partnerships (State 6, yellow).

In general, the clusters were well separated from each other by the timing and occurrence of labour market participation and family states. The two largest clusters composing of half of the respondents were characterized by (mostly) short education and family. The biggest difference was in the timing of partnership and parenthood transitions which occurred either earlier in life (cluster A) or later (cluster B). The third largest cluster (cluster C) mostly consisted of individuals, more often men, who had long education and later family transitions. Another cluster with early family transitions (cluster D) consisted of mostly women and was characterized with a long career break for mostly taking care of children.

Two clusters were characterized by no or very late parenthood. They differed in the timing of the partnerships; the larger cluster (cluster E) had earlier first partnerships while in the smaller cluster (cluster F) partnerships were delayed or omitted altogether. The two smallest clusters consisted of parents living divorced or separated (cluster H) or single parents (cluster G).

#### **7 Discussion**

When analysing complex sequence data with multiple channels, describing and visualizing the data can be a challenge. By combining sequence analysis and hidden Markov models the information in data can be compressed into hidden states (life stages) and clusters (general patterns in life courses). Hidden states were able to capture general life stages that included not only rather stable episodes such as being employed and married with children (e.g., State 7 for cluster A) but also life stages characterized by change, e.g., moving between unemployment and short-term employment (State 3 for cluster F).

We presented two different ways of HMM-based visualizations that give complementary information but could also be shown alone—it is up to the researcher to decide which one is more informative in each case. The HMM graphs show the structure of the hidden states and the transitions between them; also parameter estimates could be easily included in the graph. The most probable paths of hidden states show individual-level information on the approximate prevalence and timing of different life stages.

Despite its usefulness as a data reduction technique, this approach comes with some challenges. A major one is the estimation of several HMMs when the number of hidden states and clusters is unknown. For these challenges, we used a few approaches. In terms of the number of clusters, we used a hierarchical clustering method which reduced the number of models to be estimated compared to nonhierarchical clustering. We then estimated a single model numerous times with randomized starting values to find the one with the highest likelihood, using parallel computation for improved efficiency.

Another issue is that we take the SA clusters as fixed. In reality, there is, of course, a lot of uncertainty which we do not take into account. Also, we do not discuss other trajectory grouping techniques besides SA. To our knowledge, there are not many methods suitable for multichannel sequence data; we experimented with latent class analysis (Collins and Wugalter 1992) which did not lead to satisfactory results. On the other hand, regarding the parameter uncertainty, in theory it is possible to compute asymptotic standard errors from the Hessian matrix obtained from the numerical optimization algorithms, but in practice the underlying asymptions are typically not met (Zucchini and MacDonald 2009).

The mixture hidden Markov model (MHMM) offers a solution to the problem of uncertainty of clustering. In the MHMM, instead of fixing individuals to the clusters defined during the SA step, we could use all data to estimate a mixture of HMMs where each individual belongs to each cluster with some probability (preferably with a large probability for one cluster and a small probability to all others). In a complex setting, SA can be used to determine the range of potential clustering structures. It can also be of aid when setting initial values for the estimation process, which is often essential when using very large models.

Although in theory the MHMM approach allows even more flexibility to the modelling and potential for more interesting ways of inference, there are some practical computational problems in the MHMM methodology. The parameter estimation of HMMs is often very sensitive to initial values, and the computational costs increase rapidly when the number of hidden states grows. These problems are even more prominent in complex MHMMs, especially when the structure of the model (in terms of the number of hidden states and/or clusters) is not known. For this study, we were not able to find stable solutions for MHMMs despite large computational resources available—the multichannel structure, long sequences, and the relatively large number of individuals in our data was too challenging a combination for parameter estimation. Nevertheless, the MHMM can be useful in other settings. It has been successfully used for simpler problems, e.g., for accounting for measurement error and unobserved heterogeneity.

An extention not covered in this paper is the inclusion of external covariates. Personal characteristics and other relevant factors, time-constant as well as timevarying, could be used to predict transition probabilities between life stages. In MHMMs, time-constant covariates may also be used to predict cluster memberships. See, e.g., Vermunt et al. (2008) for a general presentation of such models.

We are currently studying algorithmic variations which can reduce the computational complexity of the MHMM estimation. Further research is also needed regarding model selection and the goodness-of-fit of left-to-right HMMs and MHMMs. Further theoretical and empirical studies are needed for detecting the reasons for the failure of BIC and for discovering selection criteria that are better suited for finding parsimonious HMMs.

Another topic for future research is the potential of hidden Markov models and Markovian models in general as mechanisms of generating social sequence data.

The aim of our study was to describe complex life sequence data and for that goal, the SA-HMM approach gave satisfactory results in a reasonable time. We were able to find meaningful and well-separating clusters and to visualize their complex life course information by using HMM graphs and the most probable paths of life stages for each individual.

**Acknowledgements** This paper uses data from the National Educational Panel Study (NEPS) Starting Cohort 6–Adults (Adult Education and Lifelong Learning), doi:10.5157/NEPS:SC6:3.0.1. From 2008 to 2013, the NEPS data were collected as part of the Framework Programme for the Promotion of Empirical Educational Research funded by the German Federal Ministry of Education and Research and supported by the Federal States. As of 2014, the NEPS survey is carried out by the Leibniz Institute for Educational Trajectories (LIfBi).

Satu Helske is grateful for support for this research from the John Fell Oxford University Press (OUP) Research Fund and the Department of Mathematics and Statistics at the University of Jyväskylä, Finland, and Jouni Helske for the Emil Aaltonen Foundation and the Academy of Finland (research grant 284513).

We also wish to thank three anonymous referees for their helpful comments and suggestions.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Part V Advances in Sequence Clustering**

## **Markovian-Based Clustering of Internet Addiction Trajectories**

203

**Zhivko Taushanov and André Berchtold**

#### **1 Introduction**

The clustering of trajectories has gained much interest in recent years from the scientific community, especially in the social sciences, because the number of longitudinal studies, as compared to cross-sectional ones, has been constantly increasing. As regards categorical data, the most common approach relies on the Optimal Matching (OM) to compute a distance between each pair of trajectories before clustering them, whereas the Growth Mixture Model (GMM) can be applied for continuous data. However, these two approaches suffer from some shortcomings, calling for the need to develop and apply alternative approaches. For instance, OM requires the choice of a substitution cost measure and other parameters. GMM gives a lot of importance to the shape of sequences. Therefore, there is a risk to overfit the data when nonlinear trajectories are considered on quite short sequences. The other issues of GMM include computational load, presence of local optima, missing data treatment, model selection criteria, the need for large sample size, and unclear Type I error rates (Wang and Bodner 2007).

In this paper, we study the use of a specific class of Markovian Models called the Hidden Mixture Transition Distribution (HMTD) model (Bolano and Berchtold 2016) for clustering purpose. Even if this model-based approach was developed as a tool for the analysis of continuous trajectories, it also allows for their clustering without a priori knowledge of cluster membership. Moreover, covariates can be easily included in the model.

The HMTD and GMM clustering approaches are applied and compared on a dataset of trajectories of the Internet Addiction Test (IAT). Excessive Internet use,

Z. Taushanov (-) · A. Berchtold

Institute of Social Sciences and NCCR LIVES, University of Lausanne, Lausanne, Switzerland e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
G. Ritschard, M. Studer (eds.), *Sequence Analysis and Related Approaches*, Life Course Research and Social Policies 10, https://doi.org/10.1007/978-3-319-95420-2\_12

especially among youths, is an emerging health issue in the medical literature, with studies showing contrasting results. Surís et al. (2014) show a significant association between problematic Internet use and several somatic disorders, including back, weight, musculoskeletal, and sleep problems. Moreover, several chronic conditions are also significantly associated with problematic Internet use. In contrast, another study finds no significant Internet use effect on the development of overweight among youths (Barrense-Dias et al. 2015).

While several alternative approaches (e.g. Skarupova et al. 2015) have been introduced over the years, the tool most often used to quantify the degree of addiction to Internet is still the Internet Addiction Test (IAT), developed by Young (1998). However, since the test's scale is based on 20 items and is quite long, its psychometric properties are matters of controversy (Faraci et al. 2013) and the test is not considered suitable for the successive measurement of the same subjects (testretest). Its use in longitudinal contexts remains problematic because of the difficulty to distinguish between the real evolution of subjects and changes due to behavior of the IAT itself.

To gain information on the behavior of the IAT in longitudinal studies, we need to compare the typical trajectories of the repeated IAT measurements with other characteristics of the subjects under study. Thus, we first cluster the IAT trajectories into a finite set of meaningful groups and then compare these groups with the known characteristics of subjects that are either time-invariant or evolve over time. Specifically, the goals of this study are (1) to separate the Internet addiction trajectories into an optimal number of meaningful categories using HMTD, (2) to explore how does the introduction of the covariates influence the previous optimal partition, and (3) to compare the HMTD clustering with an equivalent GMM clustering in order to gain information on the respective strengths of both approaches. We hypothesize that (1) the IAT scores computed for the same person can vary considerably over time, implying that the trajectories are difficult to classify; (2) a classification using covariates is easier to interpret than a classification without any additional information on the clustered variable itself; and (3) the HMTD approach can lead to more sound and easier-to-use solutions as compared to the solutions obtained using GMM. However, we must stress that it is impossible to conclude that one method is superior to another, especially using real data, without knowing the true cluster membership. So this work must be considered as a first step in the comparison of HMTD and GMM as clustering tools.

#### **2 Data and Methods**

#### *2.1 Data*

The data we considered are from<EMAIL_ADDRESS>(Surís et al. 2012), a longitudinal study on the use of Internet among youths in the Swiss canton of Vaud (the largest canton in the French-speaking part of Switzerland). The data were collected five times every six months from Spring 2012 (*T*0, baseline) to Spring 2014 (*T*4) using an online questionnaire. The data for the first time were collected from schools during the computerlab periods. Then, the students who agreed to participate in the study were contacted again by email from *T*1 to *T*4 to answer follow-up questionnaires on their home computer. A convenience sample of *n* = 185 adolescents who answered all five questionnaires is used for the present study (67% females; mean age at *T*0: 14.1 years). For more details on the overall design of the study and data collection, (see Surís et al. 2012; Piguet et al. 2016).

The main outcome is the IAT score measured at each wave for each subject. The IAT developed by Young (1998) and validated in French by Khazaal et al. (2008) is a scale ranging from 0 to 100, based on the answers to 20 items whose possible answers range from Never (coded 0) to Always (5). Examples of items are, *How often do you find yourself staying online longer than you intended?* and *How often do you fear that life without the Internet would be boring, empty, and joyless?*

In addition to the IAT, we also considered several important characteristics of the subjects, either fixed in time—gender, age at baseline, and education track at baseline (extended requirements vs. basic requirements)—or evolving over time emotional well-being (measured by the WHO-5 index) and Body Mass Index (BMI, computed from auto-reported measures of height and weight). Note that the WHO-5 index was not evaluated on the third wave of the study, and so for the present paper, we imputed values as the simple mean between the values of the second and fourth waves. Similarly, we imputed the BMI for the second wave of the study as the mean between the values of the first and third waves.

#### *2.2 Clustering Using the HMTD Model*

We used a specific class of Markovian Models, the HMTD model, to cluster the longitudinal sequences of continuous data. This model combines a latent and an observed level (Bolano and Berchtold 2016). The visible level is a Mixture Transition Distribution (MTD) model that was first introduced by Raftery in 1985 as an approximation of high-order Markov chains Raftery (1985) and then developed by Berchtold (2001, 2003) and Berchtold and Raftery (2002). Here, we used a Gaussian version of the MTD model, where the mean of the Gaussian distribution is a function of past observations. Because of the small size of each sequence of the observed outcome (five data points, from *T*0 to *T*4), long dependencies between successive observations could not be considered, and therefore we fix the dependence order for the mean of the Gaussian distributions of each component to one:

$$
\mu\_{\mathfrak{g},t} = \varphi\_{\mathfrak{g},0} + \varphi\_{\mathfrak{g},1} \chi\_{t-1},
$$

where *ϕg,*0 is the constant for the mean for component *g* and *ϕg,*1 is the autoregressive parameter indicating the dependence from the previous observation *xt*−1. Similarly the variance of each component can be written as a function of the past periods variability: *σ*2 *g,t* = *θg,*0 + -*S s*=1 *θg,s x*2 *t*−*s*. However given the small number of time periods in our dataset, and for the sake of simplicity, we decided to treat the variance as a constant: *σ*2 *g,t* = *θg,*0.

In the HMTD model, the latent level is a homogeneous Markov chain. Each state of the chain is associated with a different Gaussian component at the visible level, with the transition matrix used to determine which component best represents the current observation. To use the HMTD model as a clustering tool, we assume the hidden transition matrix to be the identity matrix. Consequently, each sequence of successive observations is associated with only one component of the model, thus generating a clustering of sequences into mutually exclusive groups. Notice that in this case, the resulting model is no more a hidden Markov model, but a mixture of Gaussian distributions. However, it is still interesting to view it as a HMTD, because it is then possible to compare the clustering model with other models, especially with semi-clustering models whose transition matrix is not the identity matrix, but a triangular matrix letting data trajectories move from one group to another in a specific order.

In addition to the clustering based on the IAT variable only, we performed a second clustering adding information from five covariates (gender, age at *T*0, education track at *T*0, WHO-5, and BMI). These covariates are introduced as additional terms in the specification of the mean of each visible component of the model, and the categorical variables are introduced as dummy variables. We then rewrite the mean of the *g*-th component as

$$\mu\_{\mathfrak{g},\mathfrak{l}} = \varphi\_{\mathfrak{g},0} + \varphi\_{\mathfrak{g},1} \mathbf{x}\_{\mathfrak{l}-1} + \varphi\_{\mathfrak{g},2} \mathbf{Gender}(male) + \varphi\_{\mathfrak{g},3} \mathbf{Age}$$

$$+ \varphi\_{\mathfrak{g},4} \text{ Education}(extended) + \varphi\_{\mathfrak{g},5} \mathbf{WHO-\mathfrak{f}} + \varphi\_{\mathfrak{g},6} \mathbf{BMI}$$

with *female* and *basic requirements* used as reference modalities for Gender and Education, respectively.

In practice, continuous covariates are centered around the sample mean before computing the clustering model in order to allow for a better convergence of the estimation algorithm. A comparison of the two specifications of the mean, with and without covariates, illustrates whether the inclusion of covariates in the model helps to improve the clustering process. It must be mentioned that, in addition to these two HMTD models, many other specifications were tried, following a hierarchical approach (Bolano and Berchtold 2016), but none of these alternative specifications seemed to give a more useful clustering of IAT trajectories.

The HMTD model is estimated by maximizing its log-likelihood. When the variance of each component of the model is constant, the log-likelihood can be derived with respect to all parameters, but in the general case of time-varying variances (Berchtold 2003), the log-likelihood is generally not differentiable, and the solution space can be very complex. A specific heuristic is then applied to obtain the solution (Taushanov and Berchtold 2017). Since this heuristic can accommodate to all possible specifications of the HMTD model, we used it for all computations. Regarding cluster assignation, we used the standard Viterbi algorithm which is able to find the best sequence of hidden states in function of the observed data and of the current model (Forney 1973). In the specific case of clustering, the Viterbi algorithm simply assign each observed trajectory to the most likely component.

We used a bootstrap procedure to obtain confidence intervals for each parameter, but since our goal here was to validate not the initial classification itself, but the parameters associated with the model describing each visible component of the model, we adopted the following approach: Instead of performing the bootstrap on the whole original sample, we divided the original sample into as many groups as can be retained in the final classification. We then applied a single-component version of the HMTD model to each sub-sample separately in order to estimate the coefficients using bootstrap. By applying the model on the sub-samples separately, instead of on the initial sample, we avoided the so-called label-switching problem that is very common in latent variable clustering. The inconvenient of separate bootstrapping is that since we rely on the validated clustering solution, we ignore the model uncertainty including the weights of each cluster. We computed the confidence intervals using 1000 bootstrap samples, and we used the results to evaluate the significance of the estimated parameters.

All computations were done using R, and a specific package should be released soon. In the meantime, a first version of the R syntaxes is available on https://github. com/ztau/5352.

#### *2.3 GMM as a Gold Standard Alternative*

To evaluate the HMTD approach as a tool for clustering sequences of continuous data, we need a gold standard alternative. We choose the Growth Mixture Model (GMM) approach for that purpose, since it is the only true longitudinal clustering tool used in the social sciences.

Growth modeling includes several similar frameworks aiming to model and discover the patterns of individual changes in a longitudinal data framework (Reinecke and Seddig 2011; McArdle and Epstein 1987). The basic growth model assumes that all trajectories belong to the same population and that they may be approximated by a single average growth trajectory using a single set of parameters. However, several models extend these assumptions; for example, the latent class growth analysis (LCGA) model, which assumes null variance-covariance for the growth trajectory within each class (Nagin 1999; Jung and Wickrama 2008), and the heterogeneity model (Verbeke and Lesaffre 1996), which goes a bit further but still imposes the same variance-covariance structure within each group of subjects. Therefore, we discuss the more flexible GMM in this section and use it in our analysis as gold standard.

The GMM developed by Muthén and Shedden (1999), Bauer and Curran (2003), and Wang and Bodner (2007) is designed to discover and describe unknown groups of sequences that share a similar pattern. This method may be represented as a mixture of mixed-effects models in which each of the unknown subpopulations follows a distinct linear mixed-effects model. Its main advantage over other similar models—like the heterogeneity model (Verbeke and Lesaffre 1996)—is that it allows for estimation of a specific variance-covariance structure within each class (Francis and Liu 2015). Within-class inter-individual variation is possible for latent variables via distinct intercept and slope variances, represented by a class-specific fixed-effects and random-effects distribution. In other words, the variation in an expected group-specific trajectory is distinct for each group (heterogeneity in growth trajectories). Because of these advantages, the model is a reference point in continuous longitudinal data modeling with various applications in criminology (Francis and Liu 2015; Reinecke and Seddig 2011), health and medicine (Muthén and Shedden 1999; Ram and Grimm 2009), psychology, and social science (Muthén 2001), among others.

The GMM approach uses both observed and latent variables. The observed variables consist of a *p*-dimensional vector of continuous dependent variables *Y* (often a variable with repeated measurements) and a *q*-dimensional vector of covariates *X*. The latent variables are represented as a continuous *m*-dimensional vector *η*. Finally, to indicate the group in which each subject is included, we use a dummy variable with multinomial distribution stored in a k-dimensional binary vector *c* (Muthén and Shedden 1999). The equation of the GMM approach for individual *i* then becomes

$$Y\_l = A\eta\_l + \epsilon\_l,\tag{1}$$

where *Λ* is a *p* ×*m* parameter matrix (or matrix with basis vectors) that can be seen as a matrix of factor loadings, *ηi* is a vector of latent continuous variables, and *i* is an error term vector with zero mean.

In our case, the latent variable parameter matrix *Λ* has one column with parameters for the latent factor accounting for the intercept and another for the latent factor accounting for the slope. The general equation for every *η* is

$$
\eta\_l = Ac\_l + \Gamma \mathbf{x}\_l + \xi\_l,\tag{2}
$$

where *A* is a matrix with columns of intercept parameters for each class, *Γ* is an *m* × *q* parameter matrix and *ζi* is an *m* × 1 vector of zero mean residuals (and covariance matrix *Ψ*).

If we assume that some time-independent covariates *z* could influence the group membership *ci*, a multinomial logistic regression can be considered (with parameters *a* and *b*) as follows:

$$P(c\_l = K | z\_l) = \frac{\exp^{a\_k + b\_k z\_l}}{\sum\_{c=1}^{K} \exp^{a\_c + b\_c z\_l}}$$

An alternative notation of the model for subject *i* as part of class *k* at time *t* is

$$Y\_{l,t|c\_l=k} = X\_{1l}(t)^T \beta + X\_{2l}(t)^T \chi\_k + V\_l(t)^T u\_{lk} + w\_l(t) + \epsilon\_{l,t},\tag{3}$$

where *X*1*i* is a vector of covariates with common fixed effects *β*, *X*2*i* is a vector of covariates with class-specific fixed effects *γk*, and *Vi* is a set of covariates with individual class-specific random effects *uik*. Finally, *wi(t)* is an autocorrelated Gaussian process with null mean and covariance equal to *cov(wi(t)wi(s))* = *σ*2 *w* exp*(*−*ρ*|*t* − *s*|*)*.

The GMM is estimated by maximizing its likelihood using an ordinary EM algorithm. The continuous latent variables *η* and group membership variables *c* are considered missing data. The R package *lcmm* (Proust-Lima et al. 2017) was used to compute the GMM.

#### *2.4 Statistical Analyses*

To start with, we used the HMTD model to identify the best clustering of the IAT dataset without covariates, considering solutions from two to five groups. The best solution was selected on the basis of the Bayesian Information Criterion (BIC) (Raftery 1995). We then added covariates to this first model and analyzed the two resulting models, with and without covariates, particularly focusing on the IAT trajectories that did change group when covariates were added to the initial model. In order to isolate the impact of the covariates from any other computational issue or local optimum, we used the optimal solution obtained without covariates as a starting point for the full model. Therefore, we observe how this new model escapes the previous optimum.

We then computed the GMM models using the same dataset, and we compared the classifications obtained with the HMTD and GMM approaches. The usefulness of each covariate for discriminating between groups was evaluated using either a chi-square test for categorical covariates, or a single factor ANOVA for continuous ones. Notice that since it is not easy to compare two solutions with different number of clusters, we chose to compute a four-cluster GMM solution with all covariates instead of finding its own optimal number of clusters.

Our results are presented as figures displaying the IAT trajectories, and as tables describing the characteristics of subjects classified into groups and giving the HMTD model parameters.

#### **3 Results**

We provide here the results of the various clustering performed using the HMTD and GMM approaches, and we compare the resulting classifications. Notice however that given the iterative nature of the optimization algorithms, it is never possible to be sure that the final models are the best possible ones. Therefore, results should never be overinterpreted.

#### *3.1 HMTD Clustering*

Without covariates, the best model identified by the BIC is a four-component model (model 1). Figure 1 shows the IAT trajectories in each group. We clearly differentiate a group with average volatility and IAT level (group 1), a group with relatively low scores and variability (group 2), a group with very low variability and a low and constantly diminishing IAT score (group 3), and a group with more complex trajectories and hence variability (group 4).

**Fig. 1** IAT trajectories associated with each group in the four-group HMTD solution without covariates (model 1)

When we include the covariates in model 1 (Fig. 2) and relabel the four groups of the solution in order to match the groups of model 1, we obtain a similar fourgroup structure (model 2). As a comparison of the two figures might show, the most important difference is with the first two groups: group 2 of model 2 lost its highervalued trajectories and focused more on a low IAT-level and stable trajectories. This change will be explored in more details later.

Table 1 provides the parameter estimation for both models. In addition to the point estimates, we also provide the 95% bootstrap confidence intervals.

As regards the first model without covariates, the *θ*0 parameters giving the variance of each component of the model confirm the first impression given by Fig. 1: Group 4 is characterized by a much higher variability than the three other groups, and group 3 has the lowest variance, indicating less variation among the successive observations of a single individual. Parameters *ϕ*0 corresponding to the constant in the modeling of the mean of each component also take expected values, with higher values associated with groups showing higher average IAT level. Finally, the autoregressive parameter *ϕ*1 takes a value closer to one for the groups with trajectories showing smoother evolutions from one wave to the next, that is groups 1 and 3. All parameters of this first model are significant at the 95% level, as demonstrated by the confidence intervals.

**Fig. 2** IAT trajectories associated with each group in the four-group HMTD solution with five covariates (model 2)

As regards model 2, even if the first three parameters (*θ*0, *ϕ*0, and *ϕ*1) take values different from those of model 1, *θ*0 and *ϕ*1 take values in the same range as of model 1. On the other hand, important differences are found for the constant parameter *ϕ*0, and this parameter is no more significant in any group. Note that *θ*0 and *ϕ*1 tend to take smaller values in model 2. This can be interpreted as the first proof of interest of the covariates included in model 2: the groups are now more homogeneous (lower intra-group variance) and the explanation of a specific trajectory relies less on the immediately preceding observation. As regards the covariates, Age is never significant and could be eventually removed from the model. This could be due to the lack of a real age difference between participants (from 13 to 15 years old at baseline). However, the four other covariates remain significant for at least one of the groups.

When we consider each component of model 2 separately, the changes occurring in the trajectories associated with the first component are found related to the wellbeing of the concerned adolescents: a higher well-being such as measured by the WHO-5 index is significantly associated with a lower IAT-level. Males tend to have a lower IAT level than females, and a higher BMI is associated with higher IAT level. In group 3, a higher WHO-5 or BMI is associated with reduced IAT level, but being in the extended requirement school track is associated with a higher IAT level. Finally, in group 4, a higher WHO-5 or BMI is associated with reduced IAT level, and males tend to show a much higher IAT level than females.

Table 2 provides the main characteristics of the subjects classified into each group. For time-dependent variables, we considered the average value of each individual. A comparison is performed for each variable separately to test whether the groups are significantly different with regard to the variable. Considering only the two HMTD models, we observe that in addition to the expected differences in IAT level, the only other variable with significantly different values across groups is the WHO-5 measure of well-being. For both models, we observe two groups (2 and 3) with lower average IAT scores. The same two groups also display higher emotional well-being, as compared to the other groups, confirming previous results (Surís et al. 2014). No differences are observed for the other covariates, even if Gender comes close to significance in model 1. Even if not significant at the 95% level, probably because of the reduced sample size, we find a gender separation at the sample level; groups 2 and 4 contain a higher proportion of boys compared to the other two groups. The education track also shows a difference at the sample level: the first two groups contain more individuals following the highest education track as compared to groups 3 and 4. On the other hand, no notable difference is observed between the groups for Age and BMI, even if BMI, used as a covariate in model 2, is statistically significant in the modeling of the mean of each component.

**Table 2** Characteristics of subjects classified into groups for different clustering. The *p*-value gives the result of the test comparing the different groups for each variable. The number of sequences classified into each group is provided in brackets after the group number

#### *3.2 Usefulness of the Covariates*

From the results of the previous section, we find that the inclusion of covariates in the first classification obtained with the HMTD model helped us better differentiate the four groups, but without entirely changing their interpretation. We would like to better understand the changes in trajectory classification that occurred between these two models. Table 3 indicates how many subjects changed groups between the initial model without covariates and model 2 with covariates. As noted earlier, most of these changes occurred between groups 1 and 2. In particular, 19 second-group subjects of model 1 were transferred to the first group in model 2, and the steady

**Table 4** The characteristics of 19 subjects moving from group 2 to group 1 (group 2→1) as compared to subjects staying in the same group (either 1 or 2) in both HMTD classifications. The means (numerical variables) or proportions (categorical variables) are provided, and differences with the subjects remaining in the same group (either 1 or 2) are assessed using *t*-tests and *χ*2-tests with continuity correction: ns: non-significant

\**p <*0.05; \*\**p <* 0.01; \*\*\**p <* 0.001

low Internet addiction profile of the second group became even more pronounced, with the higher Internet addiction subjects joining the first group. However, since some trajectories simultaneously left group 1 for the three other groups, the average IAT level of group 1 also decreased. Overall, the inclusion of covariates appears beneficial for the differentiation of trajectory features among groups.

The 19 individuals who switched from group 2 to group 1 represent the main difference between the two models, with all the other changes concerning at the most seven subjects. Thus, it is interesting to explore how these individuals differed from those who remained in the first or second group in both classifications. Table 4 summarizes our findings using t-tests and *χ*2-tests to compare the different variables. The average IAT scores are quite different between the three considered sub-groups, and, as expected, the "moving" sub-group shows an Internet dependence level between the two "stable" sub-groups. Thus, the moving individuals were among the most Internet-problematic members of the full second group of model 1, and even if the average IAT score is not the only indicator of group affiliation, a visualization of the trajectories would confirm the ambiguous nature of these individuals. The moving subgroup is also significantly different from the group of individual staying in group 1 as regards the WHO-5 index of well-being and the gender ratio, with a higher well-being and higher proportion of males among the moving subgroup. No other significant differences are observed.

#### *3.3 GMM Clustering*

Without covariates, the best GMM solution in terms of BIC is a two-group solution, but given the high difference in number of trajectories associated to each group (169 vs 16), this solution in not really interpretable and hence less useful than the four-group solution given by the HMTD approach. Therefore, we also estimated a four-group GMM (Fig. 3).

In the two-group solution, a large majority of trajectories are associated with group 1, and only 16 sequences are associated with group 2. The average IAT level is higher in group 2, but both groups exhibit an important variability, as indicated in Table 2. Moreover, in terms of interpretation, one can only say that IAT sequences with a clear increasing trend are separated from the other sequences. In the four-group solution, even if the number of groups is the same as in the HMTD models, there is no a priori correspondence between the HMTD and GMM groups. In the four-group GMM solution, the number of subjects per group shows much more variability than that observed with the HMTD group, with the majority of individuals classified in groups 1 or 3, and only three subjects in group 4.

**Fig. 3** IAT sequences associated with each group in the four-group GMM solution without covariates

**Fig. 4** IAT sequences associated with each group in the four-group GMM solution with four covariates

Finally, as with the HMTD approach, we enhanced the four-group solution by adding covariates. Four of the five covariates used in the HMTD approach appeared useful in the GMM solution as well. Figure 4 displays the resulting groups obtained after adding Gender and Education as predictors for group membership (multinomial regression on *ci*), and WHO5 and BMI as fixed effect. On the other hand, Age was not included in the final model because the estimation process would then lead to a one-group solution. Another important issue with the GMM approach is the results' sensitivity to the order in which the covariates are included in the model. Various covariate combinations were tested before we chose the abovementioned combination as the best one in terms of clustering results. For instance, *classmb = gender + education sector* does not give the same results as *classmb = education sector + gender*. This surprising result may be due to a bug in the *lcmm* R package, but in our opinion the reason could rather be the optimization procedure. It is well known that EM-type algorithms converge to the nearest local optimum, and that this optimum is not always the global one. Therefore, the solution depends on the initial values of the parameters, especially when the solution space is complex, which is the case here.

As Fig. 4 shows, the number of trajectories associated with each group is quite variable, with the large majority assigned to group 1. The first two groups are characterized by low variability and an overall low IAT level. The trajectories in these two groups seem very similar, but since this four-group solution might be suboptimal and is computed only for the purpose of comparison with HMTD models 1 and 2, a three-group solution could merge these two groups into one group. The last two groups have a higher average IAT level, both exhibiting a general linear trend over time, decreasing in group 3 and increasing in group 4.

Table 2 gives the characteristics of individuals classified in each group of the GMM models and compares the groups for each variable. Note that given the large differences in group size, the test results for the GMM models should be interpreted with caution. As observed earlier in the HMTD case, significant differences exist between groups for both the IAT and WHO-5 variables. A significant difference exists also for BMI in the four-group GMM model without covariates. More interestingly, the Age and Education track at baseline also show significantly different values across groups, with one of the variables (Education track) being included in the model as covariable, but not the other. This difference between the HMTD and GMM clustering points to the fact that the solutions provided by both approaches are not identical or interchangeable, and that the two models used information in a different manner to provide usable data sequence clusterings.

#### **4 Comparison of HMTD and GMM**

When used for clustering purposes, the HMTD and GMM models share some characteristics: They both represent a kind of mixture model, they can include covariates of any type at the visible level, and they can also include covariates at the latent level and use them to estimate the initial probability of each cluster. However, HMTD and GMM also have several differences. First of all, since GMM is a mixture of mixed models, it is able to accept both fixed and random effects. Another difference is the possibility of HMTD to include an autoregressive specification for the variance and thus to allow for the clustering of longitudinal sequences whose variance evolves in time. For instance, sequences becoming more instable over time can more easily be grouped together. However, to exploit this feature, it is necessary to work with long data sequences, what was not the case here with the IAT example.

Another feature of HMTD that is worth stressing is the possibility of using it to perform different kind of clustering (Bolano and Berchtold 2016). The transition between components is driven by the hidden transition matrix *A*. In this paper, *A* was constrained to be a diagonal identity matrix, implying that each sequence was assigned to one and only one group, and all sequences assigned to the same group were described by the same visible model. However, there are several alternatives. For instance, different latent states may be required to alternate over time in order to find the optimal modeling of a given sequence. If *A* is constrained to have the following structure:

$$A = \begin{pmatrix} a\_1 \ 1 - a\_1 \ 0 & 0 \\ a\_2 \ 1 - a\_2 \ 0 & 0 \\ 0 & 0 & a\_3 \ 1 - a\_3 \\ 0 & 0 & a\_4 \ 1 - a\_4 \end{pmatrix}$$

where *a*1*, a*2*, a*3 and *a*4 are transition probabilities, then one performs at the same time a modeling and a clustering of the data sequences. The first two states are used to model the first cluster, and states 3 and 4 are used to model the second cluster. In other words, data sequences are clustered into two groups, but inside each group there are two different visible models allowing for a better representation of these sequences when their behavior evolves over time.

Another specification of *A* would allow some sequences to remain always in the same cluster, whereas other ones could transit at some point in time from the first to the second cluster:

$$A = \begin{pmatrix} a\_1 & 1 - a\_1 & 0 & 0 \\ a\_{21} & a\_{22} & 1 - a\_{21} - a\_{22} & 0 \\ 0 & 0 & a\_3 & 1 - a\_3 \\ 0 & 0 & a\_4 & 1 - a\_4 \end{pmatrix}$$

#### **5 Conclusion**

Hidden Markovian models are known to be valuable tools to analyze the dynamics in longitudinal continuous data and in life course data (e.g. Helske et al. 2018). The present study demonstrates that the sequences of continuous longitudinal data can also be classified into as many groups as required, and that the HMTD model can be used as a valid alternative to GMM. The inclusion of covariates has beneficial effects on clustering, because the resulting groups have lower intra-variability compared to the solution without covariates.

In a comparative study involving the use of GMM for clustering, our first finding is that the HMTD approach is a good alternative to GMM, because in terms of interpretability its results are at least as interesting as the results given by GMM. However, on the basis of just one practical example, we obviously cannot conclude that one approach is better than the other; moreover, this is not the purpose of this study. What we can conclude is that the HMTD approach is not only theoretically, but also practically useful to classify sequences of continuous data in mutually excluding groups.

In the literature, excessive Internet use has been found to be highly related to several somatic conditions, sleep disturbance in particular. However, in this paper, our main objective is not to explain IAT trajectories, but to find ways to classify such trajectories into meaningful groups. Moreover, there is still an ongoing debate on the direction of the relationship between Internet use and sleep disturbance, not to speak of causality. Therefore, we chose not to consider sleep disturbance in this analysis, but to concentrate on other covariates that are more neutral to IAT scores. Nevertheless, even with this restriction, the results obtained with the HMTD model are highly significant and allow for a sound interpretation. The four resulting groups differ in terms of average value and variability. The relationship observed between IAT and the emotional well-being of subjects suggests that both concepts are linked and that a higher risk of Internet addiction is related to poorer well-being. Gender is also a discriminating factor between groups, with a lower proportion of males in the first and third groups, but, given the small sample size, the differences are not significant at the population level.

The main strength of this study is the demonstration of the usefulness of the HMTD approach as a valuable alternative to the GMM approach for clustering continuous data sequences. Researchers would be advised to consider both approaches to take full advantage of the information in their data. However, some weaknesses of this study are to be mentioned. At the theoretical level, we include covariates in the HMTD model only at the visible level, but it is also possible to include them at the latent level as well in order to enhance the prior probabilities of each cluster. As regards the application of the model to IAT trajectories, we used a rather small convenience sample; this is not representative of the population of adolescents living in the canton of Vaud. More analyses need to be conducted with larger databases to define a real typology of IAT trajectories.

Overall, in spite of some shortcomings, the HMTD model can be considered as a complete framework for the analysis of continuous data sequences. It is an explanatory tool as well as a clustering tool, and by adding covariates, constraints on the transition matrix, and autoregressive modeling of the mean and variance of each component, the model goes well beyond the traditional Markovian models such as homogeneous Markov chains or hidden Markov models.

**Acknowledgements** This publication benefited from the support of the Swiss National Centre of Competence in Research LIVES – Overcoming vulnerability: Life course perspectives, which is financed by the Swiss National Science Foundation (grant number: 51NF40-160590). The<EMAIL_ADDRESS>study was financed by the Public Health Department of the Vaud canton and by the Swiss National Science Foundation (grant number: 105319\_140354). The authors are grateful to both institutions for their financial support. The funding bodies had no role in the design and conduct of the study; collection, analysis and interpretation of data; or preparation, review, or approval of the manuscript.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Divisive Property-Based and Fuzzy Clustering for Sequence Analysis**

**Matthias Studer**

#### **1 Introduction**

In this paper, we introduce property-based and *fuzzy* clustering and discuss their usefulness in the context of sequence analysis. We also present some tools available in R code by which to conduct these analyses. These two clustering methods aim to overcome some of the limitations of the more "traditional" ones, such as partitioning around medoids or agglomerative clustering.

Most of the clustering methods used in sequence analysis are polythetic, with the notable exception of model-based clustering. This means that cluster memberships are defined according to a broad set of properties and by comparing a given sequence to all the other ones. For this reason, the rules that define which sequences belong to which cluster are implicit. In other words, we do not know exactly on which grounds a sequence is assigned to a given cluster. Having a typology based on implicit rules of cluster membership has two disadvantages.

First, the resulting clustering is sample-dependent, which means it cannot be compared to another typology that is created in a subsample or another sample, for instance. Even comparing two typologies that appear to be similar might be problematic, since their underlying implicit clustering rules might differ.1 On the other hand, having explicit rules would allow one to validate a typology in other

M. Studer (-)

1We cannot be sure that the clustering rules are sufficiently similar, because we do not know them. They are only implicit.

NCCR LIVES and Geneva School of Social Sciences, University of Geneva, Geneva, Switzerland e-mail<EMAIL_ADDRESS>
samples and reproduce a given (validated) typology in other studies. Implicit rules therefore hinder the reproducibility of life-course research, as well as the possibility of undertaking a large-scale literature review.

Second, implicit rules make the interpretation of clustering more difficult. The first step after creating a typology is usually to try to interpret and recover these implicit rules by using different kinds of graphics and analyses. Explicit clustering rules would make the interpretation of clustering results much easier.

Monothetic clustering—which here we call "property-based clustering"—aims to create a sequence typology defined by explicit classification rules. In this paper, we introduce a method that is based on the "DIVCLUS-T" algorithm proposed by Chavent et al. (2007). Following the work of Piccarreta and Billari (2007), we also discuss its use in sequence analysis and extend their work by proposing for consideration new sets of state sequence features. Finally, we make the analytical results broadly available in the WeightedCluster package.

In this paper, we also discuss the use of *fuzzy* clustering, which aims to overcome another limitation of "traditional clustering." In sequence analysis, we usually use *crisp* clustering (i.e., each sequence is assigned to only one sequence type). In *fuzzy* clustering, however, each sequence belongs to one or more clusters, with a certain degree or strength (D'Urso 2016).

The *fuzzy* approach has several advantages over the more usual *crisp* one (D'Urso 2016). First, sometimes some sequences are between two sequence types. In *crisp* clustering, these sequences would be assigned to one of the two types; in *fuzzy* clustering, however, these sequences would be considered a *hybrid-type* or a *mixture* of the two types (D'Urso 2016). From a statistical viewpoint, *fuzzy* clustering might lead to better results if some sequences are between two (or more) sequence types. This case might occur frequently, according to Warren et al. (2015), who argue that exact cluster membership should not be trusted. From a sociological perspective, this approach is of special interest when the trajectories are not strongly structured into types, and when we can think that some individuals can be influenced by several sequence types.

Second, in *fuzzy* clustering, membership is thought to be *gradual*. Some sequences are more central (typical) of a given type than are others. This is also an interesting property from a sociological viewpoint. From the Weberian *ideal– typical* perspective, nobody is the perfect incarnation of an ideal type, but some are closer to it than others. This sociological perspective is similar to the gradual membership approach inherent in *fuzzy* clustering.

For all these reasons, the use of the *fuzzy* clustering approach is promising in life-course research and sequence analysis. However, to the best of our knowledge, it has been only seldom used in sequence analysis. One of the likely reasons for this paucity is the lack of proper tools by which to analyze sequences in conjunction with a membership matrix, instead of a categorical covariate (as in *crisp* clustering). For instance, Salem et al. (2016) used *fuzzy* clustering, but they ultimately assigned each sequence to the cluster with the highest membership in all subsequent analyses, thus turning in fact back to *crisp* clustering. In this paper, we propose different tools to fill this gap and make use of the full information of the membership matrix.

This paper is organized as follows. It starts by presenting property-based clustering and the set of sequence properties whose consideration we propose. We then turn to the presentation of *fuzzy* clustering and introduce several ways of representing the results, before interpreting them. We also discuss how to properly analyze cluster membership according to some explanatory covariates. Before concluding, we briefly show how to run the proposed analysis in R, using our WeightedCluster library.

#### **2 Sample Issue**

The usefulness of the proposed methodology is illustrated by using the data and the research question from McVicar and Anyadike-Danes (2002), who studied schoolto-work transition in Northern Ireland. Their analysis was undertaken in two steps. They started by identifying ideal-typical trajectories, before explaining clustering membership by using information such as qualification at the end of compulsory schooling, family background, and demographic characteristics. Their aim was to "identify the 'at-risk' young people at age 16 years and to characterize their postschool career trajectories" (p. 317). To build our clustering, we use optimal matching with constant cost.

#### **3 Property-Based Clustering**

The aim in using property-based clustering is to build a sequence typology by identifying well-defined clustering rules that are based on the most relevant properties of the analyzed object. In the literature, these clustering methods are called *monothetic divisive* clustering methods (Chavent et al. 2007), and they were first introduced in sequence analysis by Piccarreta and Billari (2007). We propose here a conceptual presentation of the "DIVCLUS-T" algorithm (a detailed presentation can be found in Chavent et al. 2007). The "DIVCLUS-T" algorithm is very similar to one proposed by Piccarreta and Billari (2007). Our choice between the two is based on availability in R. Here, we mainly extend the work of Piccarreta and Billari (2007), by proposing new sets of sequence features for consideration.

#### *3.1 Principle*

Property-based clustering uses two kinds of information. First, the sequence properties are used to build the rules. Second, it uses a dissimilarity matrix, which is used to measure variation in the sequence and how much of this variation can be explained by a given property; this is in line with the discrepancy analysis framework (Studer et al. 2011).

The method then works in two steps: tree building and splits ordering. In the tree building phase, all sequences are grouped in an initial node. Then, this node is split according to one of the object properties, into two subnodes or clusters. This property and the associated split are chosen in such a way that the split "explains" the biggest share of the sequence discrepancy (Studer et al. 2011; Chavent et al. 2007; Piccarreta and Billari 2007). The process is then repeated on each new node until a stopping criterion is found. First, the algorithm might stop because there are no further relevant properties by which to make a split. Second, nodes with only one observation are obviously not split. This first step is roughly equivalent to the procedure that Piccarreta and Billari (2007) propose.

As in Piccarreta and Billari (2007), our implementation of the "DIVCLUS-T" algorithm also makes it possible to specify a minimum number of observations per node. Splits that would lead to at least one node with fewer than this minimum number of observations are discarded. This is a useful extension, if we want to ensure that all clusters represent at least a given percentage (for instance, 5%) of the sequences. One might also restrict to "significant" splits by using permutation tests, as in discrepancy analysis (Studer et al. 2011). However, the usefulness of the latter approach is subject to discussion, as the concept of significance is not very well defined in cluster analysis. This first step of the procedure can be seen as a decision tree, where the splits are chosen according to the explanatory power of the considered properties. In fact, our implementation is based on the tree-structured discrepancy analysis.

Once the whole tree is built, the splits are ordered according to their overall "relevance." More precisely, this "relevance" is measured by calculating the increase in the share of the overall discrepancy that is explained by adding a split. This procedure has the advantage of maximizing a global criterion. Ultimately, the result of this procedure is a series of nested partitions ranging from one group to a number of groups, any of which depend on the stopping criteria or a maximal number of groups to consider.2 This second step is the major difference from the procedure proposed by Piccarreta and Billari (2007), where the final clustering and the stopping criteria depend on a pruning procedure.

Figure 1 graphically represents the procedure, using our illustrative example of school-to-work transition in Northern Ireland, for the first nine splits. The order of the splits is presented on the right, with the associated number of clusters. We start at the top of the tree with a single cluster. At this stage, the most relevant feature in splitting this top node into two is to have spent more (or less) than 17 months in higher education (property "duration.HE"). Stopping at this level would result in a cluster in two groups (presented on the right of Fig. 1). The clustering in three groups is obtained by splitting the node "less than or equal to 17 months in higher education" in two groups: having spent more (or less) than 33 months in employment (criterion "duration.EM"). This last split was preferred to the solution of splitting the node "more than 17 months in higher education," because it has greater explanatory power regarding sequence variation at a global level. The procedure then continues until some stopping criteria are met.

2In other words, new groups are added by dividing one of the existing groups into two subgroups.

**Fig. 1** Clustering tree

As with any agglomerative or divisive clustering, the choice of the number of groups can be grounded on cluster quality measures (see Studer 2013, for a review). Here, various clustering quality measures agree on choosing either the solutions in four (*ASW* = 0*.*41) or nine (*ASW* = 0*.*38) groups as the best clustering. According to the empirical evaluation of Chavent and Lechevallier (2006), for a small number of groups, the "DIVCLUS-T" algorithm tends to produce better clustering (measured on the basis of statistical criteria) than Ward clustering. However, the reverse becomes true as the number of groups increases.

#### *3.2 Property Extraction*

The results of property-based clustering are highly dependent on the properties of the object to be clustered. In the empirical evaluation made by Chavent and Lechevallier (2006), having a large number of meaningful properties was one of the key elements that led to good-quality clustering. In this section, we propose a set of properties worthy of consideration. In our implementation of the algorithm, these properties are automatically extracted.

Within the life course paradigm, three main dimensions of the trajectories are of central interest: the timing, the duration, and the sequencing of the states (Studer and Ritschard 2016). We propose the automatic extraction of various properties that measure each of these dimensions.

To measure the timing of the state, we consider the state at each time position *t*. If we consider the sequence of length , we therefore end with categorical covariates that measure the situation over time. Piccarreta and Billari (2007) propose another way of measuring the timing, by considering the spells that form the sequences. They generate one property *As,k* that stores the age at the beginning of the *k*th spell in state *s*. If the property is not observed (i.e., there is no *k*th spell in state *s*), they propose setting it to + 1, where is the length of the sequence. However, we take here a slightly different strategy: we set it to a missing value. Missing values are then treated as a special case when defining a split. Although we use the numeric information, when available, to define the split by using an inequality, the missing values are attributed to one of the two groups—whichever gives the best result.

We propose the use of two sets of properties to measure duration. First, we consider the duration in the successive spells. This is achieved by generating one property *Ds,k* that stores the duration of the *k*th spell in state *s*. Second, we consider the overall time spent in each state that results in one property per state.

We use frequent subsequence mining to extract the properties that measure the state sequencing. With this method, the aim is to uncover frequent subsequences within a set of sequences (Studer et al. 2010; Agrawal and Srikant 1995; Zaki 2001). A subsequence *s* is defined as a subsequence of *x* if all the states of *s* occur in the same order as in *x*. For instance, the sequence *A* − *C* is a subsequence of sequence *A* − *B* − *C* because *A* and *C* occur in the same order. A subsequence is said to be frequent if it is found in more than a predefined percentage of sequences. Using this framework, several sets of sequence properties can be extracted. First, we look for frequent (in our case, 1%) subsequences in the sequence of distinct successive states (DSS). This step generates one property (i.e., variable) per identified frequent subsequence, and stores the number of times that the subsequence is found in each sequence. For instance, the subsequence *A*−*C* occurs twice in the sequence *A*−*B*− *C* −*B* −*C*. We also consider the age at the first occurrence of the pattern (i.e., when the pattern starts). Second, we look for frequent subsequences within the transition sequences. This is achieved by specifying a distinct state for each transition. For instance, the DSS *A* − *B* − *C* will be recoded as *A* − *AB* − *BC* before running the analysis. Here again, the number of occurrences and the age at the first occurrence are stored as properties.

Finally, the user can consider and add other properties. The algorithm computes different sequence complexity measures. Piccarreta and Billari (2007) suggest adding information about covariates such as education level. Application-specific sequence properties could also be of interest. For instance, one might have an interest in adding the time spent in a state of joblessness within the last 12 months of our sequences, if professional integration at the end of the sequence is of primary concern.

Although all these properties are automatically extracted, they need to be carefully chosen according to the issue under investigation. For instance, for a study that mainly concerns itself with timing differences, we suggest restricting attention to timing-related properties. For this reason, in our implementation of the algorithm, one can specify the sets of properties to be considered. Table 1 summarizes the properties considered in this study. The first column contains the name of the property used in our R implementation, and the second provides brief descriptions of the sets of properties.

**Table 1** Sequence properties considered in the clustering algorithms

#### *3.3 Running the Analysis in R*

Property-based clustering can be run using the seqpropclust function available in the WeightedCluster package. Aside from the state sequence object myseq, one needs to specify the distance matrix (argument diss), the properties (by default, all properties are computed), and the maximum number of clusters under consideration.

```
R> ## Clustering using properties "state" and "duration"
R> pclust <- seqpropclust(myseq, diss=diss, maxcluster=5,
       properties=c("state", "duration"))
R> ## Displaying the resulting clustering
R> seqtreedisplay(pclust, type="d", border=NA, showdepth=TRUE)
R> ## Computing clustering quality and cluster membership
R> pclustqual <- as.clustrange(pclust, diss=diss, ncluster=5)
```
The clustering membership can be extracted by using the as.clustrange function, which also computes various cluster quality measures. See Studer (2013) for a detailed presentation of this procedure.

#### **4 Fuzzy Clustering**

In *crisp* clustering, each sequence is assigned to exactly one sequence type. The result is a categorical covariate that summarizes the typology. In *fuzzy* clustering, each sequence can belong to more than one cluster; this is achieved by computing the degree or strength of membership of each sequence to each identified sequence type (D'Urso 2016). This is of central interest when the sequences are not thought to be strongly structured, or when some sequences could have been influenced by more than one type.

More precisely, the result of *fuzzy* clustering is a membership matrix **u** comprising one row per individual and one column per cluster. Each value *uiv* of this matrix measures the membership strength of an individual *i* to each cluster *v*. These membership degrees, which usually sum to 1, are also called "probabilities," and they lead to two slightly different interpretations. The concept of membership "strength" or "degree" refers to the closeness of each sequence to each type. The notion of "membership probabilities" refers to the chances that the underlying sequence was generated according to one of the types.

In this section, we start by presenting the algorithm used herein. We then propose different approaches to describing and visualizing the results of *fuzzy* clustering, using the membership matrix. Finally, we discuss possible strategies by which to analyze how explanatory covariates are linked to cluster membership; again, this is based on the membership matrix. We hope that the availability of these tools will lead to more widespread use of *fuzzy* clustering in sequence analysis.

#### *4.1 Fanny Algorithm*

We use the *Fanny* algorithm proposed by Kaufman and Rousseeuw (1990) and later adapted by Maechler et al. (2005). This algorithm aims to minimize the following function:

$$\sum\_{v=1}^{k} \frac{\sum\_{i=1}^{n} \sum\_{j=1}^{n} u\_{lv}^r u\_{jv}^r d(i, j)}{2 \sum\_{j=1}^{n} u\_{jv}^r} \quad , \tag{1}$$

where *n* is the number of observations, *k* a predefined number of groups, *uiv* the membership value of individual *i* to cluster *v*, and *d(i, j )* the distance between sequences *i* and *j* . The exponent *r* is a *fuzziness* parameter that needs to be set by the user. A value of 2 is often used, but values between 1.5 and 2.5 are usually recommended (D'Urso 2016). The standard procedure is to start with *r* = 2 and use a smaller value if the algorithm does not converge.

Using our sample data, we used a value of 1*.*5 for the fuzziness parameter. We kept a solution in seven groups, based on the interpretability of the resulting clustering and the aim of the study.

#### *4.2 Plotting and Describing a Fuzzy Typology*

Once the clustering has been computed, typically, the first step is to describe each cluster and give it a first interpretation. In this section, we propose several approaches to doing so by using the membership matrix.

#### **4.2.1 Most Typical Members**

A first way to label the clusters and interpret them is to identify typical sequences. In "traditional sequence analysis clustering," this can be done by identifying the medoid or a representative sequence based on other criteria, such as neighborhood density (Gabadinho et al. 2011). A natural way to do it with *fuzzy* clustering is to choose the sequence with the highest membership strength in each cluster. The first row of Table 2 presents this information. Using this strategy, we can have a first look at our clustering membership matrix. We found a first cluster related to full employment, then three patterns of education (i.e., training or further education) followed by employment: two patterns leading to higher education, and one pattern of training followed by joblessness.

Table 2 also presents descriptive statistics of membership strength for each cluster. As one will recall, in *fuzzy* clustering, each sequence has a measure of its membership strength in each cluster. Hence, it is not possible to compute a percentage of sequences belonging to each cluster, as we would do for *crisp*

*s*1 *u*13 3

clustering. However, average cluster membership provides similar information: it can be interpreted as the relative frequency of each cluster, if sequences are weighted according to their membership strength.

The maximal value is also interesting, as it provides an estimation of the quality of the chosen representative. The higher the membership, the better the representative (i.e., a value of 1 would identify a sequence that fully belongs to that cluster). In some clusters, the maximum is quite low, if we consider that the maximal possible value is 1. For instance, in the training–joblessness cluster, the maximum equals 0.8. Hence, our representative is also close to other clusters—perhaps cluster 2. To describe the clusters, we therefore need to take into account more information than just the sequences with the highest membership.

#### **4.2.2 Weight-Based Presentation**

Our second proposition in analyzing the *fuzzy* cluster is to weigh the sequences according to their membership strength or probabilities. We augment the dataset by repeating the sequence *si* of individual *i k* times (i.e., once per cluster). We therefore have *k* sequences for individual *i*, denoted as *si*1 ···*sik*. We weight these sequences according to their membership degree *ui*1 ··· *uik*. Hence, even if the same sequence were repeated *k* times, its weights will sum to 1. We then create a new categorical covariate in this augmented dataset, and it specifies the cluster (ranging from 1 to *k*) of the associated membership degree.

Table 3 presents a small example, to make the presentation clearer. Suppose we have three clusters named 1*,* 2*,* and 3. For individual 1 with the sequence *s*1, we have three observations (i.e., one per cluster). The first observation is weighted according to the strength of membership to the first cluster, and the associated cluster covariate is set to 1. We then repeat the same procedure for each observation.

**Fig. 2** Membership-weighted plots of the sequence for clusters "(SC,25)-(HE,45)" and "(TR,22)- (JL,48)"

This weighting strategy allows us to use any tools available for weighted sequence data. For instance, we can use a sequence distribution plot. Figure 2 proposes several plots of these membership-weighted sequence data for the clusters "School–Higher Education" and "Training–Joblessness." Subfigures (a) and (e) present sequence distribution plots for these clusters. If the cluster "School–Higher Education" seems to be quite well defined, the "Training–Joblessness" one shows more discrepancy, as we already noted.

This weighting strategy is also supported from a more statistical perspective. Minimizing Eq. 1 is equivalent to minimizing a residual sum of squares in a discrepancy analysis of this augmented dataset (Studer et al. 2011). More precisely, it minimizes the residual sums of squares of this augmented dataset, where each sequence is weighted *ur i*1 ··· *ur ik* and the explanatory categorical covariate would be the cluster 1 *...k*.

Following this reasoning, we can weigh the sequences by using the exponent (here, *r* = 1*.*5). The result might be closer to the underlying algorithm. However, the interpretation is more difficult and, as we will see, it is also interesting in the following analysis to rely on the membership strength. For this reason, this approach should be used mostly to describe the clusters. The result of this strategy is shown using a d-plot in subfigures (b) and (f); in both cases, the cluster seems to be better defined.

By using index plots, we can take a closer look at the longitudinal patterns. In this case, we additionally suggest ordering the sequences according to membership degree. The result is shown in subfigures (c) and (g). The most typical sequence lies at the top of the subfigures, with a high membership degree; meanwhile, the bottom shows less-characteristic patterns. Interpretation should be made with caution, as it depends on the maximal membership degree. In the cluster "School–Higher Education," this maximum is close to 1, while in the other one it reaches only 0*.*8.

It can be interesting to focus on the sequences with the highest membership. In subfigures (d) and (h), we discarded sequences with a membership degree below 0.4. Interestingly, among the sequences with the highest membership in the "Training–Joblessness" cluster, we find sequences starting in "Further Education" or "Employment," for instance.

We propose several methods by which to visualize and describe a *fuzzy* typology; these methods allow us to properly interpret this typology. However, most sequence analysis applications go beyond the typology description by studying the factors that influence the kinds of trajectories being followed. We now turn to this kind of analysis for *fuzzy* clustering.

#### *4.3 Analyzing Cluster Membership Using Dirichlet Regression*

In typical sequence analysis, one often relies on multinomial regression to explain cluster membership (Abbott and Tsay 2000). The aim is to identify how covariates explain the trajectory type that is followed. This cannot be done with *fuzzy* clustering, because our typology is described by a membership matrix and not by a categorical variable. Assigning each sequence to the cluster with the highest membership strength is not a solution either, for in doing so, we would lose all the added value inherent in *fuzzy* clustering.

Several models are available to analyze membership matrices, which can be seen as "compositional data" (Morais et al. 2016; Pawlowsky-Glahn and Buccianti 2011). Here, for two reasons, we suggest relying on Dirichlet regressions (Maier 2014), which are extensions of beta regression (Ferrari and Cribari-Neto 2004). First, interpretations of them are very similar to those of multinomial models, if we use the so-called alternative parametrization. Second, good performance is reported in Maier (2014) and in Morais et al. (2016), even under some violations of the statistical assumptions. In the current model, one of the categories (i.e., clusters) is chosen as the reference; we then estimate how explanatory factors influence the likelihood of being fully classified in a category, rather than in the reference.

Interpretations of the coefficients are very similar, then, to the multinomial ones, and they can be interpreted in the usual log-odds scale. Their exponents can therefore also be interpreted as "odds-ratio" values on cluster membership. In a Dirichlet regression, one can also estimate the effect of covariates on a "precision" parameter that measures the precision of estimation. (This parameter is named "precision" because it takes a high value when the residual variance of the dependent variable tends to be lower.) This can be used to take into account possible heteroscedasticity.

Table 4 presents the coefficients of the Dirichlet regression. We used the employment cluster as the reference. To simplify our presentation, we included only three covariates: gcseq5eq (the qualifications gained by the end of compulsory

**Table 4** Dirichlet regression of cluster membership

Log likelihood = 6402*.*14; Num. obs. = 712; Precision= 1*.*22∗∗∗(0*.*02); ∗∗∗ *p <* 0*.*001, ∗∗ *p <* 0*.*01, ∗ *p <* 0*.*05

education: five or more GCSEs at Grades A to C, or equivalent), Grammar (grammar school secondary education), and funemp (father unemployed at the time of the survey).

Let us interpret the coefficient of our covariate on the membership degree (adequacy with) or chance (if we think about probabilities) to follow the more "at-risk" pattern "(TR,22)-(JL,48)," as opposed to employment. We observe no significant difference between those having had a grammar school education and those who had not. However, individuals with unemployed fathers did tend to have a higher membership in this cluster (significant positive coefficient) than in the employment cluster—or, if we take the "probability" interpretation, they have a lower chance of following this pattern than the reference (employment). The same applies to those who had the five-grade qualification (variable gcse5eq) probably because they are very unlikely to follow the reference employment trajectory.

Additionally, it is often useful to understand the distinctive features of each cluster. For *crisp* clustering, this can be achieved by running a logistic regression on a dummy variable that measures cluster membership. Here, we can make use of beta regression, which aims to model a dependent variable that lies in the [0*,* 1] interval (Ferrari and Cribari-Neto 2004).3 The interpretation of the coefficient is similar to that of the Dirichlet regression. The exponent of the coefficients can be interpreted as an "odds-ratio" on cluster membership. Here, again, a "precision" parameter allows us to take into account over- or under-dispersion. The results lead to similar conclusions but further highlight that those who had the five-grade qualification (variable gcse5eq) are very unlikely to follow the employment trajectory type of sequence.

3Unlike logistic (binomial) regression, beta regression does not assume that the dependent variable is a proportion (i.e., the result of a count of 0 and 1). Furthermore, it can cope with over- or underdispersion.

#### *4.4 Running the Analysis in R*

The Fanny algorithm is available in the cluster package, through the fanny function. Aside from the distance matrix diss, one needs to specify the number of groups (argument k=7) and set the argument diss=TRUE to specify that we provided a distance matrix and not a dataset. Finally, the value of the fuzziness parameter *r* can be set through the memb.exp argument (default value of 2). The returned object provides the membership matrix (fclust\$membership) and additional information such as quality measures or related *crisp* clustering.

```
R> ## Fuzzy clustering in 7 groups using r=1.5
R> fclust <- fanny(diss, k=7, diss=TRUE, memb.exp=1.5)
R> ## Display the resulting clustering with membership
     threshold of 0.4
R> fuzzyseqplot(myseq, group=fclust$membership, type="I",
         membership.threshold=0.4, sortv="membership")
R> ##Estimation of Dirichlet Regression
R> ##Dependent variable formatting
R> fmember <- DR_data(fclust$membership)
R> ## Estimation
R> bdirig <- DirichReg(fmember~var1+var2|1,
         data=mydata, model="alternative")
R> ## Displaying results of Dirichlet regression
R> summary(bdirig)
R> ## Estimation of beta regression
R> breg1 <- betareg(fclust$membership[, 1]~var1+var2, data=mydata)
R> ## Displaying results
R> summary(breg1)
```
All the visualizations proposed here are available in the WeightedCluster package, through the fuzzyseqplot function. The function works in the same ways as the usual seqplot function available in TraMineR, except that the group argument should be a membership matrix or a fanny object. Furthermore, one can specify a membership threshold (for instance, 0.4) and whether graphics should be weighted by membership strength. If one wants to weights the sequences using the fuzziness parameter, one should set memb.exp to the correct value. By default, the fuzziness parameter is not used; hence, the memb.exp=1. When using index plots (type="I"), one can additionally set sortv="membership" to sort the sequences in each plot according to their membership strength.

Dirichlet regression can be estimated using the DirichReg function in the DirichletReg package (Maier 2014), while the beta regression can be computed with the function betareg available in the betareg package (Cribari-Neto and Zeileis 2010). For the former, the dependent variable should first be formatted using the function DR\_data before estimating the model using DirichReg. For beta regressions, a separate regression should be estimated for each cluster. One needs to specify the cluster membership strength on the right-hand side of the R formula, while adding covariates on the left-hand side as usual. In both cases, one can set a data frame where covariates should be found, using the usual data argument.

#### **5 Conclusion**

In this paper, we introduced two alternative clustering methods, each of which has its own strengths. We believe that property-based clustering is a very promising sequence analysis tool. Having clustering membership rules allows one to reproduce and validate a typology; furthermore, it significantly simplifies the interpretation of the clustering.

Property-based clustering is also useful in understanding the underlying criteria used by a dissimilarity measure to compare trajectories. For instance, in our example application, all splits were made according to the overall time spent in different states. This prevalence of duration illustrates once again that optimal matching tends to favor duration while comparing sequences. The use of other distance measures such as those reviewed in Studer and Ritschard (2016) or those introduced in this bundle in Collas (2018) or Bison and Scalcon (2018) would lead to the selection of other properties. For instance, sequence pattern properties would probably be selected by using a distance sensitive to sequencing, such as SVRspell (Elzinga and Studer 2015).

On the other hand, *fuzzy* clustering has been seldom used in sequence analysis. Nonetheless, the method should be useful in many situations. First, in many cases, exact cluster membership is doubtful (Warren et al. 2015). *Fuzzy* clustering allows one to relax the assumption that cluster memberships have been correctly retrieved by the cluster analysis; it does so by allowing multiple cluster memberships. This is also an interesting perspective from a sociological viewpoint, as trajectories might be influenced by several trajectory types. Second, in *fuzzy* clustering, membership is thought to be gradual; this too is interesting from a social science perspective. Some trajectories might be more typical of a type than others.

The aim of this study was to develop tools by which to facilitate the use, interpretation, and analysis of both clustering methods. However, further application of these methods is still needed to fully assess their strengths and weaknesses with regards to sequence analysis. We believe that this study is a first step in that direction.

**Acknowledgements** The author warmly thanks the two reviewers for their insightful comments on an earlier version of the manuscript. This publication benefited from the support of the Swiss National Centre of Competence in Research LIVES Overcoming vulnerability: Life course perspectives (NCCR LIVES), which is financed by the Swiss National Science Foundation (grant number: 51NF40–160590). The author is grateful to the Swiss National Science Foundation for its financial assistance.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **From 07.00 to 22.00: A Dual-Earner Couple's Typical Day in Italy**

### **Old Questions and New Evidence from Social Sequence Analysis**

**Ivano Bison and Alessandro Scalcon**

#### **1 Introduction**

How do dual-earner couples organize their workdays and how do they (de)synchronize their daily activities? These are the questions that we address in this paper using a multichannel sequence analysis approach. Our purpose is to consider the couples' division of work-family activities in holistic terms by setting it within the context of everyday life, that is, the overall temporal pattern of combination of His and Her multiple activities.1

Our multichannel approach is based on a Lexicographic Index (Bison 2011) that seeks to overcome some optimal matching limits of sequence analysis (Bison 2009). The case-study concerns Italian dual-earner couples and uses data from the Italian Time Use Survey 2008 (Istat 2011).

We know that for dual-earner couples the risk of experiencing a certain "lack of family time" is higher than for other couples (Saraceno 2012), due to the combination and the rigidity of His and Her work constraints. The spouses of these couples face various challenges: according to their working schedules, they are required to find the right amount of time for their family—i.e. housework, childcare and other non-paid work—as well as with their family—i.e. desirable

1We immediately point out that, in this paper, we only consider heterosexual couples, because of the limitations of the Italian Time Use Survey questionnaire. Moreover, we stress that the choice of using the male pronoun before the female one is perfectly conscious: to make the reading easier, we needed to follow a single criterion and we decided to cite the spouses following the order of records in our data files, i.e. male-female.

I. Bison (-) · A. Scalcon

Department of Sociology and Social Research, University of Trento, Trento, Italy e-mail<EMAIL_ADDRESS>
and shared activities, such as free time—integrating collective needs with individual ones. In other words, the time scarcity of dual-earner couples obliges them to adopt a *complementary strategy* (Mansour and McKinnish 2014) in order to reconcile their multiple family needs of "production" and "consumption", while preserving their personal satisfaction with work and daily life. The two individual careers have to coexist with a third one: family life.

In this scenario, an important component of a dual-earner couple's strategy is synchronization/desynchronization. For instance, previous studies show that a certain degree of desynchronization of working schedules may be a useful solution for partners because it can promote a more equal division of housework and child care (Presser 1994; Chenu and Robinson 2002; Lesnard 2008; Naldini and Saraceno 2011). At the same time, it is recognized that a certain degree of synchronization in work commitments can encourage the partners to spend time together in other desirable activities (Hamermesh 2002; Lesnard 2008).

More generally, looking at (de)synchronizations is important because they reveal a latent behavioral pattern of different work-family specializations and suggest new explanations for the continuing persistence of gender inequalities in the division of work-family activities. The study of (de)synchronizations could enable identification of multiple equilibria (Esping-Andersen et al. 2013). According to Esping-Andersen and colleagues, such study could reveal different behavioral patterns of work-family specializations—(i) egalitarian, (ii) unstable and (iii) traditional resulting from systematic co-action by different generative mechanisms, both symbolic-cultural (e.g. Berk 1985; West and Zimmerman 1987) and economicmaterial (e.g. Becker 1964; Coverman 1985; Manser and Brown 1980).

How do we measure (de)synchronizations? There are two radically different main approaches: time budgets—the dominant approach in the time use literature and sequence analysis.

In the former approach, we measure (de)synchronizations as the amount of time in which both the spouses have or have not done an activity at the same time or have or have not spent that time together in the same place. Thus, we obtain synchronicity ratios or percentages. However, in this way nothing is known about when the activity schedules overlap. This is a crucial limitation for two main reasons.

First, time is socially structured, and so too are social rhythms and constraints. Hence, being simultaneously at work at 10 a.m. or 10 p.m. has radically different impacts on a couple's daily life. Furthermore, different timings in working schedules may have radically different impacts on daily life if combined with other time demanding features, like for instance the institutional constraints of children's schedules (e.g. school hours).

Second, take the case of a full-time shift perfectly synchronized with a part-time afternoon shift: by considering only the duration of the overlap, we will mistakenly classify it as a highly desynchronized working schedule. However, such kinds of structural desynchronization—due simply to differences in duration—should not be compared with hypothetical other kinds with the same off-scheduling amount but different and more complex organization during the day (Nock and Kingston 1984).

According to Lesnard (2008), if we know little about how family time is daily balanced with work time for both spouses, this is mainly because we are used to adopting the too simplistic approach of the dominant time-budget perspective. To date, scholars have underestimated the importance of daily scheduling, while paying more attention to total amounts of time (Lesnard 2008). They have traditionally acquired time budget information related to different daily activities, but these should be seen in a holistic perspective that makes it possible to study the couple's days as a whole, avoiding the manipulation of time as if it were clay.

An alternative to the time-budget approach is sequence analysis (Lesnard 2008). According to Hallberg (2003), "while the traditional time allocation model typically studies the total time spent in, e.g., market work, over a day or a week, it provides little or no insight into the temporal pattern of time-use and therefore, potentially, misses a vital part of the mechanisms underlying empirical observations". A sequence analysis of time-use would evidence the routine aspects of daily life, as well as the couples' projects (Hagerstrand 1982), the performance of their complementary strategy across several daily constraints and unexpected events (Hellgren 2014). Finally, the analysis of time-use temporal patterns—instead of time budgets—seems more relevant in the study of the daily strategies and behaviors of a couple (Hallberg 2003).

We have pointed out that in order to understand the complexity of work-family balance strategies, it is necessary to study the couple's daily time-use pattern as a whole and in a more holistic way by adopting a multichannel sequence analysis approach. In the following section, we introduce the Lexicographic Index used to measure the resemblance between multinomial sequences. Section 3 sets out the data and methods. Section 4 presents the main statistical and graphical results of this study. Finally, Sect. 5 is devoted to summing up the main findings.

#### **2 The Lexicographic Index**

There are three main problems with current techniques used to compute the distances among sequences. One derives from the way in which similarity between two sequences is defined (Abbott and Tsay 2000; Wu 2000; Dijkstra and Taris 1995; Elzinga 2003; Bison 2009); the second is how to handle multinomial sequences (Abbott 1990); the third is how to treat a multichannel sequence as a whole (Gauthier et al. 2010).

Here we present an alternative method for computing distances among sequences. The lexicographic index (Bison 2011) is based on the sorting order of two different modes of observing events in a binary sequence. The first order is given by duration, that is the quantity of time, and is therefore based on the total number of observed events *u* in the sequence *x*. The second order is timing, that is when this event happens, i.e. the 'places' *sk* (*k* = 1*,...,u*) in the sequence when 1 occurs.2 For instance, we may have only three binary sequences of length 3 and *u* = 1. They differ according to when the event occurred; at time *t*1, *t*2 or *t*3. Hence, we may order these sequences [100], [010] and [001] according to the time order of events. Because the nature of the sorting order is double, the proposed index consists of two distinct parts.

The first part, *d (x)*, ranging from 0 to 1, takes account of the duration and therefore the different amounts of realization *u* recorded in the sequence:

$$d'(\mathbf{x}) = \boldsymbol{\mu}/T \text{ for } \boldsymbol{\mu} > 0 \qquad \text{and} \qquad 0 \text{ for } \boldsymbol{\mu} = 0 \tag{1}$$

where *T* is the length of the sequence.

The second part, *d(x)*, ranging from 0 to 1, takes account of timing and therefore the different numbers of combinations displayed by the sequences with variation in the amount of time. It is

$$d''(\mathbf{x}) = \begin{cases} 0 & \text{for } \boldsymbol{\mu} = \mathbf{0} \\ \frac{\binom{T}{\boldsymbol{u}}}{\binom{T}{\boldsymbol{u}} + 1} \frac{1 + \binom{T}{\boldsymbol{u}} - \left[ \binom{B\_{\boldsymbol{u}}}{C\_{\boldsymbol{u}}} - \sum\_{k=1}^{\boldsymbol{u}} \binom{B\_{\boldsymbol{k}}}{C\_{\boldsymbol{k}}} - \binom{A\_{\boldsymbol{k}}}{C\_{\boldsymbol{k}}} \right]}{\binom{T}{\boldsymbol{u}}} & \text{for } \boldsymbol{0} < \boldsymbol{u} < T \\ 1 & \text{for } \boldsymbol{u} = T \end{cases} (2)$$

where *Ak* is the exact position of *sk* in the sequence, *Bk* is the last position that *sk* can occupy within the sequence, and *Ck* is the first position that *sk* can occupy. For example, for sequence [0101], with *T* = 4 and *u* = 2, we have: for *s*1 the exact place of the first 1 is *A*1 = 2, the last position is *B*1 = 3 and the first position is *C*1 = 1; for *s*2 the exact place of the second 1 is *A*2 = 4, the last position is *B*2 = 4 and the first position is *C*2 = 2, the value of *d(*[0101]*)* is 0.285714 and is obtained as

$$\frac{\frac{4!}{2!(4-2)!}}{\frac{4!}{2!(4-2)!}+1} \frac{1+\frac{4!}{2!(4-2)!} - \left(\frac{4!}{2!(4-2)!} - \left(\frac{3!}{1!(3-1)!} - \frac{2!}{1!(2-1)!} + \frac{4!}{2!(4-2)!} - \frac{4!}{2!(4-2)!}\right)\right)}{\frac{4!}{2!(4-2)!}} \ .$$

In turn, these two indices are the coordinates of the sequence in a bi-dimensional space and the distance between two binary sequences *(xi, x)* is the Euclidean distance between a couple of lexicographic indices

$$r(\mathbf{x}\_{l}, \mathbf{x}\_{\ell}) = \sqrt{\left(d'(\mathbf{x}\_{l}) - d'(\mathbf{x}\_{\ell})\right)^{2} + \left(d''(\mathbf{x}\_{l}) - d''(\mathbf{x}\_{\ell})\right)^{2}}\,. \tag{3}$$

2Duration and timing are two of the three aspects identified by Studer and Ritschard (2016) as mattering in sequence comparison. Here, we do not pay attention to the third one, sequencing, which is not a concern for studying (de)synchronization.

Passing from a binary sequence to a multinomial sequence is easy. Just as a qualitative variable of *m* modality can be represented by *m* dummy variables, so a multinomial sequence of alphabet *Q* can be represented by |*Q*| binary sequences *xq* with values 0–1. For example, the sequence *x* = [123321] and alphabet *Q* = {1*,* 2*,* 3} can be represented by the following three binary sequences *x*1 = [100001]; *x*2 = [010010]; *x*3 = [001100]. To each of these binary sequences it is possible to apply the lexicographic index and compute the coordinates {*d q (xq )*; *d q (xq )*}. The multinomial sequence *x* is therefore described by a vector of real numbers. The distance between two multinomial sequences *(xi, x)* is the Euclidean distance between their transformations {*d (xiq )*; *d(xiq )*} and {*d (xq )*; *d(xq )*}. Formally, it is:

$$r'(\mathbf{x}\_l, \mathbf{x}\_\ell) = \sqrt{\sum\_{q=1}^{|\mathcal{Q}|} \left( d'(\mathbf{x}\_{lq}) - d'(\mathbf{x}\_{\ell q}) \right)^2 + \left( d''(\mathbf{x}\_{lq}) - d''(\mathbf{x}\_{\ell q}) \right)^2} \,. \tag{4}$$