Indicateurs synthétiques

A titre d’exemple, on construit une matrice de distance à partir d’indicateurs décrivant le nombre d’épisodes dans les différents états.

indics <- seqinepi(seqact)
head(indics)

  nepi1 nepi2 nepi3 nepi4 nepi5 nepi6
1     0     1     0     0     0     0
2     1     1     2     0     0     0
3     1     2     1     0     2     0
4     1     2     0     0     0     1
5     1     1     0     0     0     0
6     1     2     0     0     1     0

La matrice peut être calculée directement à partir des indicateurs ou après une étape d’analyse en composantes principales (ACP), ici en retenant les 5 premières dimensions.

# matrice de distance à partir des indicateurs
dissim <- dist(indics, method='euclidean') %>% as.matrix

# matrice de distance à partir des résultats d'une ACP
acp_coords <- PCA(indics, scale.unit=FALSE, ncp=5, graph=FALSE)$ind$coord
dissim <- dist(acp_coords, method='euclidean') %>% as.matrix

D’autres indicateurs synthétiques (durées, états visités, etc.) peuvent être calculés simplement à partir des fonctions seqistatd, seqi1epi, seqifpos, seqindic ou seqpropclust.

Disjonctif complet et AHQ

Dans le cas du codage sous forme de disjonctif complet, la matrice de distance peut être calculée directement, avec la distance euclidienne ou la distance du chi2, ou après une étape d’analyse en composantes principales (ACP) ou d’analyse des correspondances multiples (ACM), ici en retenant les 5 premières dimensions.

# codage disjonctif complet
disjo <- dichotom(seqact)
disjo <- disjo[,colSums(disjo)>0]

# distance euclidienne
dissim <- dist(disjo, method='euclidean') %>% as.matrix

# distance du chi2
dissim <- map_df(disjo, as.factor) %>%
          dudi.acm(scannf=FALSE, nf=ncol(disjo)) %>%
          dist.dudi() %>%
          as.matrix

# après une ACP
acp_coords <- PCA(disjo, scale.unit=FALSE, ncp=5, graph=FALSE)$ind$coord
dissim <- dist(acp_coords, method='euclidean') %>% as.matrix

# après une ACM
acm_res <- purrr::map_df(disjo, as.factor) %>%
           MCA(ncp=5, graph=FALSE)
dissim <- dist(acm_res$ind$coord, method='euclidean') %>% as.matrix

Pour l’analyse harmonique qualitative (AHQ), le calcul de la matrice de distance peut se faire directement (distance du chi2) ou après une analyse factorielle des correspondances (AFC), ici en retenant les 5 premières dimensions.

# codage AHQ
ahq <- seq2qha(seqact, c(1,3,7,10,15,20,28))
ahq <- ahq[,colSums(ahq)>0]

# distance du chi2
dissim <- dudi.coa(ahq, scannf=FALSE, nf=ncol(ahq)) %>%
          dist.dudi() %>%
          as.matrix

# après une AFC
afc_coord <- CA(ahq, ncp=5, graph=FALSE)$row$coord
dissim <- dist(afc_coord, method='euclidean') %>% as.matrix

Optimal Matching et alternatives

Pour l’Optimal Matching, la construction d’une matrice de distance entre les séquences s’effectue avec la fonction seqdist. Cela implique de définir également une matrice de coûts de substitution entre les états (avec la fonction seqsubm). Ici, les coûts de substitution sont constants et égaux à 2 et le coût indel est égal à 1,5.

# construction de la matrice de distance
couts <- seqsubm(seqact, method="CONSTANT", cval=2)
dissim <- seqdist(seqact, method="OM", sm=couts, indel=1.5)

D’expérience, l’Optimal Matching avec le paramétrage adopté ici constitue un choix permettant de prendre en compte les différentes dimensions de la temporalité des séquences - ordonnancement (sequencing), calendrier (timing), durée (dans les différents états ou des épisodes, duration et spell duration). Si on souhaite privilégier l’une de ces dimensions, on peut suivre les recommandations de Studer & Ritschard (2016, voir en particulier pages 507-509), et choisir l’une des nombreuses autres métriques implémentées dans l’extension TraMineR.

# sequencing
dissim <- seqdist(seqact, method="OMstran", otto=0.1, sm=couts, indel=1)
dissim <- seqdist(seqact, method="OMspell", expcost=0, sm=couts, indel=1)
dissim <- seqdist(seqact, method="SVRspell", tpow=0)

# timing
dissim <- seqdist(seqact, method="HAM", sm=couts)
dissim <- seqdist(seqact, method="CHI2", step=1)

# duration (distribution aver the entire period)
dissim <- seqdist(seqact, method="EUCLID", step=37)

# duration (spell lengths)
dissim <- seqdist(seqact, method="OMspell", expcost=1, sm=couts, indel=1)
dissim <- seqdist(seqact, method="LCS")

A noter, les méthodes passant par le disjonctif complet ou l’AHQ sont également implémentées dans la fonction seqdist (méthodes “EUCLID” et “CHI2”).

Construction d’une typologie

On réalise ensuite une classification ascendante hiérarchique (CAH) avec le critère d’agrégation de Ward, à l’aide de la fonction agnes de l’extension cluster.

NB : Avec un nombre élevé de séquences, la CAH peut nécessiter un temps de calcul important. Il existe cependant une implémentation nettement plus rapide dans l’extension fastcluster (fonction hclust).

# classification ascendante hiérarchique
agnes <- as.dist(dissim) %>% agnes(method="ward", keep.diss=FALSE)

Pour explorer les solutions d’une classification ascendante hiérarchique, on commence généralement par examiner le dendrogramme.

# dendrogramme
as.dendrogram(agnes) %>% plot(leaflab="none")

Le graphique suivant combine dendrogramme et index plot: les séquences de l’index plot sont triées selon leur position dans le dendrogramme, lui-même représenté en marge du graphique.

# heatmap (dendrogramme + index plot)
seq_heatmap(seqact, agnes)

L’examen des sauts d’inertie peut également être utile pour déterminer le nombre de classes de la typologie. On voit par exemple qu’il y a une différence d’inertie notable entre les partitions en 5 et 6 classes.

# graphique des sauts d'inertie
plot(sort(agnes$height, decreasing=TRUE)[1:20], type="s", xlab="nombre de classes", ylab="inertie")

Il existe aussi un certain nombre d’indicateurs de qualité des partitions (silhouette, Calinski-Harabasz, pseudo-R2, etc.).

# indicateurs de qualité
wardRange <- as.clustrange(agnes, diss=dissim)
summary(wardRange, max.rank=2)

     1. N groups     1.  stat 2. N groups     2.  stat
PBC            5   0.79612318           4   0.79012613
HG             5   0.92899997           4   0.92259604
HGSD           5   0.92672338           4   0.92027518
ASW            3   0.53491169           2   0.52958820
ASWw           3   0.53862662           2   0.53182999
CH             2 156.86409716           3 103.74539004
R2            20   0.60648740          19   0.59473652
CHsq           2 324.53910100           3 245.89231855
R2sq          20   0.82862437          19   0.82073439
HC             5   0.04010082           4   0.04392402

On représente ici graphiquement la qualité des partitions pour différents nombres de classes pour les indicateurs silhouette, pseudo-R2 et Calinski-Harabasz.

plot(wardRange, stat=c('ASW','R2','CH'), norm="zscore")

On opte au final pour une partition en 5 classes, en “coupant l’arbre” de la CAH à l’aide de la fonction cutree.

# choix de la partition en 5 classes
nbcl <- 5
part <- cutree(agnes, nbcl)

Il est possible de “consolider” la partition à l’aide de l’algorithme PAM (Partition Around Medoids) et de la fonction wcKMedoids de l’extension WeightedCluster. On aboutit ainsi à une distribution des séquences entre les classes très similaire (cf le tri croisant les classes avant et après consolidation) mais la qualité de la partition consolidée est légèrement supérieure (le R² passe de 61 à 64%).

# consolidation de la partition
newpart <- wcKMedoids(dissim, k=nbcl, initialclust=part, cluster.only=TRUE)
table(part, newpart)

    newpart
part  28  34 103 307 480
   1 104 241   1   0   0
   2   9  24   2   7   0
   3   0   6   6  63   1
   4   0   0  17   0   0
   5   8   0   0   1  10

wcClusterQuality(dissim, part)$stats['R2sq'] %>% round(3)

 R2sq 
0.607 

wcClusterQuality(dissim, newpart)$stats['R2sq'] %>% round(3)

R2sq 
0.64 

Si on souhaite conserver la partition consolidée :

part <- as.numeric(as.factor(newpart))

NB: Autre option, la classification “floue”, ici avec l’algorithme Fanny (extension cluster). A la différence de la CAH ou de PAM, chaque séquence n’appartient pas à une seule classe mais se caractérise par des degrés d’appartenance aux différentes classes. Le tableau suivant présente les degrés d’appartenance aux 6 classes des 3 premières séquences du corpus. La première séquence appartient à 99% à la classe 1, mais la seconde est plus “partagée”, principalement entre les classes 2 et 5.

# classification floue (fuzzy clustering)
fanny <- as.dist(dissim) %>% fanny(k=5, metric='euclidean', memb.exp=1.2)
fanny$membership %>% round(2) %>% .[1:3,]

     [,1] [,2] [,3] [,4] [,5]
[1,] 0.99 0.00 0.00 0.01 0.00
[2,] 0.06 0.51 0.02 0.13 0.28
[3,] 0.00 0.11 0.80 0.00 0.08

Description de la typologie: graphiques

Les représentations graphiques permettent de se faire une idée rapide et intuitive de la nature des classes de la typologie. Le type de graphique le plus utilisé est le chronogramme (state distribution plot).

# chronogrammes de la typologie
seqdplot(seqact, group=part, xtlab=14:50, border=NA, cex.legend=0.8)

Les index plots (ou “tapis”) sont également très répandus.

# index plots de la typologie
seqIplot(seqact, group=part, xtlab=14:50, yaxis=FALSE, cex.legend=0.8)

Les index plots sont souvent plus faciles à interpréter lorsqu’on trie les séquences, en particulier à partir d’une procédure de multidimensional scaling.

# index plots de la typologie, triés par multidimensional scaling
mds.order <- cmdscale(dissim,k=1)
seqIplot(seqact, sortv=mds.order, group=part, xtlab=14:50, yaxis=FALSE, cex.legend=0.8)

Ils peuvent également être “lissés”, pour les rendre plus lisibles.

Méthode des “smoothed MDS sequence plots” (Piccarreta, 2012):

smoothed <- seqsmooth(seqact, dissim, k=30)$seqdata
seqIplot(smoothed, sortv=mds.order, group=part, xtlab=14:50, yaxis=FALSE, cex.legend=0.8)

Méthode des “relative frequency sequence plots” (Fasang & Liao, 2013):

# relative frequency sequence plots
seqplot.rf(seqact, diss=dissim, group=part, xtlab=14:50)

Les sequence frequency plots représentent, pour chaque classe, les 10 séquences les plus fréquentes (avec une épaisseur proportionnelle à leur fréquence).

# frequency plots
seqfplot(seqact, group=part, ylab="", xtlab=14:50, cex.legend=0.8)

Les modal state sequence plots représentent, pour chaque classe, la séquence des états modaux pour chaque position dans le temps. A chaque position dans le temps, la hauteur de la barre est proportionnelle à la fréquence de l’état modal.

# modal state plots
seqmsplot(seqact, group=part, xtlab=14:50, cex.legend=0.8)

Description de la typologie: statistiques

La première étape de description statistique de la typologie consiste généralement à présenter le poids des classes. La classe 2 rassemble plus de la moitié des individus, alors que les classes 3 et 5 sont de très petite taille.

# effectifs
table(part)

part
  1   2   3   4   5 
121 271  26  71  11 

# pourcentages
100*(table(part))/length(part)

part
   1    2    3    4    5 
24.2 54.2  5.2 14.2  2.2 

Il est utile d’évaluer l’homogénéité des classes. Cela peut se faire à partir des distances intra-classes. Les classes 1 et 2 sont les plus homogènes.

# distances intra-classes
Dintra <- integer(length=nbcl)
for(i in 1:nbcl) Dintra[i] <- round(mean(dissim[part==i,part==i]),1)
Dintra

[1] 18.3 15.2 30.2 28.4 24.3

Les résultats sont convergents à partir des distances moyennes aux centres de classe :

# distances moyennes au centre de la classe
dissassoc(dissim, part)$groups

        n discrepancy
1     121    9.171231
2     271    7.617135
3      26   15.078402
4      71   14.208490
5      11   12.173554
Total 500   15.876208

De même avec l’entropie transversale :

# entropie transversale moyenne par classe
entropie <- vector()
for(i in 1:nbcl) entropie[i] <- round(mean(seqstatd(seqact[part==i,])$Entropy),2)
entropie

[1] 0.33 0.27 0.47 0.46 0.36

Pour donner une vision plus détaillée de la forme des trajectoires de chaque classe, on calcule des indicateurs synthétiques, puis leur moyenne selon la classe de la typologie. La durée passée dans les états :

# durées dans les états selon la classe
dur <- seqistatd(seqact)
durees <- round(aggregate(dur, by=list(part), FUN=mean), 1)
rownames(durees) <- NULL
durees

  Group.1    1    2    3    4    5   6
1       1 10.2 23.1  0.7  0.7  1.3 1.0
2       2  3.5 29.9  0.4  0.9  1.7 0.6
3       3  6.3  9.2 16.4  1.0  4.0 0.0
4       4  5.0  7.1  1.2  1.6 22.0 0.1
5       5  7.0  3.6  0.6 22.8  2.7 0.2

La part d’individus ayant connu au moins un épisode dans les états :

# au moins un épisode dans les états, selon la classe
epi <- seqi1epi(seqact)
episodes <- round(aggregate(epi, by=list(part), FUN=mean), 2)
rownames(episodes) <- NULL
episodes

  Group.1 epi1 epi2 epi3 epi4 epi5 epi6
1       1 1.00 0.99 0.16 0.14 0.32 0.36
2       2 0.86 1.00 0.06 0.17 0.31 0.35
3       3 0.92 0.73 1.00 0.15 0.62 0.04
4       4 0.86 0.75 0.24 0.24 1.00 0.03
5       5 0.91 0.82 0.27 1.00 0.36 0.09

Ensuite, on croise la typologie avec les caractéristiques des individus. On commence par une analyse bivariée du type de trajectoire selon le sexe. Le V² de Cramer est de 0,16, indiquant une association notable. Les femmes sont très sur-représentées dans la classe 4 et secondairement dans la classe 3, alors que les hommes sont sur-représentés dans les classes 1 et 2 (les valeurs “phi” correspondent aux attractions ou répulsions entre modalités).

assoc.twocat(factor(part), socdem$sexe)

$freq
    homme femme Sum
1      72    49 121
2     145   126 271
3       2    24  26
4       3    68  71
5       4     7  11
Sum   226   274 500

$prop
    homme femme   Sum
1    14.4   9.8  24.2
2    29.0  25.2  54.2
3     0.4   4.8   5.2
4     0.6  13.6  14.2
5     0.8   1.4   2.2
Sum  45.2  54.8 100.0

$rprop
    homme femme Sum
1    59.5  40.5 100
2    53.5  46.5 100
3     7.7  92.3 100
4     4.2  95.8 100
5    36.4  63.6 100
Sum  45.2  54.8 100

$cprop
    homme femme   Sum
1    31.9  17.9  24.2
2    64.2  46.0  54.2
3     0.9   8.8   5.2
4     1.3  24.8  14.2
5     1.8   2.6   2.2
Sum 100.0 100.0 100.0

$v2.cramer
[1] 0.1615622

$phi
   homme  femme
1  0.162 -0.162
2  0.182 -0.182
3 -0.177  0.177
4 -0.335  0.335
5 -0.027  0.027

On examine ensuite les sur-représentations pour chaque classe de la typologie à partir de l’ensemble des caractéristiques individuelles présentes dans le tableau de données socdem. On constate tout d’abord que seules trois variables ne semblent pas notablement liées au type de trajectoire (mereactive, nbunion, nationalite). Pour ne prendre que l’exemple de la classe 4, on voit que les femmes y sont sur-représentées, ainsi que les individus ayant trois enfants ou plus, sans diplôme ou de PCS inconnue (ce qui n’est certainement pas sans lien avec l’inactivité).

catdesc(factor(part), socdem, min.phi=0.1)

$variables
     variable measure assoc p.value    criterion
1         PCS  cramer 0.362 0.00000 0.000000e+00
2     diplome  cramer 0.330 0.00000 0.000000e+00
3        sexe  cramer 0.402 0.00000 0.000000e+00
4     PCSpere  cramer 0.202 0.00000 3.423300e-06
5      annais    eta2 0.033 0.00376 3.767329e-03
6       nbenf  cramer 0.146 0.01237 1.245175e-02
7  mereactive  cramer 0.130 0.51600 7.256666e-01
8     nbunion  cramer 0.094 0.98143 3.986295e+00
9 nationalite  cramer 0.057 1.00000 1.476386e+01

$bylevel
$bylevel$`1`
$bylevel$`1`$categories
       categories pct.ycat.in.xcat pct.xcat.in.ycat pct.xcat.global    phi
31   diplome.>bac            0.598            0.628           0.254  0.480
15       PCS.pcis            0.579            0.512           0.214  0.405
22   PCSpere.pcis            0.519            0.331           0.154  0.270
10     sexe.homme            0.319            0.595           0.452  0.158
1     nbenf.aucun            0.367            0.182           0.120  0.100
12         PCS.NA            0.111            0.066           0.144 -0.119
26 PCSpere.onqual            0.058            0.025           0.104 -0.139
28  diplome.aucun            0.075            0.041           0.134 -0.147
11     sexe.femme            0.179            0.405           0.548 -0.158
17       PCS.empl            0.107            0.107           0.242 -0.172
29   diplome.<bac            0.057            0.107           0.454 -0.389

$bylevel$`1`$continuous.var
  variables median.x.in.ycat median.x.global sd.x.in.ycat sd.x.global corr.coef
1    annais             1945            1943     5.743567    6.028132 0.0804953


$bylevel$`2`
$bylevel$`2`$categories
      categories pct.ycat.in.xcat pct.xcat.in.ycat pct.xcat.global    phi
29  diplome.<bac            0.722            0.605           0.454  0.326
18     PCS.oqual            0.878            0.133           0.082  0.194
10    sexe.homme            0.642            0.535           0.452  0.178
17      PCS.empl            0.678            0.303           0.242  0.149
25 PCSpere.oqual            0.670            0.262           0.212  0.128
11    sexe.femme            0.460            0.465           0.548 -0.178
12        PCS.NA            0.306            0.081           0.144 -0.189
15      PCS.pcis            0.308            0.122           0.214 -0.240
22  PCSpere.pcis            0.195            0.055           0.154 -0.292
31  diplome.>bac            0.228            0.107           0.254 -0.363

$bylevel$`2`$continuous.var
  variables median.x.in.ycat median.x.global sd.x.in.ycat sd.x.global
1    annais             1942            1943      5.98217    6.028132
   corr.coef
1 -0.1019808


$bylevel$`3`
$bylevel$`3`$categories
   categories pct.ycat.in.xcat pct.xcat.in.ycat pct.xcat.global    phi
11 sexe.femme            0.088            0.923           0.548  0.167
10 sexe.homme            0.009            0.077           0.452 -0.167

$bylevel$`3`$continuous.var
  variables median.x.in.ycat median.x.global sd.x.in.ycat sd.x.global corr.coef
1    annais             1947            1943     4.254229    6.028132 0.1394812


$bylevel$`4`
$bylevel$`4`$categories
         categories pct.ycat.in.xcat pct.xcat.in.ycat pct.xcat.global    phi
12           PCS.NA            0.542            0.549           0.144  0.461
11       sexe.femme            0.248            0.958           0.548  0.329
4  nbenf.trois ou +            0.235            0.535           0.324  0.177
28    diplome.aucun            0.239            0.225           0.134  0.101
16         PCS.pint            0.065            0.113           0.248 -0.121
15         PCS.pcis            0.009            0.014           0.214 -0.191
10       sexe.homme            0.013            0.042           0.452 -0.329

$bylevel$`4`$continuous.var
  variables median.x.in.ycat median.x.global sd.x.in.ycat sd.x.global
1    annais             1942            1943     6.502561    6.028132
    corr.coef
1 -0.07858364


$bylevel$`5`
$bylevel$`5`$categories
   categories pct.ycat.in.xcat pct.xcat.in.ycat pct.xcat.global   phi
13  PCS.agric                1            0.091           0.002 0.146

$bylevel$`5`$continuous.var
  variables median.x.in.ycat median.x.global sd.x.in.ycat sd.x.global
1    annais             1948            1943     5.989385    6.028132
   corr.coef
1 0.08723237

Description de la typologie: parangons

Pour “incarner” la typologie, on recourt parfois aux parangons (medoids) des classes, dont on retrace de manière détaillée les trajectoires à partir d’informations non prises en compte dans le codage des trajectoires.

# parangon de chaque classe (numéros de ligne dans le fichier de données)
medoids(dissim, part)

[1]  28  34 103 307 480

Distance à une séquence de référence

On définit ici une séquence “de référence” correspondant à une trajectoire d’emploi à temps plein en continu à partir de 18 ans. On calcule ensuite, pour chaque séquence, sa distance à la séquence de référence: dans quelle mesure s’écartent-elles de cette référence ?

ref <- seqdef(as.matrix("(1,4)-(2,33)"), informat="SPS", alphabet=alphabet(seqact))
distref <- seqdist(seqact, refseq = ref, method="OM", sm=couts, indel=1.5)

On observe ensuite la distribution de ces écarts en fonction du sexe et du nombre d’enfants. On constate que les trajectoires des femmes s’écartent plus de la trajectoire d’emploi continu à temps plein que celles des hommes lorsqu’elles ou ils ont un ou plusieurs enfants. L’écart est particulièrement fort chez les femmes ayant trois enfants ou plus.

socdem %>% select(sexe,nbenf) %>%
           tibble(distref=distref) %>%
           ggplot(aes(x=nbenf, y=distref)) + 
             geom_boxplot(aes(fill=sexe), notch=T) +
             xlab("nombre d'enfants") +
             ylab("distance à la référence") +
             theme_bw()

Distances inter ou intra-groupes

On compare les distances entre les trajectoires féminines et les distances entre trajectoires masculines. Les trajectoires masculines sont nettement plus homogènes ou, pour le dire autrement, celles des femmes sont plus diversifiées.

# distances intra-classes selon le sexe
sapply(levels(socdem$sexe), function(x) round(mean(dissim[socdem$sexe==x,socdem$sexe==x]),1))

homme femme 
 18.8  38.1 

La matrice de distance entre trajectoires peut être résumé graphiquement à partir d’un multidimensional scaling (MDS). On représente ici le nuage des individus dans le plan formé par les deux premiers axes, en colorant les points selon leur classe, puis on projette le sexe en variable supplémentaire. On observe par exemple que, sur le premier axe, les classes 1 et 2 s’opposent aux autres, et que les hommes sont du côté des classes 1 et 2.

mds <- cmdscale(dissim, k=2)
plot(mds, type='n', xlab="axe 1", ylab="axe 2")
abline(h=0, v=0, lty=2, col='lightgray')
points(mds, pch=20, col=part)
legend('topleft', paste('classe',1:nbcl), pch=20, col=1:nbcl, cex=0.8)
text(aggregate(mds, list(socdem$sexe), mean)[,-1], levels(socdem$sexe), col='orange', cex=1, font=2)

Indicateurs synthétiques

On peut étudier la distribution d’indicateurs synthétisant les trajectoires en fonction d’autres caractéristiques des individus. Par exemple, le temps passé dans les différents états selon le sexe : on voit apparaître le poids de l’inactivité dans les trajectoires féminines.

# durées dans les états selon le sexe
dur <- seqistatd(seqact)
durees_sexe <- aggregate(dur, by=list(socdem$sexe), function(x) round(mean(x),1))
rownames(durees_sexe) <- NULL
colnames(durees_sexe) <- c("classe",labs)
durees_sexe

  classe études temps plein temps partiel petits boulots inactivité
1  homme    5.5        27.7           0.4            1.2        0.9
2  femme    5.6        19.8           2.3            1.6        7.7
  serv. militaire
1             1.3
2             0.0

On donne maintenant un exemple d’indicateur de complexité des trajectoires, la turbulence, croisé avec l’année de naissance. On ne constate pas d’évolution notable.

# turbulence
turbu <- aggregate(seqST(seqact), list(socdem$annais), mean)
plot(turbu, type='l', ylim=c(0,10), xlab='Année de naissance')

D’autres indicateurs de complexité des trajectoires peuvent être calculés facilement : indice de complexité (fonction seqici), entropie individuelle (fonction seqient), nombre de transitions (fonction seqtransn).

Analyse de variance

Avec une variable explicative unique, on mesure la part de la variance des dissimilarités expliquée par la variable (à l’aide d’un pseudo-R2), ainsi qu’une mesure de la variabilité des trajectoires pour chacune des modalités de la variable, i.e. dans chaque sous-population. Dans notre exemple, le sexe explique 7,4% de la variance des distances entre trajectoires d’emploi, et la variabilité des trajectoires est nettement plus élevée chez les femmes que chez les hommes (19,1 contre 9,4).

# analyse de variance avec le sexe comme facteur
dissassoc(dissim, socdem$sexe)

Pseudo ANOVA table:
             SS  df       MSE
Exp    591.0611   1 591.06114
Res   7347.0429 498  14.75310
Total 7938.1040 499  15.90802

Test values  (p-values based on 1000 permutation):
                     t0 p.value
Pseudo F    40.06352676   0.001
Pseudo Fbf  42.94304887   0.001
Pseudo R2    0.07445873   0.001
Bartlett    29.68308216   0.001
Levene     123.83930458   0.001

Inconclusive intervals: 
0.00383  <  0.01  <  0.0162
0.03649  <  0.05  <  0.0635 

Discrepancy per level:
        n discrepancy
homme 226    9.383918
femme 274   19.074005
Total 500   15.876208

Il est possible de détailler ces indicateurs pour chaque position dans le temps des trajectoires. Ainsi, la part de variance expliquée par le sexe est presque nulle en début de trajectoire, elle croît ensuite fortement entre 18 ans et environ 30 ans (elle est alors de 14%), puis diminue pour n’être plus que de 5% à 50 ans.

# analyse de variance selon la position dans le temps
diff <- seqdiff(seqact, group=socdem$sexe)
rownames(diff$stat) <- rownames(diff$discrepancy) <- 14:49
plot(diff, stat="Pseudo R2")

La variabilité des trajectoires des femmes et des hommes est faible en début de trajectoire et augmente jusqu’à l’âge de 21 ans. Les résultats divergent ensuite selon le sexe: la variabilité des trajectoires féminines se maintient jusqu’à 50 ans, alors que celle des hommes diminue fortement après 21 ans et atteint un niveau très faible entre 30 et 50 ans.

pal <- brewer.pal(ncol(diff$discrepancy), "Set2")
plot(diff, stat="discrepancy", legend.pos=NA, col=pal, lwd=1.5)
legend('topright', fill=pal, legend=colnames(diff$discrepancy), cex=0.7)

Avec plusieurs variables explicatives, on obtient la part de variance expliquée par l’ensemble des variables et la décomposition de cette part entre les variables. Ici, l’année de naissance, le sexe, le niveau de diplôme et le nombre d’enfants expliquent ensemble 16,2% de la variance des dissimilarités entre les trajectoires d’emploi: 7,7% pour le sexe, 6,0% pour le diplôme, 2,1% pour le nombre d’enfants et 0,3% pour l’année de naissance.

# analyse de la variance avec facteurs multiples
dissmfacw(dissim ~ annais+nbenf+sexe+diplome, data=socdem)

  Variable   PseudoF    PseudoR2 p_value
1   annais  1.939020 0.003308502   0.042
2    nbenf  4.147849 0.021232117   0.001
3     sexe 45.251273 0.077211134   0.001
4  diplome 11.654366 0.059656674   0.001
5    Total 11.883983 0.162218781   0.001

Inconclusive intervals: 
0.00383  <  0.01  <  0.0162
0.03649  <  0.05  <  0.0635 

Enfin, l’analyse de variance peut servir à construire un arbre de décision, dit aussi arbre d’induction.

# arbre d'induction
arbre <- seqtree(seqact ~ annais+nbenf+sexe+diplome, data=socdem, diss=dissim, min.size=0.1, max.depth=3)

L’arbre peut être présenté sous forme textuelle ou graphique. A noter, la représentation graphique nécessite l’installation du logiciel GraphViz (cf l’aide de la function seqtreedisplay).

Comme précédemment, on constate que la variable qui explique la plus grande part de variance des dissimilarités est le sexe. Ensuite, pour les femmes, la variable la plus discriminante est le nombre d’enfants et, plus précisément, le fait d’avoir ou non trois enfants ou plus. L’inactivité et, secondairement, le temps partiel sont plus présents dans les trajectoires des femmes ayant au moins trois enfants que dans celles des autres. Chez les hommes, en revanche, c’est le niveau de diplôme qui est le plus discriminant (les hommes entrent plus tard sur le marché du travail lorsqu’ils ont un diplôme supérieur ou égal au baccalauréat).

# résultats de l'arbre sous forme textuelle
print(arbre)

Dissimilarity tree:
 Parameters: min.size=50, max.depth=3, R=1000, pval=0.01 
 Formula: seqact ~ annais + nbenf + sexe + diplome 
 Global R2: 0.12934 

 Fitted tree: 

 |-- Root  (n: 500 disc: 15.876) 
   |-> sexe 0.074459 
             |-- [ homme ]  (n: 226 disc: 9.3839) 
               |-> diplome 0.093027 
                   |-- [ aucun,<bac ]    (n: 131 disc: 8.5907)[(1,4)-(2,3)-(6,1)-(2,29)] * 
                   |-- [ bac,>bac ]    (n: 95 disc: 8.401)[(1,9)-(2,28)] * 
             |-- [ femme ]  (n: 274 disc: 19.074) 
               |-> nbenf 0.045614 
                   |-- [ trois ou + ]    (n: 86 disc: 21.168)[(1,4)-(2,10)-(5,13)-(2,10)] * 
                   |-- [ aucun,un,deux ]    (n: 188 disc: 16.848)[(1,6)-(2,31)] * 

# résultats de l'arbre sous forme graphique
seqtreedisplay(arbre,type="d",border=NA,show.depth=TRUE)

Statistiques implicatives

Pour étudier la manière dont les trajectoires diffèrent entre plusieurs sous-populations, Studer propose d’utiliser les statistiques implicatives. Il s’agit de reconstituer, pour chaque population, la séquence des états typiques (Studer, 2012; Struffolino et al, 2016).

Ici, le service militaire est typique des trajectoires des hommes autour de 20 ans, puis c’est l’emploi à temps plein à partir de 25 ans (Figure 10). L’inactivité est caractéristique des trajectoires des femmes, et cela dès l’âge de 18 ans. L’emploi à temps partiel l’est également, mais de manière moins marquée et à partir de 30 ans.

# statistiques implicatives
implic <- seqimplic(seqact, group=socdem$sexe)
# par(mar=c(2,2,2,2))
plot(implic, xtlab=14:50, lwd=2, conf.level=c(0.95, 0.99), cex.legend=0.7)

Mise en route

On commence par charger les données correspondant à trois dimensions des trajectoires de 500 individus tirés au sort dans l’enquête Biographies et entourage - matrimoniale, parentale et résidentielle - et observées de 14 à 35 ans.

data(seqmsa)
trajmat <- seqmsa %>% select(starts_with('smat'))
str(trajmat)

'data.frame':   500 obs. of  22 variables:
 $ smat14: num  1 1 1 1 1 1 1 1 1 1 ...
 $ smat15: num  1 1 1 1 1 1 1 1 1 1 ...
 $ smat16: num  1 1 1 1 1 1 1 1 1 1 ...
 $ smat17: num  1 1 1 1 1 1 1 1 1 1 ...
 $ smat18: num  1 1 1 1 1 1 1 1 1 1 ...
 $ smat19: num  1 1 1 1 1 1 1 1 1 1 ...
 $ smat20: num  1 3 1 1 3 1 1 1 1 1 ...
 $ smat21: num  1 3 1 1 3 3 1 3 1 1 ...
 $ smat22: num  1 3 1 3 3 3 1 3 1 1 ...
 $ smat23: num  1 3 1 3 3 3 1 3 1 1 ...
 $ smat24: num  1 3 1 3 3 3 1 3 1 1 ...
 $ smat25: num  1 3 3 3 3 3 1 3 2 1 ...
 $ smat26: num  2 3 3 3 3 3 1 3 2 1 ...
 $ smat27: num  2 3 3 3 3 3 1 3 2 1 ...
 $ smat28: num  2 3 3 3 3 3 1 3 3 1 ...
 $ smat29: num  2 3 3 3 3 3 1 3 3 1 ...
 $ smat30: num  2 3 3 3 3 3 1 3 3 1 ...
 $ smat31: num  3 3 3 3 3 3 1 3 3 1 ...
 $ smat32: num  3 3 3 3 2 3 1 3 3 1 ...
 $ smat33: num  3 3 3 3 2 3 3 3 3 1 ...
 $ smat34: num  3 3 3 3 2 3 3 3 3 1 ...
 $ smat35: num  3 3 3 3 2 3 3 3 3 1 ...

trajenf <- seqmsa %>% select(starts_with('nenf'))
str(trajenf)

'data.frame':   500 obs. of  22 variables:
 $ nenf14: num  0 0 0 0 0 0 0 0 0 0 ...
 $ nenf15: num  0 0 0 0 0 0 0 0 0 0 ...
 $ nenf16: num  0 0 0 0 0 0 0 0 0 0 ...
 $ nenf17: num  0 0 0 0 0 0 0 0 0 0 ...
 $ nenf18: num  0 0 0 0 0 0 0 0 0 0 ...
 $ nenf19: num  0 0 0 0 0 0 0 0 0 0 ...
 $ nenf20: num  0 0 0 0 1 0 0 0 0 0 ...
 $ nenf21: num  0 1 0 0 1 0 0 0 0 0 ...
 $ nenf22: num  0 1 0 0 1 0 0 0 0 0 ...
 $ nenf23: num  0 1 0 0 1 0 0 0 0 0 ...
 $ nenf24: num  0 1 0 0 1 1 0 0 0 0 ...
 $ nenf25: num  0 1 0 1 1 1 0 1 0 0 ...
 $ nenf26: num  0 1 0 1 1 1 0 1 2 0 ...
 $ nenf27: num  0 1 0 1 1 1 0 1 2 0 ...
 $ nenf28: num  0 1 1 1 1 1 0 2 3 0 ...
 $ nenf29: num  0 2 1 2 1 1 0 2 3 0 ...
 $ nenf30: num  0 2 2 2 1 1 0 3 3 0 ...
 $ nenf31: num  1 2 2 2 1 1 0 3 3 0 ...
 $ nenf32: num  1 2 2 2 1 1 0 3 3 0 ...
 $ nenf33: num  1 2 2 2 2 1 0 3 3 0 ...
 $ nenf34: num  1 2 2 2 2 1 0 3 3 0 ...
 $ nenf35: num  1 2 2 2 2 1 0 3 3 0 ...

trajlog <- seqmsa %>% select(starts_with('slog'))
str(trajlog)

'data.frame':   500 obs. of  22 variables:
 $ slog14: num  0 0 0 0 0 0 0 0 0 0 ...
 $ slog15: num  0 0 0 0 0 0 0 0 0 0 ...
 $ slog16: num  0 0 0 0 0 0 0 0 0 0 ...
 $ slog17: num  0 0 0 0 0 0 0 0 1 0 ...
 $ slog18: num  0 0 0 0 0 0 0 0 1 0 ...
 $ slog19: num  0 0 0 0 0 1 0 0 1 0 ...
 $ slog20: num  0 0 0 0 1 1 0 0 0 0 ...
 $ slog21: num  0 0 0 0 1 0 0 0 0 0 ...
 $ slog22: num  0 1 0 1 1 0 0 0 1 0 ...
 $ slog23: num  0 1 0 1 1 1 0 1 1 0 ...
 $ slog24: num  0 1 0 1 1 1 0 1 1 0 ...
 $ slog25: num  0 1 1 1 1 1 0 1 1 0 ...
 $ slog26: num  1 1 1 1 1 1 0 1 1 0 ...
 $ slog27: num  1 1 1 1 1 1 0 1 1 1 ...
 $ slog28: num  1 1 1 1 1 1 0 1 1 1 ...
 $ slog29: num  1 1 1 1 1 1 0 1 1 1 ...
 $ slog30: num  1 1 1 1 1 1 0 1 1 1 ...
 $ slog31: num  1 1 1 1 1 1 0 1 1 1 ...
 $ slog32: num  1 1 1 1 1 1 0 1 1 1 ...
 $ slog33: num  1 1 1 1 1 1 0 1 1 1 ...
 $ slog34: num  1 1 1 1 1 1 0 1 1 1 ...
 $ slog35: num  1 1 1 1 1 1 1 1 1 1 ...

On définit ensuite trois objets séquences (un par dimension).

# définition de la trajectoire matrimoniale
labs <- c("jamais","union libre","marié","separé")
palette <- brewer.pal(length(labs), 'Set2')
seqmat <- seqdef(trajmat, labels=labs, cpal=palette)

# définition de la trajectoire parentale
labs <- c("0","1","2","3+")
palette <- brewer.pal(length(labs), 'YlOrRd')
seqenf <- seqdef(trajenf, labels=labs, cpal=palette)

# définition de la trajectoire d'indépendance résidentielle
labs <- c("non indépendant","indépendant")
palette <- brewer.pal(3, 'Set1')[1:2]
seqlog <- seqdef(trajlog, labels=labs, cpal=palette)

Association entre dimensions

On calcule la matrice de distance de chacune des dimensions, en utisant l’Optimal Matching (coût de substitution unique et égal à 2, coût indel de 1,5).

# matrices de distances des différentes dimensions
dmat <- seqdist(seqmat, method="OM", indel=1.5, sm=seqsubm(seqmat,"CONSTANT",cval=2))
denf <- seqdist(seqenf, method="OM", indel=1.5, sm=seqsubm(seqenf,"CONSTANT",cval=2))
dlog <- seqdist(seqlog, method="OM", indel=1.5, sm=seqsubm(seqlog,"CONSTANT",cval=2))

On calcule également une matrice de distance à partir d’une analyse de séquences multiples (dite aussi “multichannel sequence analysis”). On décide ici de garder les mêmes coûts pour chacune des dimensions, mais ce n’est pas obligatoire.

# matrice de distances entre séquences multidimensionnelles
dissim.MSA <- seqdistmc(list(seqmat,seqenf,seqlog), method="OM", indel=1.5, sm=as.list(rep("CONSTANT",3)), cval=2)

On étudie ensuite dans quelle mesure les différentes dimensions sont liées, ce qui peut être réalisé de plusieurs manières. D’après les corrélations, on voit que la dimension logement (“log”) est la moins liée à la matrice de distance entre séquences multiples (appelée ici “jsa” pour Joint Sequence Analysis). Elle est de plus très peu associée à la dimension parentale (“enf”). Les dimensions matrimoniale (“mat”) et parentales sont les plus liées. L’examen des alpha de Cronbach confirme ces résultats.

asso <- assoc.domains(list(dmat,denf,dlog), c('mat','enf','log'), dissim.MSA)
asso

$correlations
$correlations$pearson
      mat   enf   log   jsa
mat 1.000 0.372 0.275 0.797
enf 0.372 1.000 0.090 0.702
log 0.275 0.090 1.000 0.604
jsa 0.797 0.702 0.604 1.000

$correlations$spearman
      mat   enf   log   jsa
mat 1.000 0.370 0.280 0.794
enf 0.370 1.000 0.096 0.705
log 0.280 0.096 1.000 0.565
jsa 0.794 0.705 0.565 1.000


$`mean squared correlations`
$`mean squared correlations`$pearson
[1] 0.498

$`mean squared correlations`$spearman
[1] 0.483


$`Cronbach's alpha`
$`Cronbach's alpha`$`(mat,enf,log)`
[1] 0.494

$`Cronbach's alpha`$`(mat,enf)`
[1] 0.542

$`Cronbach's alpha`$`(mat,log)`
[1] 0.431

$`Cronbach's alpha`$`(enf,log)`
[1] 0.166

On peut encore se faire une idée de la structure des associations entre dimensions en réalisant une analyse en composantes principales (ACP) de la matrice de corrélation entre dimensions.

# ACP à partir des corrélations entre dimensions
matcor <- asso$correlations$pearson[1:3,1:3]
PCA <- PCA(matcor, scale.unit=F, graph=F)
plot.PCA(PCA, choix='varcor')

On calcule ensuite plusieurs matrices de distance entre séquences multiples en excluant à chaque l’une des dimensions.

dnolog <- seqdistmc(list(seqmat,seqenf), method="OM", indel=1.5, sm=as.list(rep("CONSTANT",2)), cval=2)
dnomat <- seqdistmc(list(seqenf,seqlog), method="OM", indel=1.5, sm=as.list(rep("CONSTANT",2)), cval=2)
dnoenf <- seqdistmc(list(seqmat,seqlog), method="OM", indel=1.5, sm=as.list(rep("CONSTANT",2)), cval=2)
dnolog <- as.numeric(as.dist(dnolog))
dnomat <- as.numeric(as.dist(dnomat))
dnoenf <- as.numeric(as.dist(dnoenf))

On examine ensuite, pour chaque dimension, la corrélation entre sa matrice de distance et la matrice de distance entre séquences multiples excluant cette dimension. Une fois encore, on voit que c’est la dimension logement qui est la moins associée aux autres.

dmat %>% as.dist %>% as.numeric %>% cor(dnomat) %>% round(3)

[1] 0.418

denf %>% as.dist %>% as.numeric %>% cor(dnoenf) %>% round(3)

[1] 0.304

dlog %>% as.dist %>% as.numeric %>% cor(dnolog) %>% round(3)

[1] 0.216

Une dernière option consiste à trier les index plots de toutes les dimensions selon les résultats d’un multidimensional scaling (MDS) réalisé pour une seule dimension. Ici, on trie les séquences à partir d’un MDS de la dimension matrimoniale. L’ordre des séquences des deux autres dimensions semble obéir à une logique interprétable, on peut donc considérer que les dimensions sont suffisamment liées entre elles pour justifier une analyse de séquences multiples.

mds.msa <- cmdscale(dmat,k=1)
par(mfrow=c(3,2), mar=c(2.1,2.1,2.1,2.1))
seqIplot(seqmat, sortv=mds.msa, xtlab=14:35, with.legend=FALSE, yaxis=FALSE, ylab="")
seqlegend(seqmat, cex=0.7)
seqIplot(seqenf, sortv=mds.msa, xtlab=14:35, with.legend=FALSE, yaxis=FALSE, ylab="")
seqlegend(seqenf, cex=0.7)
seqIplot(seqlog, sortv=mds.msa, xtlab=14:35, with.legend=FALSE, yaxis=FALSE, ylab="")
seqlegend(seqlog, cex=0.7)

Typologie de séquences via Multidimensional Sequence Analysis (MSA)

Pour obtenir une typologie, on réalise une CAH à partir de la matrice de distance issue de l’analyse de séquences multiples.

# classification ascendante hiérarchique
agnes.MSA <- agnes(as.dist(dissim.MSA), method="ward", keep.diss=FALSE)
plot(as.dendrogram(agnes.MSA), leaflab="none")

On opte pour une typologie en 5 classes.

# choix d'une typologie en 5 classes
nbcl.MSA <- 5
part.MSA <- cutree(agnes.MSA, nbcl.MSA) %>% factor

Pour obtenir les chronogrammes des classes de la typologie (un graphique par classe et par dimension) :

# chronogrammes de la typologie
# ngroups <- nlevels(part.MSA)
par(mfrow=c(3,nbcl.MSA+1), mar=c(2.5, 2.1, 2.1, 2.1))
for(i in 1:nbcl.MSA) seqdplot(seqmat[part.MSA==i,], xtlab=14:35, border=NA, with.legend=FALSE, main=paste('classe',i))
seqlegend(seqmat, cex=0.5)
for(i in 1:nbcl.MSA) seqdplot(seqenf[part.MSA==i,], xtlab=14:35, border=NA, with.legend=FALSE)
seqlegend(seqenf, cex=0.5)
for(i in 1:nbcl.MSA) seqdplot(seqlog[part.MSA==i,], xtlab=14:35, border=NA, with.legend=FALSE)
seqlegend(seqlog, cex=0.5)

Typologie de séquences via Globally Interdependent Multiple Sequence Analysis (GIMSA)

Une alternative à l’analyse de séquences multiples est la Globally Interdependent Multiple Sequence Analysis (GIMSA, voir Robette et al, 2015).

On commence par charger les données et les coder sous forme de séquences. Il s’agit de : - 400 trajectoires professionnelles de mères, observées entre 14 et 60 ans avec les états suivants: indépendante, catégorie socioprofessionnelle moyenne/supérieure, catégorie socioprofessionnelle populaire, inactivité, études. - 400 trajectoires d’emploi de leurs filles, observées au cours des 15 premières années après la fin des études, avec les états suivants: études, inactivité, temps partiel , temps plein.

# chargement des données
data(seqgimsa)
trajfilles <- seqgimsa %>% select(starts_with('f'))
str(trajfilles)

'data.frame':   400 obs. of  15 variables:
 $ f1 : num  4 4 4 4 3 4 4 4 4 4 ...
 $ f2 : num  4 4 4 4 4 4 4 4 4 4 ...
 $ f3 : num  4 4 4 4 4 4 4 4 4 2 ...
 $ f4 : num  4 4 4 4 4 4 4 4 4 2 ...
 $ f5 : num  4 4 4 4 4 4 4 2 4 2 ...
 $ f6 : num  4 4 4 4 4 4 4 2 4 2 ...
 $ f7 : num  4 4 4 4 4 4 4 2 4 2 ...
 $ f8 : num  4 4 4 2 4 4 4 2 4 2 ...
 $ f9 : num  4 4 4 4 4 4 4 2 4 2 ...
 $ f10: num  4 4 4 4 4 4 4 2 4 2 ...
 $ f11: num  4 4 4 4 4 4 4 2 4 2 ...
 $ f12: num  4 4 4 4 4 4 4 2 4 2 ...
 $ f13: num  4 4 4 4 4 4 4 2 4 2 ...
 $ f14: num  4 4 4 4 4 4 4 2 4 2 ...
 $ f15: num  4 4 4 4 4 4 4 2 4 2 ...

trajmeres <- seqgimsa %>% select(starts_with('m'))
str(trajmeres)

'data.frame':   400 obs. of  47 variables:
 $ m14: num  5 5 9 9 8 8 9 9 9 9 ...
 $ m15: num  5 5 9 9 8 8 9 9 9 9 ...
 $ m16: num  5 5 9 9 8 8 9 9 9 9 ...
 $ m17: num  5 5 9 9 8 8 9 9 9 9 ...
 $ m18: num  5 5 9 3 8 8 9 9 3 5 ...
 $ m19: num  5 5 9 3 8 8 9 9 3 5 ...
 $ m20: num  5 5 9 3 8 8 9 9 3 5 ...
 $ m21: num  5 5 9 3 8 8 9 5 3 5 ...
 $ m22: num  5 5 9 3 8 8 9 5 3 5 ...
 $ m23: num  5 5 9 3 8 8 9 5 3 5 ...
 $ m24: num  5 5 9 3 8 8 9 5 3 8 ...
 $ m25: num  5 5 9 3 8 8 1 5 3 8 ...
 $ m26: num  5 5 9 3 8 8 1 5 3 8 ...
 $ m27: num  5 5 9 3 8 8 1 5 3 8 ...
 $ m28: num  5 5 9 3 8 8 1 5 3 8 ...
 $ m29: num  5 5 9 3 8 8 1 5 3 8 ...
 $ m30: num  5 5 9 8 8 8 1 5 3 8 ...
 $ m31: num  5 5 9 8 8 8 1 5 3 8 ...
 $ m32: num  5 5 9 8 8 8 1 5 3 8 ...
 $ m33: num  5 5 9 8 8 8 1 5 3 8 ...
 $ m34: num  5 5 8 8 8 8 1 5 3 8 ...
 $ m35: num  5 5 8 8 8 8 1 5 3 8 ...
 $ m36: num  5 5 8 8 8 8 1 5 3 8 ...
 $ m37: num  5 5 8 8 8 8 1 5 3 8 ...
 $ m38: num  5 5 8 8 8 8 1 5 3 8 ...
 $ m39: num  5 5 8 8 8 8 1 5 3 8 ...
 $ m40: num  5 5 8 8 8 8 1 5 3 8 ...
 $ m41: num  5 5 5 8 8 8 1 5 3 8 ...
 $ m42: num  5 5 5 8 8 8 1 5 3 8 ...
 $ m43: num  5 5 5 8 8 8 1 5 3 8 ...
 $ m44: num  5 5 5 8 8 8 1 5 3 8 ...
 $ m45: num  5 5 5 8 8 8 1 5 3 8 ...
 $ m46: num  5 5 5 8 8 8 1 5 3 8 ...
 $ m47: num  5 5 5 8 8 8 1 5 3 8 ...
 $ m48: num  5 5 5 8 8 8 1 5 3 8 ...
 $ m49: num  8 5 5 9 8 8 1 5 3 8 ...
 $ m50: num  8 5 5 9 8 8 1 5 3 8 ...
 $ m51: num  8 5 5 9 8 8 1 5 3 8 ...
 $ m52: num  8 5 5 9 8 8 1 5 3 8 ...
 $ m53: num  8 5 5 9 8 8 1 5 8 8 ...
 $ m54: num  8 5 5 9 8 8 1 5 8 8 ...
 $ m55: num  8 5 5 9 8 8 1 5 8 8 ...
 $ m56: num  8 5 5 9 8 8 1 5 8 8 ...
 $ m57: num  8 5 5 9 8 8 1 5 8 8 ...
 $ m58: num  8 5 5 9 8 8 1 5 8 8 ...
 $ m59: num  8 5 5 9 8 8 1 5 8 8 ...
 $ m60: num  8 8 5 9 8 8 1 8 8 8 ...

# définition des séquences
lab.meres <- c("indép","moyen/sup","popu","inactivité","études")
pal.meres <- brewer.pal(5, "Set1")
seqmeres <- seqdef(trajmeres,lab=lab.meres, cpal=pal.meres)

lab.filles <- c("études","inactivité","temps partiel","temps plein")
pal.meres <- brewer.pal(4, "Set1")
seqfilles <- seqdef(trajfilles,lab=lab.filles, cpla=pal.filles)

La première étape consiste à calculer une matrice de distance pour les mères et une pour les filles. On utilise la métrique LCS pour les mères et la distance de Hamming pour les filles.

# étape 1 : mesure de dissimilarité
dmeres <- seqdist(seqmeres,method="LCS")
cout.filles <- seqsubm(seqfilles, method="CONSTANT", cval=2)
dfilles <- seqdist(seqfilles, method="HAM", sm=cout.filles)

Dans la seconde étape, on résume les matrices de distance à partir d’un MDS.

# étape 2 : multidimensional scaling
mds.meres <- cmdscale(dmeres, k=20, eig=TRUE)
mds.filles <- cmdscale(dfilles, k=20, eig=TRUE)

On sélectionne le nombre de facteurs à retenir pour chacun des MDS.

# choix du nombre de dimensions à retenir pour les mères
par(mfrow=c(1,2))
# mesure de stress
seqmds.stress(dmeres, mds.meres) %>% plot(type='l', xlab='nombre de facteurs', ylab='stress')
# part de variance expliquée
(mds.meres$eig[1:10]/mds.meres$eig[1]) %>% plot(type='s', xlab='nombre de facteurs', ylab='part de variance expliquée')

# choix du nombre de dimensions à retenir pour les filles
par(mfrow=c(1,2))
# mesure de stress
seqmds.stress(dfilles, mds.filles) %>% plot(type='l', xlab='nombre de facteurs', ylab='stress')
# part de variance expliquée
(mds.filles$eig[1:10]/mds.filles$eig[1]) %>% plot(type='s', xlab='nombre de facteurs', ylab='part de variance expliquée')

Dans la troisième étape, on résume les relations entre les résultats du MDS des mères et ceux du MDS des filles, à l’aide d’une PLS symétrique.

# étape 3 : PLS symétrique
a <- mds.meres$points[,1:5]
b <- mds.filles$points[,1:4]
pls <- symPLS(a,b)

Dans la quatrième et dernière étape, on calcule une matrice de distance unique entre les dyades de trajectoires mères-filles. Il est possible de pondérer les dimensions des mères et des filles, pour équilibrer leur contribution aux résultats finaux.

# étape 4 : distance et classification

# pas de pondération
F <- pls$F
G <- pls$G

# pondération par la variance des composantes de la PLS (w1)
F <- apply(pls$F,2,scale,center=FALSE)
G <- apply(pls$G,2,scale,center=FALSE)

# pondération par le nombre de séquences distinctes (w2)
F <- pls$F/nrow(seqtab(seqmeres,tlim=0))
G <- pls$G/nrow(seqtab(seqfilles,tlim=0))

# pondération par la 1ère valeur propre du MDS (w3)
F <- pls$F/mds.meres$eig[1]
G <- pls$G/mds.filles$eig[1]

On utilise ici la pondération w1.

# pondération par la variance des composantes de la PLS (w1)
F <- apply(pls$F,2,scale,center=FALSE)
G <- apply(pls$G,2,scale,center=FALSE)

# calcul de distance
diff2 <- function(X) return(as.matrix(dist(X,upper=T,diag=T)^2,nrow=nrow(X)))
D <- (diff2(F)+diff2(G))^0.5

On réalise ensuite une classification automatique (ici une CAH).

# classification
dist.GIMSA <- as.dist(D)
agnes.GIMSA <- agnes(dist.GIMSA, method="ward", keep.diss=FALSE)
plot(as.dendrogram(agnes.GIMSA), leaflab="none")

On choisit une partition en 5 classes.

# partition en 5 classes
nbcl.GIMSA <- 5
part.GIMSA <- cutree(agnes.GIMSA, nbcl.GIMSA) %>% factor

Pour finir, on représente graphiquement la typologie à l’aide de chronogrammes.

par(mfrow=c(3,nbcl.GIMSA), mar=c(2.5, 2.1, 2.1, 2.1))
for(i in 1:nbcl.GIMSA) seqdplot(seqmeres[part.GIMSA==i,], xtlab=14:60, border=NA, with.legend=FALSE, main=paste('classe',i))
for(i in 1:nbcl.GIMSA) seqdplot(seqfilles[part.GIMSA==i,], xtlab=1:15, border=NA, with.legend=FALSE)
seqlegend(seqmeres, cex=0.6)
seqlegend(seqfilles, cex=0.6)