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")