Ce chapitre est évoqué dans le webin-R #13 (exemples de graphiques avancés) sur YouTube.
+
Tidy data
Comme indiqué dans l’introduction au tidyverse, les extensions du tidyverse comme dplyr ou ggplot2 partent du principe que les données sont “bien rangées” sous forme de tidy data.
@@ -10763,7 +10766,7 @@
Fichiers volumineux
-
+
@@ -10960,7 +10963,7 @@
Combinaison des résultats
-
+
@@ -11089,7 +11092,7 @@
Sauvegarder des objets
-
+
@@ -11194,7 +11197,7 @@
Export avec ggplot2
-
+
@@ -12602,7 +12605,7 @@
Exporter les graphiques obtenus
-
+
@@ -14998,7 +15001,7 @@
Tableaux faciles avec gtsummary
-
+
@@ -15044,6 +15047,9 @@
Introduction à ggplot2, la grammaire des graphique
Ce chapitre est tiré d’une séance de cours de François Briatte et destinée à des étudiants de L2 sans aucune connaissance de R. Cette séance de cours est elle-même inspirée d’un exercice tiré d’un cours de Cosma Shalizi.
+
+
Ce chapitre est évoqué dans le webin-R #08 (ggplot2 et la grammaires des graphiques) sur YouTube.
+
R possède un puissant moteur graphique interne, qui permet de dessiner dans un graphique en y rajoutant des segments, des points, du texte, ou toutes sortes d’autres symboles. Toutefois, pour produire un graphique complet avec les fonctions basiques de R, il faut un peu bricoler : d’abord, ouvrir une fenêtre ; puis rajouter des points ; puis rajouter des lignes ; tout en configurant les couleurs au fur-et-à-mesure ; puis finir par fermer la fenêtre graphique.
L’extension ggplot21, développée par Hadley Wickham et mettant en œuvre la grammaire graphiquethéorisée par Leland Wilkinson, devient vite indispensable lorsque l’on souhaite réaliser des graphiques plus complexes2. On renvoit le lecteur intéressé à l’ouvrage de Winston Chang, R Graphics Cookbook, disponible gratuitement dans sa deuxième édition à l’adresse suivante : https://r-graphics.org.
Ce chapitre, articulé autour d’une étude de cas, présente ggplot2 à partir d’un exemple simple de visualisation de séries temporelles, puis rentre dans le détail de sa syntaxe. Pour une présentation plus formelle, on pourra se référer au chapitre dédié de la section Approfondir.
@@ -15286,7 +15292,7 @@
Extensions de ggplot2
-
+
@@ -15315,6 +15321,7 @@
Graphiques univariés et bivariés avec ggplot2
Ce chapitre est évoqué dans le webin-R #03 (statistiques descriptives avec gtsummary et esquisse) sur YouTube.
+
Ce chapitre est évoqué dans le webin-R #09 (Graphiques uni- et bivariés avec ggplot2) sur YouTube.
Après avoir introduit l’extension ggplot2 au travers d’une étude de cas, nous reprenons ici les graphiques produits dans les chapitres statistique univariée et statistique bivariée et montrons comment les réaliser avec ggplot2.
@@ -15835,7 +15842,7 @@
Exporter les graphiques obtenus
-
+
@@ -15861,6 +15868,9 @@
Données pondérées
+
+
Ce chapitre est évoqué dans le webin-R #10 (Données pondérées, plan d’échantillonnage complexe & survey) sur YouTube.
+
S’il est tout à fait possible de travailler avec des données pondérées sous R, cette fonctionnalité n’est pas aussi bien intégrée que dans la plupart des autres logiciels de traitement statistique. En particulier, il y a plusieurs manières possibles de gérer la pondération. Cependant, lorsque l’on doit également prendre un compte un plan d’échantillonnage complexe (voir section dédiée ci-après), R fournit tous les outils nécessaires, alors que dans la plupart des logiciels propriétaires, il faut disposer d’une extension adéquate, pas toujours vendue de base avec le logiciel.
Dans ce qui suit, on utilisera le jeu de données tiré de l’enquête Histoire de vie et notamment sa variable de pondération poids1.
library(questionr)
@@ -16887,7 +16897,7 @@
Conclusion
-
+
@@ -17066,7 +17076,7 @@
Données pondérées et l’extension survey
-
+
@@ -18548,7 +18558,7 @@
Tests dans les tableaux de gtsummary
-
+
@@ -18565,6 +18575,9 @@
Définir un plan d’échantillonnage complexe
+
+
Ce chapitre est évoqué dans le webin-R #10 (Données pondérées, plan d’échantillonnage complexe & survey) sur YouTube.
+
L’extension survey ne permet pas seulement d’indiquer une variable de pondération mais également de prendre les spécificités du plan d’échantillonnage (strates, grappes, …). Le plan d’échantillonnage ne joue pas seulement sur la pondération des données, mais influence le calcul des variances et par ricochet tous les tests statistiques. Deux échantillons identiques avec la même variable de pondération mais des designs différents produiront les mêmes moyennes et proportions mais des intervalles de confiance différents.
On notera la présence d’une colonne supplémentaire, y.level. De fait, la fonction ggcoef_model de GGally ne peut s’appliquer directement, car les coefficients vont se supperposer. On aura dès lors recours à sa variante ggcoef_multinom dédiée aux modèles multinomiax.
-
-Représentation alternative des effets marginaux du modèle multinomial
-
-
-
-
Régression logistique ordinale
-
La régression logistique ordinale s’applique lorsque la variable à expliquer possède trois ou plus modalités qui sont ordonnées (par exemple : modéré, moyen, fort).
-
L’extension la plus utilisée pour réaliser des modèles ordinaux est ordinal et sa fonction clm. Il est même possible de réaliser des modèles ordinaux avec des effets aléatoires (modèles mixtes) à l’aide de la fonction clmm.
On va reprendre l’exemple précédent puisque la variable trav.satisf est une variable ordonnée.
-
freq(d$trav.satisf)
-
-
-
-
ATTENTION : Dans le cas d’une régression logistique ordinale, il importante que les niveaux du facteur soient classés selon leur ordre hiéarchique (du plus faible au plus fort). On va dès lors recoder notre variable à expliquer.
i Multinomial models have a different underlying structure than the models
+gtsummary was designed for. Other gtsummary functions designed to work with
+tbl_regression objects may yield unexpected results.
+
Warning: The `.dots` argument of `group_by()` is deprecated
+as of dplyr 1.0.0.
-
+
Characteristic
-
exp(Beta)
+
+OR1
+
95% CI1
p-value
+
+
Satisfaction
+
-
Sexe
+
Niveau d'étude
-
Femme
+
Primaire
—
—
-
Homme
-
0.85
-
0.67, 1.08
-
0.2
+
Secondaire
+
1.05
+
0.63, 1.76
+
0.9
+
+
+
Technique/Professionnel
+
1.08
+
0.67, 1.73
+
0.7
+
+
+
Supérieur
+
2.01
+
1.24, 3.27
+
0.005
+
+
+
manquant
+
0.58
+
0.18, 1.86
+
0.4
+
+
+
Importance du travail
+
+
+
+
+
+
Le plus
+
—
+
—
+
+
+
+
Aussi
+
1.29
+
0.56, 2.98
+
0.5
+
+
+
Moins
+
0.84
+
0.37, 1.88
+
0.7
+
+
+
Peu
+
0.55
+
0.18, 1.64
+
0.3
+
+
+
Insatisfaction
Niveau d'étude
@@ -23986,27 +23926,27 @@
Régression logistique ordinale
Secondaire
-
1.06
-
0.67, 1.67
+
0.91
+
0.41, 1.99
0.8
Technique/Professionnel
-
1.03
-
0.68, 1.57
-
0.9
+
1.09
+
0.54, 2.19
+
0.8
Supérieur
-
1.85
-
1.20, 2.84
-
0.005
+
1.08
+
0.51, 2.29
+
0.8
manquant
-
0.63
-
0.24, 1.66
-
0.4
+
0.96
+
0.18, 4.97
+
>0.9
Importance du travail
@@ -24022,21 +23962,21 @@
Régression logistique ordinale
Aussi
-
1.42
-
0.66, 3.03
-
0.4
+
0.80
+
0.24, 2.69
+
0.7
Moins
-
1.03
-
0.49, 2.13
-
>0.9
+
0.59
+
0.18, 1.88
+
0.4
Peu
-
0.19
-
0.08, 0.47
-
<0.001
+
3.82
+
1.05, 13.9
+
0.042
@@ -24045,75 +23985,658 @@
Régression logistique ordinale
1
- CI = Confidence Interval
+ OR = Odds Ratio, CI = Confidence Interval
+
car::Anova(regm2)
+
+
-
-
Données pondérées et l’extension survey
-
Lorsque l’on utilise des données pondérées, on aura recours à l’extension survey6.
-
Préparons des données d’exemple :
-
library(survey)
-dw <-svydesign(ids =~1, data = d, weights =~poids)
-
-
Régression logistique binaire
-
L’extension survey fournit une fonction svyglm permettant de calculer un modèle statistique tout en prenant en compte le plan d’échantillonnage spécifié. La syntaxe de svyglm est proche de celle de glm. Cependant, le cadre d’une régression logistique, il est nécessaire d’utiliser family = quasibinomial() afin d’éviter un message d’erreur indiquant un nombre non entier de succès :
-
reg <-svyglm(sport ~ sexe + age + relig + heures.tv, dw, family =binomial())
-
Warning in eval(family$initialize): non-integer #successes
-in a binomial glm!
-
reg <-svyglm(sport ~ sexe + age + relig + heures.tv, dw, family =quasibinomial())
-reg
-
Independent Sampling design (with replacement)
-svydesign(ids = ~1, data = d, weights = ~poids)
-
-Call: svyglm(formula = sport ~ sexe + age + relig + heures.tv, design = dw,
- family = quasibinomial())
-
-Coefficients:
- (Intercept)
- 1.53590
- sexeHomme
- 0.36526
- age
- -0.04127
- religPratiquant occasionnel
- 0.05577
- religAppartenance sans pratique
- 0.16367
-religNi croyance ni appartenance
- 0.03988
- religRejet
- -0.14862
- religNSP ou NVPR
- -0.22682
- heures.tv
- -0.18204
-
-Degrees of Freedom: 1994 Total (i.e. Null); 1986 Residual
- (5 observations deleted due to missingness)
-Null Deviance: 2672
-Residual Deviance: 2378 AIC: NA
-
Le résultat obtenu est similaire à celui de glm et l’on peut utiliser sans problème les fonctions coef, confint, odds.ratio, predict ou encore tidy, tidy_plus_plus et ggcoef_model.
-
odds.ratio(reg)
+
De même, il est possible de calculer la matrice de confusion :
On notera la présence d’une colonne supplémentaire, y.level. De fait, la fonction ggcoef_model de GGally ne peut s’appliquer directement, car les coefficients vont se supperposer. On aura dès lors recours à sa variante ggcoef_multinom dédiée aux modèles multinomiax.
+
+Représentation alternative des effets marginaux du modèle multinomial
+
+
+
+
Régression logistique ordinale
+
La régression logistique ordinale s’applique lorsque la variable à expliquer possède trois ou plus modalités qui sont ordonnées (par exemple : modéré, moyen, fort).
+
L’extension la plus utilisée pour réaliser des modèles ordinaux est ordinal et sa fonction clm. Il est même possible de réaliser des modèles ordinaux avec des effets aléatoires (modèles mixtes) à l’aide de la fonction clmm.
On va reprendre l’exemple précédent puisque la variable trav.satisf est une variable ordonnée.
+
freq(d$trav.satisf)
+
+
ATTENTION : Dans le cas d’une régression logistique ordinale, il importante que les niveaux du facteur soient classés selon leur ordre hiéarchique (du plus faible au plus fort). On va dès lors recoder notre variable à expliquer.
Lorsque l’on utilise des données pondérées, on aura recours à l’extension survey6.
+
Préparons des données d’exemple :
+
library(survey)
+dw <-svydesign(ids =~1, data = d, weights =~poids)
+
+
Régression logistique binaire
+
L’extension survey fournit une fonction svyglm permettant de calculer un modèle statistique tout en prenant en compte le plan d’échantillonnage spécifié. La syntaxe de svyglm est proche de celle de glm. Cependant, le cadre d’une régression logistique, il est nécessaire d’utiliser family = quasibinomial() afin d’éviter un message d’erreur indiquant un nombre non entier de succès :
+
reg <-svyglm(sport ~ sexe + age + relig + heures.tv, dw, family =binomial())
+
Warning in eval(family$initialize): non-integer #successes
+in a binomial glm!
+
reg <-svyglm(sport ~ sexe + age + relig + heures.tv, dw, family =quasibinomial())
+reg
+
Independent Sampling design (with replacement)
+svydesign(ids = ~1, data = d, weights = ~poids)
+
+Call: svyglm(formula = sport ~ sexe + age + relig + heures.tv, design = dw,
+ family = quasibinomial())
+
+Coefficients:
+ (Intercept)
+ 1.53590
+ sexeHomme
+ 0.36526
+ age
+ -0.04127
+ religPratiquant occasionnel
+ 0.05577
+ religAppartenance sans pratique
+ 0.16367
+religNi croyance ni appartenance
+ 0.03988
+ religRejet
+ -0.14862
+ religNSP ou NVPR
+ -0.22682
+ heures.tv
+ -0.18204
+
+Degrees of Freedom: 1994 Total (i.e. Null); 1986 Residual
+ (5 observations deleted due to missingness)
+Null Deviance: 2672
+Residual Deviance: 2378 AIC: NA
+
Le résultat obtenu est similaire à celui de glm et l’on peut utiliser sans problème les fonctions coef, confint, odds.ratio, predict ou encore tidy, tidy_plus_plus et ggcoef_model.
+
odds.ratio(reg)
+
+
+
+
tidy_plus_plus(reg, exponentiate =TRUE)
+
+
+
+
tbl_regression(reg, exponentiate =TRUE)
+
-
+
-
Caractéristique
+
Characteristic
-OR1
+1, N = 1,0311
-95% CI1
+2, N = 1641
+
+
+3, N = 5531
+
+
+4, N = 2051
+
+
+5, N = 471
-
p-value
-
Sexe
+
grpage
+
+
-
Femme
-
—
-
—
+
[16,25)
+
0 (0%)
+
164 (100%)
+
0 (0%)
+
0 (0%)
+
5 (11%)
+
+
+
[25,45)
+
590 (57%)
+
0 (0%)
+
20 (3.6%)
+
80 (39%)
+
16 (34%)
+
+
+
[45,65)
+
441 (43%)
+
0 (0%)
+
185 (34%)
+
97 (47%)
+
22 (47%)
+
+
+
[65,93]
+
0 (0%)
+
0 (0%)
+
346 (63%)
+
28 (14%)
+
4 (8.5%)
+
+
+
Unknown
+
0
+
0
+
2
+
0
+
0
+
+
+
sexe
+
+
+
+
Homme
-
1,55
-
1,26 – 1,91
-
<0,001
+
429 (42%)
+
74 (45%)
+
220 (40%)
+
161 (79%)
+
15 (32%)
-
Groupe d'âges
+
Femme
+
602 (58%)
+
90 (55%)
+
333 (60%)
+
44 (21%)
+
32 (68%)
+
+
+
etud
+
+
-
[16,25)
-
—
-
—
-
+
Primaire
+
7 (0.7%)
+
1 (1.3%)
+
402 (73%)
+
50 (25%)
+
6 (14%)
-
[25,45)
-
0,66
-
0,42 – 1,03
-
0,065
+
Secondaire
+
267 (26%)
+
15 (19%)
+
65 (12%)
+
35 (17%)
+
5 (12%)
-
[45,65)
-
0,34
-
0,21 – 0,54
-
<0,001
+
Technique/Professionnel
+
406 (40%)
+
48 (62%)
+
48 (8.7%)
+
87 (43%)
+
5 (12%)
-
[65,99]
-
0,25
-
0,15 – 0,43
-
<0,001
+
Supérieur
+
334 (33%)
+
14 (18%)
+
36 (6.5%)
+
30 (15%)
+
27 (63%)
-
Niveau d'étude
+
Unknown
+
17
+
86
+
2
+
3
+
4
+
+
+
peche.chasse
+
+
-
Primaire
-
—
-
—
+
Non
+
1,030 (100%)
+
151 (92%)
+
553 (100%)
+
0 (0%)
+
42 (89%)
+
+
+
Oui
+
1 (<0.1%)
+
13 (7.9%)
+
0 (0%)
+
205 (100%)
+
5 (11%)
+
+
+
cinema
+
+
+
+
-
Secondaire
-
2,59
-
1,77 – 3,83
-
<0,001
+
Non
+
497 (48%)
+
32 (20%)
+
486 (88%)
+
141 (69%)
+
18 (38%)
-
Technique/Professionnel
-
2,86
-
1,98 – 4,17
-
<0,001
+
Oui
+
534 (52%)
+
132 (80%)
+
67 (12%)
+
64 (31%)
+
29 (62%)
-
Supérieur
-
6,63
-
4,55 – 9,80
-
<0,001
+
cuisine
+
+
+
+
+
-
Manquant
-
8,59
-
4,53 – 16,6
-
<0,001
+
Non
+
539 (52%)
+
85 (52%)
+
361 (65%)
+
111 (54%)
+
23 (49%)
-
Nombre d'heures passées devant la télévision par jour
-
0,89
-
0,83 – 0,95
-
<0,001
+
Oui
+
492 (48%)
+
79 (48%)
+
192 (35%)
+
94 (46%)
+
24 (51%)
-
Pratique religieuse
+
bricol
+
+
-
Pratiquant regulier
-
—
-
—
+
Non
+
537 (52%)
+
103 (63%)
+
397 (72%)
+
89 (43%)
+
21 (45%)
+
+
+
Oui
+
494 (48%)
+
61 (37%)
+
156 (28%)
+
116 (57%)
+
26 (55%)
+
+
+
sport
+
+
+
+
-
Pratiquant occasionnel
-
0,98
-
0,68 – 1,42
-
>0,9
+
Non
+
588 (57%)
+
57 (35%)
+
479 (87%)
+
130 (63%)
+
23 (49%)
-
Appartenance sans pratique
-
0,99
-
0,71 – 1,40
-
>0,9
+
Oui
+
443 (43%)
+
107 (65%)
+
74 (13%)
+
75 (37%)
+
24 (51%)
-
Ni croyance ni appartenance
-
0,81
-
0,55 – 1,18
-
0,3
+
lecture.bd
+
+
+
+
+
-
Rejet
-
0,68
-
0,39 – 1,19
-
0,2
+
Non
+
1,031 (100%)
+
164 (100%)
+
553 (100%)
+
205 (100%)
+
0 (0%)
-
NSP ou NVPR
-
0,92
-
0,40 – 2,02
-
0,8
+
Oui
+
0 (0%)
+
0 (0%)
+
0 (0%)
+
0 (0%)
+
47 (100%)
-
+
1
- OR = rapport de cotes, CI = intervalle de confiance
+ n (%)
-
Selon les résultats de notre modèle, les hommes pratiquent plus un sport que les femmes et la pratique du sport diminue avec l’âge. Pour représenter les effets différentes variables, on peut avoir recours à la fonction allEffects de l’extension effects.
-
library(effects)
-
Loading required package: carData
-
Registered S3 methods overwritten by 'lme4':
- method from
- cooks.distance.influence.merMod car
- influence.merMod car
- dfbeta.influence.merMod car
- dfbetas.influence.merMod car
-
lattice theme set by effectsTheme()
-See ?effectsTheme for details.
-
plot(allEffects(mod))
-
-Représentation graphique des effets du modèle
-
Cependant, l’effet de l’âge est-il le même selon le sexe ? Nous allons donc introduire une interaction entre l’âge et le sexe dans notre modèle, ce qui sera représenté par sexe * grpage dans l’équation du modèle.
-
mod2 <-glm(sport ~ sexe * grpage + etud + heures.tv + relig, data = d, family =binomial())
-tbl_regression(mod2, exponentiate =TRUE)
+
On peut également avoir recours à ggtable de GGally pour représenter les résidus du Chi² et mieux repérer les différences.
Warning in chisq.test(xtabs(weight ~ y + x, data = data)):
+Chi-squared approximation may be incorrect
+
+Warning in chisq.test(xtabs(weight ~ y + x, data = data)):
+Chi-squared approximation may be incorrect
+
+Warning in chisq.test(xtabs(weight ~ y + x, data = data)):
+Chi-squared approximation may be incorrect
+
+Warning in chisq.test(xtabs(weight ~ y + x, data = data)):
+Chi-squared approximation may be incorrect
+
+Warning in chisq.test(xtabs(weight ~ y + x, data = data)):
+Chi-squared approximation may be incorrect
+
+Warning in chisq.test(xtabs(weight ~ y + x, data = data)):
+Chi-squared approximation may be incorrect
+
+Résidus du Chi²
+
+
+
+
CAH avec l’extension FactoMineR
+
L’extension FactoMineR fournit une fonction HCPC permettant de réaliser une classification hiérarchique à partir du résultats d’une analyse factorielle réalisée avec la même extension (voir la section dédiée du chapitre sur l’ACM).
+
HCPC réalise à la fois le calcul de la matrice des distances, du dendrogramme et le partitionnement de la population en classes. Par défaut, HCPC calcule le dendrogramme à partir du carré des distances du F² et avec la méthode de Ward.
+
Par défaut, l’arbre est affiché à l’écran et l’arbre sera coupé selon la partition ayant la plus grande perte relative d’inertie (comme avec best.cutree). Utilisez graph=FALSE pour ne pas afficher le graphique et l’argument nb.clust pour indiquer le nombre de classes désirées.
On pourra représenter le dendrogramme avec plot et l’argument choice="tree".
+
plot(cah, choice ="tree")
+
+Dendrogramme obtenu avec HCPC (5 axes)
+
Il apparait que le dendrogramme obtenu avec HCPC diffère de celui que nous avons calculé précédemment en utilisant la matrice des distances fournies par dist.dudi. Cela est dû au fait que HCPC procède différement pour calculer la matrice des distances en ne prenant en compte que les axes retenus dans le cadre de l’ACM. Pour rappel, nous avions retenu que 5 axes dans le cadre de notre ACM :
HCPC n’a donc pris en compte que ces 5 premiers axes pour calculer les distances entre les individus, considérant que les autres axes n’apportent que du « bruit » rendant la classification instable. Cependant, comme le montre summary(acm2), nos cinq premiers axes n’expliquent que 54 % de la variance. Il usuellement préférable de garder un plus grande nombre d’axes afin de couvrir au moins 80 à 90 % de la variance10. De son côté, dist.dudi prends en compte l’ensemble des axes pour calculer la matrice des distances. On peut reproduire cela avec FactoMineR en indiquant ncp=Inf lors du calcul de l’ACM.
+
+Dendrogramme obtenu avec HCPC (tous les axes)
+
D’autres graphiques sont disponibles, en faisant varier la valeur de l’argument choice :
+
plot(cah, choice ="3D.map")
+
+Représentation en 3 dimensions du dendrogramme
+
plot(cah, choice ="bar")
+
+Gains d’inertie
+
plot(cah, choice ="map")
+
+Projection des catégories sur le plan factoriel
+
L’objet renvoyé par HCPC contient de nombreuses informations. La partition peut notamment être récupérée avec cah$data.clust$clust. Il y a également diverses statistiques pour décrire les catégories.
+
cah
+
**Results for the Hierarchical Clustering on Principal Components**
+ name
+1 "$data.clust"
+2 "$desc.var"
+3 "$desc.var$test.chi2"
+4 "$desc.axes$category"
+5 "$desc.axes"
+6 "$desc.axes$quanti.var"
+7 "$desc.axes$quanti"
+8 "$desc.ind"
+9 "$desc.ind$para"
+10 "$desc.ind$dist"
+11 "$call"
+12 "$call$t"
+ description
+1 "dataset with the cluster of the individuals"
+2 "description of the clusters by the variables"
+3 "description of the cluster var. by the categorical var."
+4 "description of the clusters by the categories."
+5 "description of the clusters by the dimensions"
+6 "description of the cluster var. by the axes"
+7 "description of the clusters by the axes"
+8 "description of the clusters by the individuals"
+9 "parangons of each clusters"
+10 "specific individuals"
+11 "summary statistics"
+12 "description of the tree"
Cette même fonction peut aussi être utilisée pour calculer une distance après une analyse en composantes principales ou une analyse mixte de Hill et Smith.
Pour une description mathématique plus détaillée de cette fonction, notamment en cas de valeur manquante, se référer à l’article original de Gower précédemment cité.
+
On pourra consulter le cours de FG Carpentier déjà cité ou bien des ouvrages d’analyse statistique.
Depuis la version 3.1 de R. L’option method = "ward.D" correspondant à la version disponible dans les versions précédentes de R. Mais il est à noter que la méthode décrite par Ward dans son article de 1963 correspond en réalité à method = "ward.D2.
Ce chapitre est évoqué dans le webin-R #07 (régression logistique partie 2) sur YouTube.
+
+
Dans un modèle statistique classique, on fait l’hypothèse implicite que chaque variable explicative est indépendante des autres. Cependant, cela ne se vérifie pas toujours. Par exemple, l’effet de l’âge peut varier en fonction du sexe. Il est dès lors nécessaire de prendre en compte dans son modèle les effets d’interaction1.
+
+
Exemple d’interaction
+
Reprenons le modèle que nous avons utilisé dans le chapitre sur la régression logistique.
Nous avions alors exploré les facteurs associés au fait de pratiquer du sport.
+
mod <-glm(sport ~ sexe + grpage + etud + heures.tv + relig, data = d, family =binomial())
+library(gtsummary)
+tbl_regression(mod, exponentiate =TRUE)
-
+
Caractéristique
@@ -26421,8 +27233,8 @@
Exemple d’interaction
Homme
-
5,20
-
2,46 – 11,7
+
1,55
+
1,26 – 1,91
<0,001
@@ -26439,21 +27251,21 @@
Exemple d’interaction
[25,45)
-
1,02
-
0,59 – 1,80
-
>0,9
+
0,66
+
0,42 – 1,03
+
0,065
[45,65)
-
0,61
-
0,34 – 1,10
-
0,10
+
0,34
+
0,21 – 0,54
+
<0,001
[65,99]
-
0,46
-
0,23 – 0,92
-
0,027
+
0,25
+
0,15 – 0,43
+
<0,001
Niveau d'étude
@@ -26469,26 +27281,26 @@
Exemple d’interaction
Secondaire
-
2,60
-
1,77 – 3,85
+
2,59
+
1,77 – 3,83
<0,001
Technique/Professionnel
-
2,88
-
1,99 – 4,21
+
2,86
+
1,98 – 4,17
<0,001
Supérieur
-
6,77
-
4,64 – 10,0
+
6,63
+
4,55 – 9,80
<0,001
Manquant
-
9,02
-
4,68 – 17,8
+
8,59
+
4,53 – 16,6
<0,001
@@ -26511,57 +27323,33 @@
Exemple d’interaction
Pratiquant occasionnel
-
0,99
-
0,68 – 1,43
+
0,98
+
0,68 – 1,42
>0,9
Appartenance sans pratique
-
1,02
-
0,73 – 1,44
-
0,9
+
0,99
+
0,71 – 1,40
+
>0,9
Ni croyance ni appartenance
-
0,83
-
0,57 – 1,22
+
0,81
+
0,55 – 1,18
0,3
Rejet
0,68
-
0,38 – 1,19
+
0,39 – 1,19
0,2
NSP ou NVPR
-
0,93
-
0,40 – 2,04
-
0,9
-
-
-
Sexe * Groupe d'âges
-
-
-
-
-
-
Homme * [25,45)
-
0,31
-
0,13 – 0,70
-
0,006
-
-
-
Homme * [45,65)
-
0,24
-
0,10 – 0,54
-
<0,001
-
-
-
Homme * [65,99]
-
0,23
-
0,09 – 0,59
-
0,003
+
0,92
+
0,40 – 2,02
+
0,8
@@ -26575,74 +27363,28 @@
Exemple d’interaction
-
Commençons par regarder les effets du modèle.
-
plot(allEffects(mod2))
-
-Représentation graphique des effets du modèle avec interaction entre le sexe et le groupe d’âge
-
Sur ce graphique, on voit que l’effet de l’âge sur la pratique d’un sport est surtout marqué chez les hommes. Chez les femmes, le même effet est observé, mais dans une moindre mesure et seulement à partir de 45 ans.
-
On peut tester si l’ajout de l’interaction améliore significativement le modèle avec anova.
-
anova(mod2, test ="Chisq")
-
-
-
-
Jetons maintenant un oeil aux coefficients du modèle. Pour rendre les choses plus visuelles, nous aurons recours à ggcoef_model de l’extension GGally.
-
-Représentation graphique des coefficients du modèle avec interaction entre le sexe et le groupe d’âge
-
Concernant l’âge et le sexe, nous avons trois séries de coefficients : trois coefficients (grpage[25,45), grpage[45,65) et grpage[65,99]) qui correspondent à l’effet global de la variable âge, un coefficient (sexeHomme)pour l’effet global du sexe et trois coefficients qui sont des moficateurs de l’effet d’âge pour les hommes (grpage[25,45), grpage[45,65) et grpage[65,99]).
-
Pour bien interpréter ces coefficients, il faut toujours avoir en tête les modalités choisies comme référence pour chaque variable. Supposons une femme de 60 ans, dont toutes lautres variables correspondent aux modalités de référence (c’est donc une pratiquante régulière, de niveau primaire, qui ne regarde pas la télévision). Regardons ce que prédit le modèle quant à sa probabilité de faire du sport au travers d’une représentation graphique
Scale for 'y' is already present. Adding another scale
-for 'y', which will replace the existing scale.
-
-Représentation graphique de l’estimation de la probabilité de faire du sport pour une femme de 60 ans
-
En premier lieu, l’intercept s’applique et permet de déterminer la probabilité de base de faire du sport (si toutes les variables sont à leur valeur de référence). Femme étant la modalité de référence pour la variable sexe, cela ne modifie pas le calcul de la probabilité de faire du sport. Par contre, il y a une modification induite par la modalité 45-65 de la variable grpage.
-
Regardons maintenant la situation d’un homme de 20 ans.
-
nouvelle_observation$sexe[1] <-"Homme"
-nouvelle_observation$grpage[1] <-"[16,25)"
-plot(
-broken(mod2, nouvelle_observation, predict.function = betas),
-trans = logit
-) +ylim(0, 1) +ylab("Probabilité de faire du sport")
-
Scale for 'y' is already present. Adding another scale
-for 'y', which will replace the existing scale.
-
-Représentation graphique de l’estimation de la probabilité de faire du sport pour un homme de 20 ans
-
Nous sommes à la modalité de référence pour l’âge par contre il y a un effet important du sexe. Le coefficient associé globalement à la variable sexe correspond donc à l’effet du sexe à la modalité de référence du groupe d’âges.
-
La situation est différente pour un homme de 60 ans.
-
nouvelle_observation$grpage[1] <-"[45,65)"
-plot(
-broken(mod2, nouvelle_observation, predict.function = betas),
-trans = logit
-) +ylim(0, 1) +ylab("Probabilité de faire du sport")
-
Scale for 'y' is already present. Adding another scale
-for 'y', which will replace the existing scale.
-
-Représentation graphique de l’estimation de la probabilité de faire du sport pour un homme de 60 ans
-
Cette fois-ci, il y a plusieurs modifications d’effet. On applique en effet à la fois le coefficient sexe = Homme (effet du sexe pour les 15-24 ans), le coefficient grpage = [45-65) qui est l’effet de l’âge pour les femmes de 45-64 ans et le coefficient sexe:grpage = Homme:[45-65) qui indique l’effet spécifique qui s’applique aux hommes de 45-64, d’une part par rapport aux femmes du même et d’autre part par rapport aux hommes de 16-24 ans. L’effet des coefficients d’interaction doivent donc être interprétés par rapport aux autres coefficients du modèle qui s’appliquent, en tenant compte des modalités de référence.
-
Il est cependant possible d’écrire le même modèle différemment. En effet, sexe * grpage dans la formule du modèle est équivalent à l’écriture sexe + grpage + sexe:grpage, c’est-à-dire à modéliser un coefficient global pour chaque variable plus un des coefficients d’interaction. On aurait pu demander juste des coefficients d’interaction, en ne mettant que sexe:grpage.
-
mod3 <-glm(sport ~ sexe:grpage + etud + heures.tv + relig, data = d, family =binomial())
-tbl_regression(mod3, exponentiate =TRUE)
+
Selon les résultats de notre modèle, les hommes pratiquent plus un sport que les femmes et la pratique du sport diminue avec l’âge. Pour représenter les effets différentes variables, on peut avoir recours à la fonction allEffects de l’extension effects.
+
library(effects)
+
Loading required package: carData
+
Registered S3 methods overwritten by 'lme4':
+ method from
+ cooks.distance.influence.merMod car
+ influence.merMod car
+ dfbeta.influence.merMod car
+ dfbetas.influence.merMod car
+
lattice theme set by effectsTheme()
+See ?effectsTheme for details.
+
plot(allEffects(mod))
+
+Représentation graphique des effets du modèle
+
Cependant, l’effet de l’âge est-il le même selon le sexe ? Nous allons donc introduire une interaction entre l’âge et le sexe dans notre modèle, ce qui sera représenté par sexe * grpage dans l’équation du modèle.
+
mod2 <-glm(sport ~ sexe * grpage + etud + heures.tv + relig, data = d, family =binomial())
+tbl_regression(mod2, exponentiate =TRUE)
-
+
Caractéristique
@@ -26987,6 +27729,54 @@
Exemple d’interaction
+
Sexe
+
+
+
+
+
+
Femme
+
—
+
—
+
+
+
+
Homme
+
5,20
+
2,46 – 11,7
+
<0,001
+
+
+
Groupe d'âges
+
+
+
+
+
+
[16,25)
+
—
+
—
+
+
+
+
[25,45)
+
1,02
+
0,59 – 1,80
+
>0,9
+
+
+
[45,65)
+
0,61
+
0,34 – 1,10
+
0,10
+
+
+
[65,99]
+
0,46
+
0,23 – 0,92
+
0,027
+
+
Niveau d'étude
@@ -27077,52 +27867,22 @@
Exemple d’interaction
-
Femme * [16,25)
-
1,82
-
0,92 – 3,58
-
0,083
-
-
-
Homme * [16,25)
-
9,46
-
4,34 – 21,8
-
<0,001
-
-
-
Femme * [25,45)
-
1,86
-
1,17 – 3,00
-
0,009
-
-
Homme * [25,45)
-
3,02
-
1,87 – 4,96
-
<0,001
-
-
-
Femme * [45,65)
-
1,11
-
0,69 – 1,80
-
0,7
+
0,31
+
0,13 – 0,70
+
0,006
Homme * [45,65)
-
1,35
-
0,84 – 2,21
-
0,2
-
-
-
Femme * [65,99]
-
0,84
-
0,47 – 1,50
-
0,6
+
0,24
+
0,10 – 0,54
+
<0,001
Homme * [65,99]
-
-
-
+
0,23
+
0,09 – 0,59
+
0,003
@@ -27136,60 +27896,74 @@
Exemple d’interaction
-
Au sens strict, ce modèle explique tout autant le phénomène étudié que le modèle précédent. On peut le vérifier facilement avec anova.
-
anova(mod2, mod3, test ="Chisq")
+
Commençons par regarder les effets du modèle.
+
plot(allEffects(mod2))
+
+Représentation graphique des effets du modèle avec interaction entre le sexe et le groupe d’âge
+
Sur ce graphique, on voit que l’effet de l’âge sur la pratique d’un sport est surtout marqué chez les hommes. Chez les femmes, le même effet est observé, mais dans une moindre mesure et seulement à partir de 45 ans.
+
On peut tester si l’ajout de l’interaction améliore significativement le modèle avec anova.
+
anova(mod2, test ="Chisq")
-
De même, les effets modélisés sont les mêmes.
-
plot(allEffects(mod3))
-
-Représentation graphique des effets du modèle avec interaction simple entre le sexe et le groupe d’âge
-
Par contre, regardons d’un peu plus près les coefficients de ce nouveau modèle. Nous allons voir que leur interprétation est légèrement différente.
-
ggcoef_model(mod3, exponentiate =TRUE)
-
-Représentation graphique des coefficients du modèle avec interaction simple entre le sexe et le groupe d’âge
-
Cette fois-ci, il n’y a plus de coefficients globaux pour la variable sexe ni pour grpage mais des coefficients pour chaque combinaison de ces deux variables.
-
plot(
-broken(mod3, nouvelle_observation, predict.function = betas),
-trans = logit
-) +ylim(0, 1) +ylab("Probabilité de faire du sport")
+
Jetons maintenant un oeil aux coefficients du modèle. Pour rendre les choses plus visuelles, nous aurons recours à ggcoef_model de l’extension GGally.
+
+Représentation graphique des coefficients du modèle avec interaction entre le sexe et le groupe d’âge
+
Concernant l’âge et le sexe, nous avons trois séries de coefficients : trois coefficients (grpage[25,45), grpage[45,65) et grpage[65,99]) qui correspondent à l’effet global de la variable âge, un coefficient (sexeHomme)pour l’effet global du sexe et trois coefficients qui sont des moficateurs de l’effet d’âge pour les hommes (grpage[25,45), grpage[45,65) et grpage[65,99]).
+
Pour bien interpréter ces coefficients, il faut toujours avoir en tête les modalités choisies comme référence pour chaque variable. Supposons une femme de 60 ans, dont toutes lautres variables correspondent aux modalités de référence (c’est donc une pratiquante régulière, de niveau primaire, qui ne regarde pas la télévision). Regardons ce que prédit le modèle quant à sa probabilité de faire du sport au travers d’une représentation graphique
Scale for 'y' is already present. Adding another scale
for 'y', which will replace the existing scale.
-
-Représentation graphique de l’estimation de la probabilité de faire du sport pour un homme de 40 ans
-
Cette fois-ci, le coefficient d’interaction fourrnit l’effet global du sexe et de l’âge, et non plus la modification de cette combinaison par rapport aux coefficients globaux. Leur sens est donc différent et il faudra les interpréter en conséquence.
-
-
-
Un second exemple d’interaction
-
Intéressons-nous maintenant à l’interaction entre le sexe et le niveau d’étude. L’effet du niveau d’étude diffère-t-il selon le sexe ?
-
mod4 <-glm(sport ~ sexe * etud + grpage + heures.tv + relig, data = d, family =binomial())
-
Regardons d’abord les effets.
-
plot(allEffects(mod4))
-
-Représentation graphique des effets du modèle avec interaction entre le sexe et le niveau d’étude
-
À première vue, l’effet du niveau d’étude semble être le même chez les hommes et chez les femmes. Ceci dit, cela serait peut être plus lisible si l’on superposait les deux sexe sur un même graphique. Nous allons utiliser la fonction ggeffect de l’extension ggeffects qui permets de récupérer les effets calculés avec effect dans un format utilisable avec ggplot2.
-
library(ggeffects)
-
-Attaching package: 'ggeffects'
-
The following object is masked from 'package:cowplot':
-
- get_title
-
plot(ggeffect(mod4, c("etud", "sexe")))
-
-Effets du niveau d’étude selon le sexe
-
Cela confirme ce que l’on suppose. Regardons les coefficients du modèle.
-
ggcoef_model(mod4, exponentiate =TRUE)
-
-Représentation graphique des coefficients du modèle avec interaction simple entre le sexe et le niveau d’étude
-
tbl_regression(mod4, exponentiate =TRUE)
+
+Représentation graphique de l’estimation de la probabilité de faire du sport pour une femme de 60 ans
+
En premier lieu, l’intercept s’applique et permet de déterminer la probabilité de base de faire du sport (si toutes les variables sont à leur valeur de référence). Femme étant la modalité de référence pour la variable sexe, cela ne modifie pas le calcul de la probabilité de faire du sport. Par contre, il y a une modification induite par la modalité 45-65 de la variable grpage.
+
Regardons maintenant la situation d’un homme de 20 ans.
+
nouvelle_observation$sexe[1] <-"Homme"
+nouvelle_observation$grpage[1] <-"[16,25)"
+plot(
+broken(mod2, nouvelle_observation, predict.function = betas),
+trans = logit
+) +ylim(0, 1) +ylab("Probabilité de faire du sport")
+
Scale for 'y' is already present. Adding another scale
+for 'y', which will replace the existing scale.
+
+Représentation graphique de l’estimation de la probabilité de faire du sport pour un homme de 20 ans
+
Nous sommes à la modalité de référence pour l’âge par contre il y a un effet important du sexe. Le coefficient associé globalement à la variable sexe correspond donc à l’effet du sexe à la modalité de référence du groupe d’âges.
+
La situation est différente pour un homme de 60 ans.
+
nouvelle_observation$grpage[1] <-"[45,65)"
+plot(
+broken(mod2, nouvelle_observation, predict.function = betas),
+trans = logit
+) +ylim(0, 1) +ylab("Probabilité de faire du sport")
+
Scale for 'y' is already present. Adding another scale
+for 'y', which will replace the existing scale.
+
+Représentation graphique de l’estimation de la probabilité de faire du sport pour un homme de 60 ans
+
Cette fois-ci, il y a plusieurs modifications d’effet. On applique en effet à la fois le coefficient sexe = Homme (effet du sexe pour les 15-24 ans), le coefficient grpage = [45-65) qui est l’effet de l’âge pour les femmes de 45-64 ans et le coefficient sexe:grpage = Homme:[45-65) qui indique l’effet spécifique qui s’applique aux hommes de 45-64, d’une part par rapport aux femmes du même et d’autre part par rapport aux hommes de 16-24 ans. L’effet des coefficients d’interaction doivent donc être interprétés par rapport aux autres coefficients du modèle qui s’appliquent, en tenant compte des modalités de référence.
+
Il est cependant possible d’écrire le même modèle différemment. En effet, sexe * grpage dans la formule du modèle est équivalent à l’écriture sexe + grpage + sexe:grpage, c’est-à-dire à modéliser un coefficient global pour chaque variable plus un des coefficients d’interaction. On aurait pu demander juste des coefficients d’interaction, en ne mettant que sexe:grpage.
+
mod3 <-glm(sport ~ sexe:grpage + etud + heures.tv + relig, data = d, family =binomial())
+tbl_regression(mod3, exponentiate =TRUE)
-
+
Caractéristique
@@ -27534,24 +28308,6 @@
Un second exemple d’interaction
-
Sexe
-
-
-
-
-
-
Femme
-
—
-
—
-
-
-
-
Homme
-
1,42
-
0,77 – 2,58
-
0,3
-
-
Niveau d'étude
@@ -27565,56 +28321,26 @@
Un second exemple d’interaction
Secondaire
-
2,25
-
1,35 – 3,82
-
0,002
+
2,60
+
1,77 – 3,85
+
<0,001
Technique/Professionnel
-
3,09
-
1,90 – 5,17
+
2,88
+
1,99 – 4,21
<0,001
Supérieur
-
6,59
-
4,02 – 11,1
+
6,77
+
4,64 – 10,0
<0,001
Manquant
-
5,21
-
2,42 – 11,4
-
<0,001
-
-
-
Groupe d'âges
-
-
-
-
-
-
[16,25)
-
—
-
—
-
-
-
-
[25,45)
-
0,66
-
0,42 – 1,04
-
0,071
-
-
-
[45,65)
-
0,34
-
0,21 – 0,55
-
<0,001
-
-
-
[65,99]
-
0,25
-
0,15 – 0,43
+
9,02
+
4,68 – 17,8
<0,001
@@ -27637,63 +28363,87 @@
Un second exemple d’interaction
Pratiquant occasionnel
-
0,98
-
0,67 – 1,42
-
0,9
+
0,99
+
0,68 – 1,43
+
>0,9
Appartenance sans pratique
-
1,01
-
0,72 – 1,43
-
>0,9
+
1,02
+
0,73 – 1,44
+
0,9
Ni croyance ni appartenance
0,83
-
0,56 – 1,21
+
0,57 – 1,22
0,3
Rejet
-
0,70
-
0,40 – 1,23
+
0,68
+
0,38 – 1,19
0,2
NSP ou NVPR
-
0,90
-
0,39 – 1,98
-
0,8
+
0,93
+
0,40 – 2,04
+
0,9
-
Sexe * Niveau d'étude
+
Sexe * Groupe d'âges
-
Homme * Secondaire
-
1,38
-
0,65 – 2,93
-
0,4
+
Femme * [16,25)
+
1,82
+
0,92 – 3,58
+
0,083
-
Homme * Technique/Professionnel
-
0,87
-
0,44 – 1,75
+
Homme * [16,25)
+
9,46
+
4,34 – 21,8
+
<0,001
+
+
+
Femme * [25,45)
+
1,86
+
1,17 – 3,00
+
0,009
+
+
+
Homme * [25,45)
+
3,02
+
1,87 – 4,96
+
<0,001
+
+
+
Femme * [45,65)
+
1,11
+
0,69 – 1,80
0,7
-
Homme * Supérieur
-
1,01
-
0,49 – 2,07
-
>0,9
+
Homme * [45,65)
+
1,35
+
0,84 – 2,21
+
0,2
-
Homme * Manquant
-
4,28
-
1,32 – 15,8
-
0,020
+
Femme * [65,99]
+
0,84
+
0,47 – 1,50
+
0,6
+
+
+
Homme * [65,99]
+
+
+
@@ -27707,411 +28457,832 @@
Un second exemple d’interaction
-
Si les coefficients associés au niveau d’étude sont significatifs, ceux de l’interaction ne le sont pas (sauf sexeHomme:etudManquant) et celui associé au sexe, précédemment significatif ne l’est plus. Testons avec anova si l’interaction est belle et bien significative.
-
anova(mod4, test ="Chisq")
+
Au sens strict, ce modèle explique tout autant le phénomène étudié que le modèle précédent. On peut le vérifier facilement avec anova.
+
anova(mod2, mod3, test ="Chisq")
-
L’interaction est bien significative mais faiblement. Vu que l’effet du niveau d’étude reste nénamoins très similaire selon le sexe, on peut se demander s’il est pertinent de la conserver.
+
De même, les effets modélisés sont les mêmes.
+
plot(allEffects(mod3))
+
+Représentation graphique des effets du modèle avec interaction simple entre le sexe et le groupe d’âge
+
Par contre, regardons d’un peu plus près les coefficients de ce nouveau modèle. Nous allons voir que leur interprétation est légèrement différente.
+
ggcoef_model(mod3, exponentiate =TRUE)
+
+Représentation graphique des coefficients du modèle avec interaction simple entre le sexe et le groupe d’âge
+
Cette fois-ci, il n’y a plus de coefficients globaux pour la variable sexe ni pour grpage mais des coefficients pour chaque combinaison de ces deux variables.
+
plot(
+broken(mod3, nouvelle_observation, predict.function = betas),
+trans = logit
+) +ylim(0, 1) +ylab("Probabilité de faire du sport")
+
Scale for 'y' is already present. Adding another scale
+for 'y', which will replace the existing scale.
+
+Représentation graphique de l’estimation de la probabilité de faire du sport pour un homme de 40 ans
+
Cette fois-ci, le coefficient d’interaction fourrnit l’effet global du sexe et de l’âge, et non plus la modification de cette combinaison par rapport aux coefficients globaux. Leur sens est donc différent et il faudra les interpréter en conséquence.
-
-
Explorer les différentes interactions possibles
-
Il peut y avoir de multiples interactions dans un modèle, d’ordre 2 (entre deux variables) ou plus (entre trois variables ou plus). Il est dès lors tentant de tester les multiples interactions possibles de manière itératives afin d’identifier celles à retenir. C’est justement le but de la fonction glmulti de l’extension du même nom. glmulti permets de tester toutes les combinaisons d’interactions d’ordre 2 dans un modèle, en retenant le meilleur modèle à partir d’un critère spécifié (par défaut l’AIC). ATTENTION : le temps de calcul de glmulti peut-être long.
-
library(glmulti)
-glmulti(sport ~ sexe + grpage + etud + heures.tv + relig, data = d, family =binomial())
-
Initialization...
-TASK: Exhaustive screening of candidate set.
-Fitting...
+
+
Un second exemple d’interaction
+
Intéressons-nous maintenant à l’interaction entre le sexe et le niveau d’étude. L’effet du niveau d’étude diffère-t-il selon le sexe ?
+
mod4 <-glm(sport ~ sexe * etud + grpage + heures.tv + relig, data = d, family =binomial())
+
Regardons d’abord les effets.
+
plot(allEffects(mod4))
+
+Représentation graphique des effets du modèle avec interaction entre le sexe et le niveau d’étude
+
À première vue, l’effet du niveau d’étude semble être le même chez les hommes et chez les femmes. Ceci dit, cela serait peut être plus lisible si l’on superposait les deux sexe sur un même graphique. Nous allons utiliser la fonction ggeffect de l’extension ggeffects qui permets de récupérer les effets calculés avec effect dans un format utilisable avec ggplot2.
+
library(ggeffects)
+
+Attaching package: 'ggeffects'
+
The following object is masked from 'package:cowplot':
-After 50 models:
-Best model: sport~1+grpage+heures.tv+sexe:heures.tv+grpage:heures.tv+etud:heures.tv
-Crit= 2284.87861987263
-Mean crit= 2406.80086471225
+ get_title
+
plot(ggeffect(mod4, c("etud", "sexe")))
+
+Effets du niveau d’étude selon le sexe
+
Cela confirme ce que l’on suppose. Regardons les coefficients du modèle.
+
ggcoef_model(mod4, exponentiate =TRUE)
+
+Représentation graphique des coefficients du modèle avec interaction simple entre le sexe et le niveau d’étude
+
tbl_regression(mod4, exponentiate =TRUE)
+
+
+
+
Caractéristique
+
+OR1
+
+
+95% CI1
+
+
p-value
+
+
+
+
Sexe
+
+
+
+
+
+
Femme
+
—
+
—
+
+
+
+
Homme
+
1,42
+
0,77 – 2,58
+
0,3
+
+
+
Niveau d'étude
+
+
+
+
+
+
Primaire
+
—
+
—
+
+
+
+
Secondaire
+
2,25
+
1,35 – 3,82
+
0,002
+
+
+
Technique/Professionnel
+
3,09
+
1,90 – 5,17
+
<0,001
+
+
+
Supérieur
+
6,59
+
4,02 – 11,1
+
<0,001
+
+
+
Manquant
+
5,21
+
2,42 – 11,4
+
<0,001
+
+
+
Groupe d'âges
+
+
+
+
+
+
[16,25)
+
—
+
—
+
+
+
+
[25,45)
+
0,66
+
0,42 – 1,04
+
0,071
+
+
+
[45,65)
+
0,34
+
0,21 – 0,55
+
<0,001
+
+
+
[65,99]
+
0,25
+
0,15 – 0,43
+
<0,001
+
+
+
Nombre d'heures passées devant la télévision par jour
+
0,89
+
0,83 – 0,95
+
<0,001
+
+
+
Pratique religieuse
+
+
+
+
+
+
Pratiquant regulier
+
—
+
—
+
+
+
+
Pratiquant occasionnel
+
0,98
+
0,67 – 1,42
+
0,9
+
+
+
Appartenance sans pratique
+
1,01
+
0,72 – 1,43
+
>0,9
+
+
+
Ni croyance ni appartenance
+
0,83
+
0,56 – 1,21
+
0,3
+
+
+
Rejet
+
0,70
+
0,40 – 1,23
+
0,2
+
+
+
NSP ou NVPR
+
0,90
+
0,39 – 1,98
+
0,8
+
+
+
Sexe * Niveau d'étude
+
+
+
+
+
+
Homme * Secondaire
+
1,38
+
0,65 – 2,93
+
0,4
+
+
+
Homme * Technique/Professionnel
+
0,87
+
0,44 – 1,75
+
0,7
+
+
+
Homme * Supérieur
+
1,01
+
0,49 – 2,07
+
>0,9
+
+
+
Homme * Manquant
+
4,28
+
1,32 – 15,8
+
0,020
+
+
+
+
+
+ 1
+
+
+ OR = rapport de cotes, CI = intervalle de confiance
+
+
+
+
+
Si les coefficients associés au niveau d’étude sont significatifs, ceux de l’interaction ne le sont pas (sauf sexeHomme:etudManquant) et celui associé au sexe, précédemment significatif ne l’est plus. Testons avec anova si l’interaction est belle et bien significative.
+
anova(mod4, test ="Chisq")
-
Reste à calculer si le rang de naissance de l’enfant est supérieur au nombre idéal d’enfants tel que défini par la mère. On aura recours à la fonction rank appliquée par groupe (ici calculé séparément pour chaque mère). L’argument ties.method permet d’indiquer comment gérer les égalités (ici les naissances multiples, e.g. les jumeaux). Comme nous voulons comparer le rang de l’enfant au nombre idéal d’enfants, nous allons retenir la méthode "max" pour obtenir, dans le cas présent, le nombre total d’enfants déjà nés1. Avant de calculer un rang, il est impératif de trier préalablement le tableau (voir le chapitre Tris).
-
setorder(enfants, id_femme, date_naissance)
-enfants[, rang :=rank(date_naissance, ties.method ="max"), by = id_femme]
-enfants[, rang_apres_ideal :="non"]
-# note: unclass() requis en raison d'un bug non corrigé dans haven empéchant de comparer haven_labelled_spss et integer
-enfants[rang >unclass(nb_enf_ideal), rang_apres_ideal :="oui"]
-enfants[, rang_apres_ideal :=factor(rang_apres_ideal)]
-enfants[, rang_apres_ideal :=relevel(rang_apres_ideal, "non")]
+
L’interaction est bien significative mais faiblement. Vu que l’effet du niveau d’étude reste nénamoins très similaire selon le sexe, on peut se demander s’il est pertinent de la conserver.
-
-
Préparation des données avec dplyr
+
+
Explorer les différentes interactions possibles
+
Il peut y avoir de multiples interactions dans un modèle, d’ordre 2 (entre deux variables) ou plus (entre trois variables ou plus). Il est dès lors tentant de tester les multiples interactions possibles de manière itératives afin d’identifier celles à retenir. C’est justement le but de la fonction glmulti de l’extension du même nom. glmulti permets de tester toutes les combinaisons d’interactions d’ordre 2 dans un modèle, en retenant le meilleur modèle à partir d’un critère spécifié (par défaut l’AIC). ATTENTION : le temps de calcul de glmulti peut-être long.
+
library(glmulti)
+glmulti(sport ~ sexe + grpage + etud + heures.tv + relig, data = d, family =binomial())
On voit qu’au bout d’un moment, l’algorithme se statibilise autour d’un modèle comportant une interaction entre le sexe et l’âge d’une part et entre le sexe et le nombre d’heures passées quotidiennement devant la télé. On voit également que la variable religion a été retirée du modèle final.
+
best <-glm(sport ~1+ sexe + grpage + etud + heures.tv + grpage:sexe + sexe:heures.tv, data = d, family =binomial())
+odds.ratio(best)
+
Waiting for profiling to be done...
+
+
+
+
ggcoef_model(best, exponentiate =TRUE)
+
+Représentation graphique des coefficients du modèle avec interaction entre le sexe, le niveau d’étude et le nombre d’heures passées devant la télévision
+
plot(allEffects(best))
+
+Représentation graphique des effets du modèle avec interaction entre le sexe, le niveau d’étude et le nombre d’heures passées devant la télévision
+
+
+
+
Pour aller plus loin
+
Il y a d’autres extensions dédiées à l’analyse des interactions d’un modèle, de même que de nombreux supports de cours en ligne dédiés à cette question.
Ce chapitre est évoqué dans le webin-R #07 (régression logistique partie 2) sur YouTube.
+
+
Dans une régression, la multicolinéarité est un problème qui survient lorsque certaines variables de prévision du modèle mesurent le même phénomène. Une multicolinéarité prononcée s’avère problématique, car elle peut augmenter la variance des coefficients de régression et les rendre instables et difficiles à interpréter. Les conséquences de coefficients instables peuvent être les suivantes :
+
+
les coefficients peuvent sembler non significatifs, même lorsqu’une relation significative existe entre le prédicteur et la réponse ;
+
les coefficients de prédicteurs fortement corrélés varieront considérablement d’un échantillon à un autre ;
+
lorsque des termes d’un modèle sont fortement corrélés, la suppression de l’un de ces termes aura une incidence considérable sur les coefficients estimés des autres. Les coefficients des termes fortement corrélés peuvent même présenter le mauvais signe.
+
+
La multicolinéarité n’a aucune incidence sur l’adéquation de l’ajustement, ni sur la qualité de la prévision. Cependant, les coefficients individuels associés à chaque variable explicative ne peuvent pas être interprétés de façon fiable.
+
+
Définition
+
Au sens strict, on parle de multicolinéarité parfaite lorsqu’une des variables explicatives d’un modèle est une combinaison linéraire d’une ou plusieurs autres variables explicatives introduites dans le même modèle. L’absence de multicolinéarité parfaite est une des conditions requises pour pouvoir estimer un modèle linéaire et, par extension, un modèle linéaire généralisé (dont les modèles de régression logistique).
+
Dans les faits, une multicolinéarité parfaite n’est quasiment jamais observée. Mais une forte multicolinéarité entre plusieurs variables peut poser problème dans l’estimation et l’interprétation d’un modèle.
+
Une erreur fréquente est de confondre multicolinéarité et corrélation. Si des variables colinéaires sont de facto fortement corrélées entre elles, deux variables corrélées ne sont pas forcément colinéaires. En termes non statistiques, il y a colinéarité lorsque deux ou plusieurs variables mesurent la même chose.
+
Prenons un exemple. Nous étudions les complications après l’accouchement dans différentes maternités d’un pays en développement. On souhaite mettre dans le modèle, à la fois le milieu de résidence (urbain ou rural) et le fait qu’il y ait ou non un médecin dans la clinique. Or, dans la zone d’enquête, les maternités rurales sont dirigées seulement par des sage-femmes tandis que l’on trouve un médecin dans toutes les maternités urbaines sauf une. Dès lors, dans ce contexte précis, le milieu de résidence prédit presque totalement la présence d’un médecin et on se retrouve face à une multicolinéarité (qui serait même parfaite s’il n’y avait pas une clinique urbaine sans médecin). On ne peut donc distinguer l’effet de la présence d’un médecin de celui du milieu de résidence et il ne faut mettre qu’une seule de ces deux variables dans le modèle, sachant que du point de vue de l’interprétation elle capturera à la fois l’effet de la présence d’un médecin et celui du milieu de résidence.
+
Par contre, si dans notre région d’étude, seule la moitié des maternités urbaines disposait d’un médecin, alors le milieu de résidence n’aurait pas été suffisant pour prédire la présence d’un médecin. Certes, les deux variables seraient corrélées mais pas colinéaires. Un autre exemple de corrélation sans colinéarité, c’est la relation entre milieu de résidence et niveau d’instruction. Il y a une corrélation entre ces deux variables, les personnes résidant en ville étant généralement plus instruites. Cependant, il existe également des personnes non instruites en ville et des personnes instruites en milieu rural. Le milieu de résidence n’est donc pas suffisant pour prédire le niveau d’instruction.
+
+
+
Mesure de la colinéarité
+
Il existe différentes mesures de la multicolinéarité. L’extension mctest en fournie plusieurs, mais elle n’est utilisable que si l’ensemble des variables explicatives sont de type numérique.
+
L’approche la plus classique consiste à examiner les facteurs d’inflation de la variance (FIV) ou variance inflation factor (VIF) en anglais. Les FIV estimenent de combien la variance d’un coefficient est augmentée en raison d’une relation linéaire avec d’autres prédicteurs. Ainsi, un FIV de 1,8 nous dit que la variance de ce coefficient particulier est supérieure de 80 % à la variance que l’on aurait dû observer si ce facteur n’est absolument pas corrélé aux autres prédicteurs.
+
Si tous les FIV sont égaux à 1, il n’existe pas de multicolinéarité, mais si certains FIV sont supérieurs à 1, les prédicteurs sont corrélés. Il n’y a pas de consensus sur la valeur au-delà de laquelle on doit considérer qu’il y a multicolinéarité. Certains auteurs, comme Paul Allison1, disent regarder plus en détail les variables avec un FIV supérieur à 2,5. D’autres ne s’inquiètent qu’à partir de 5. Il n’existe pas de test statistique qui permettrait de dire s’il y a colinéarité ou non2.
+
L’extension car fournit une fonction vif permettant de calculer les FIV à partir d’un modèle. Elle implémente même une version généralisée permettant de considérer des facteurs catégoriels et des modèles linéaires généralisés comme la régression logistique.
+
Reprenons, pour exemple, un modèle logistique que nous avons déjà abordé dans d’autres chapitres.
Dans notre exemple, tous les FIV sont proches de 1. Il n’y a donc pas de problème potentiel de colinéarité à explorer.
+
+
+
La multicolinéarité est-elle toujours un problème ?
+
Là encore, il n’y a pas de consensus sur cette question. Certains analystes considèrent que tout modèle où certains prédicteurs seraient colinéaires n’est pas valable. Dans un billet sur internet, Paul Allison évoque quant à lui des situations où la multicolinéarité peut être ignorée en toute sécurité. Le texte ci-dessous est une traduction du billet de Paul Allison.
+
1. Les variables avec des FIV élevés sont des variables de contrôle, et les variables d’intérêt n’ont pas de FIV élevés.
+
Voici le problème de la multicollinéarité : ce n’est un problème que pour les variables qui sont colinéaires. Il augmente les erreurs-types de leurs coefficients et peut rendre ces coefficients instables de plusieurs façons. Mais tant que les variables colinéaires ne sont utilisées que comme variables de contrôle, et qu’elles ne sont pas colinéaires avec vos variables d’intérêt, il n’y a pas de problème. Les coefficients des variables d’intérêt ne sont pas affectés et la performance des variables de contrôle n’est pas altérée.
+
Voici un exemple tiré de ces propres travaux : l’échantillon est constitué de collèges américains, la variable dépendante est le taux d’obtention de diplôme et la variable d’intérêt est un indicateur (factice) pour les secteurs public et privé. Deux variables de contrôle sont les scores moyens au SAT et les scores moyens à l’ACT pour l’entrée en première année. Ces deux variables ont une corrélation supérieure à ,9, ce qui correspond à des FIV d’au moins 5,26 pour chacune d’entre elles. Mais le FIV pour l’indicateur public/privé n’est que de 1,04. Il n’y a donc pas de problème à se préoccuper et il n’est pas nécessaire de supprimer l’un ou l’autre des deux contrôles, à condition que l’on ne cherche pas à interpréter ou comparer l’un par rapport à l’autre les coefficients de ces deux variables de contrôle.
+
2. Les FIV élevés sont causés par l’inclusion de puissances ou de produits d’autres variables.
+
Si vous spécifiez un modèle de régression avec x et x2, il y a de bonnes chances que ces deux variables soient fortement corrélées. De même, si votre modèle a x, z et xz, x et z sont susceptibles d’être fortement corrélés avec leur produit. Il n’y a pas de quoi s’inquiéter, car la valeur p de xz n’est pas affectée par la multicollinéarité. Ceci est facile à démontrer : vous pouvez réduire considérablement les corrélations en centrant les variables (c’est-à-dire en soustrayant leurs moyennes) avant de créer les puissances ou les produits. Mais la valeur p pour x2 ou pour xz sera exactement la même, que l’on centre ou non. Et tous les résultats pour les autres variables (y compris le R2 mais sans les termes d’ordre inférieur) seront les mêmes dans les deux cas. La multicollinéarité n’a donc pas de conséquences négatives.
+
3. Les variables avec des FIV élevés sont des variables indicatrices (factices) qui représentent une variable catégorielle avec trois catégories ou plus.
+
Si la proportion de cas dans la catégorie de référence est faible, les variables indicatrices auront nécessairement des FIV élevés, même si la variable catégorielle n’est pas associée à d’autres variables dans le modèle de régression.
+
Supposons, par exemple, qu’une variable de l’état matrimonial comporte trois catégories : actuellement marié, jamais marié et anciennement marié. Vous choisissez anciennement marié comme catégorie de référence, avec des variables d’indicateur pour les deux autres. Ce qui se passe, c’est que la corrélation entre ces deux indicateurs devient plus négative à mesure que la fraction de personnes dans la catégorie de référence diminue. Par exemple, si 45 % des personnes ne sont jamais mariées, 45 % sont mariées et 10 % sont anciennement mariées, les valeurs du FIV pour les personnes mariées et les personnes jamais mariées seront d’au moins 3,0.
+
Est-ce un problème ? Eh bien, cela signifie que les valeurs p des variables indicatrices peuvent être élevées. Mais le test global selon lequel tous les indicateurs ont des coefficients de zéro n’est pas affecté par des FIV élevés. Et rien d’autre dans la régression n’est affecté. Si vous voulez vraiment éviter des FIV élevés, il suffit de choisir une catégorie de référence avec une plus grande fraction des cas. Cela peut être souhaitable pour éviter les situations où aucun des indicateurs individuels n’est statistiquement significatif, même si l’ensemble des indicateurs est significatif.
Ce chapitre est évoqué dans le webin-R #15 (analyse de survie) sur YouTube.
+
+
+
Ressources en ligne
+
L’extension centrale pour l’analyse de survie est survival.
+
Un très bon tutoriel (en anglais et en 3 étapes), introduisant les concepts de l’analyse de survie, des courbes de Kaplan-Meier et des modèles de Cox et leur mise en oeuvre pratique sous R est disponible en ligne :
A noter, il est possible d’utiliser la fonction step sur un modèle de Cox, pour une sélection pas à pas d’un meilleur modèle basé sur une minimisation de l’AIC (voir le chapitre sur la régression logistique).
+
L’excellente extension broom peut également être utilisée sur des modèles de survie (Kaplan-Meier ou Cox) pour en convertir les résultats sous la forme d’un tableau de données.
+
Pour approfondir les possibilités offertes par l’extension survival, on pourra également consulter les différentes vignettes fournies avec l’extension (voir https://cran.r-project.org/package=survival).
+
+
+
Un exemple concret : mortalité infanto-juvénile
+
Dans cet exemple, nous allons utiliser le jeu de données fecondite fourni par l’extension questionr. Ce jeu de données comporte trois tableaux de données : menages, femmes et enfants.
+
Nous souhaitons étudier ici la survie des enfants entre la naissance et l’âge de 5 ans. Dans un premier temps, nous comparerons la survie des jeunes filles et des jeunes garçons. Dans un second temps, nous procéderons à une analyse multivariée en prenant en compte les variables suivantes :
+
+
sexe de l’enfant
+
milieu de résidence
+
niveau de vie du ménage
+
structure du ménage
+
niveau d’éducation de la mère
+
âge de la mère à la naissance de l’enfant
+
enfin, une variable un peu plus compliquée, à savoir si le rang de naissance de l’enfant (second, troisième, quatrième, etc.) est supérieur au nombre idéal d’enfants selon la mère.
+
+
Nous allons préparer les données selon deux approches : soit en utilisant l’extension data.table (voir le chapitre dédié à data.table), soit en utilisant l’extension dplyr (voir le chapitre sur dplyr).
+
Chargeons les données en mémoire et listons les variables disponibles.
Les tableaux de données sont déjà au format tibble (c’est-à-dire sont de la classe tbl_df)2 et les variables catégorielles sont du type labelled (voir le chapitre sur les vecteurs labellisés). Ce format correspond au format de données si on les avait importées depuis SPSS avec l’extension haven (voir le chapitre sur l’import de données).
-
Nous allons charger en mémoire l’extension labelled pour la gestion des vecteurs labellisés en plus de dplyr.
-
library(dplyr)
-library(labelled)
+
Les tableaux de données sont au format tibble (c’est-à-dire sont de la classe tbl_df) et les variables catégorielles sont du type haven_labelled (voir le chapitre sur les vecteurs labellisés). Ce format correspond au format de données si on les avait importées depuis SPSS avec l’extension haven (voir le chapitre sur l’import de données).
+
En premier lieu, il nous faut convertir les tableaux de données au format data.table, ce qui peut se faire avec la fonction setDT. Par ailleurs, nous allons également charger en mémoire l’extension labelled pour la gestion des vecteurs labellisés.
En premier lieu, il nous faut calculer la durée d’observation des enfants, à savoir le temps passé entre la date de naissance (variable du fichier enfants) et la date de passation de l’entretien (fournie par le tableau de données femmes). Pour récupérer des variables du fichier femmes dans le fichier enfants, nous allons procéder à une fusion de table (voir le chapitre dédié). Pour le calcul de la durée d’observation, nous allons utiliser le package lubridate (voir le chapitre calculer un âge et celui sur la gestion des dates). Nous effectuerons l’analyse en mois (puisque l’âge au décès est connu en mois). Dès lors, la durée d’observation sera calculée en mois.
-
ATTENTION : les étiquettes de valeurs sont le plus souvent perdus lors des fusions de tables avec dplyr. Pour les récupérer après fusion, nous allons conserver une version originale du tableau de données enfants et utiliser la fonction copy_labels_from de l’extension labelled.
enfants <-merge(
+ enfants,
+ femmes[, .(id_femme, date_entretien)],
+by ="id_femme",
+all.x =TRUE
+)
+
+# duree observation en mois
+library(lubridate, quietly =TRUE)
+enfants[, duree_observation :=time_length(interval(date_naissance, date_entretien), unit ="months")]
+
ATTENTION : il y 11 enfants soi-disant nés après la date d’enquête ! Quelle que soit l’enquête, il est rare de ne pas observer d’incohérences. Dans le cas présent, il est fort possible que la date d’entretien puisse parfois être erronnée (par exemple si l’enquêteur a inscrit une date sur le questionnaire papier le jour du recensement du ménage mais n’ai pu effectué le questionnaire individuel que plus tard). Nous décidons ici de procéder à une correction en ajoutant un mois aux dates d’entretien problématiques. D’autres approches auraient pu être envisagées, comme par exemple exclure ces observations problématiques. Cependant, cela aurait impacté le calcul du range de naissance pour les autres enfants issus de la même mère. Quoiqu’il en soit, il n’y a pas de réponse unique. À vous de vous adapter au contexte particulier de votre analyse.
Regardons maintenant comment les âges au décès ont été collectés.
+
freq(enfants$age_deces)
+
+
+
+
Les âges au décès sont ici exprimés en mois révolus. Les décès à un mois révolu correspondent à des décès entre 1 et 2 mois exacts. Par ailleurs, les durées d’observation que nous avons calculées avec time_length sont des durées exactes, c’est-à-dire avec la partie décimale. Pour une analyse de survie, on ne peut mélanger des durées exactes et des durées révolues. Trois approches peuvent être envisagées :
+
+
faire l’analyse en mois révolus, auquel cas on ne gardera que la partie entière des durées d’observations avec la fonction trunc ;
+
considérer qu’un âge au décès de 3 mois révolus correspond en moyenne à 3,5 mois exacts et donc ajouter 0,5 à tous les âges révolus ;
+
imputer un âge au décès exact en distribuant aléatoirement les décès à 3 mois révolus entre 3 et 4 mois exacts, autrement dit en ajoutant aléatoirement une partie décimale aux âges révolus.
+
+
Nous allons ici adopter la troisième approche en considérant que les décès se répartissent de manière uniforme au sein d’un même mois. Nous aurons donc recours à la fonction runif qui permets de générer des valeurs aléatoires entre 0 et 1 selon une distribustion uniforme.
Pour définir notre objet de survie, il nous faudra deux variables. Une première, temporelle, indiquant la durée à laquelle survient l’évènement étudié (ici le décès) pour ceux ayant vécu l’évènement et la durée d’observation pour ceux n’ayant pas vécu l’évènement (censure à droite). Par ailleurs, une seconde variable indiquant si les individus ont vécu l’évènement (0 pour non, 1 pour oui). Or, ici, la variable survie est codée 0 pour les décès et 1 pour ceux ayant survécu. Pour plus de détails, voir l’aide de la fonction Surv.
Occupons-nous maintenant des variables explicatives que nous allons inclure dans l’analyse. Tout d’abord, ajoutons à la table enfants les variables nécessaires des tables femmes et menages. Notons qu’il nous faudra importer id_menage de la table femmes pour pouvoir fusionner ensuite la table enfants avec la table menages. Par ailleurs, pour éviter une confusion sur la variable date_naissance, nous renommons à la volée cette variable de la table femmes en date_naissance_mere.
Les variables catégorielles sont pour l’heure sous formes de vecteurs labellisés. Or, dans un modèle, il est impératif de les convertir en facteurs pour qu’elles soient bien traitées comme des variables catégorielles (autrement elles seraient traitées comme des variables continues). On aura donc recours à la fonction to_factor de l’extension labelled.
Regardons plus attentivement, la variable structure.
+
freq(enfants$structure)
+
+
+
+
Tout d’abord, la modalité pas d’adulte n’est pas représentée dans l’échantillon. On aura donc recours à l’argument drop_unused_labels pour ne pas conserver cette modalité. Par ailleurs, nous considérons que la situation familiale à partir de laquelle nous voudrons comparer les autres dans notre modèle, donc celle qui doit être considérée comme la modalité de référence, est celle du ménage nucléaire. Cette modalité (deux adultes de sexe opposé) n’étant pas la première, nous aurons recours à la fonction relevel.
La modalité supérieur est peu représentée dans notre échantillon. Nous allons la fusionner avec la modalité secondaire (voir la section Regrouper les modalités d’une variable du chapitre Recodage).
Reste à calculer si le rang de naissance de l’enfant est supérieur au nombre idéal d’enfants tel que défini par la mère. On aura recours à la fonction rank appliquée par groupe (ici calculé séparément pour chaque mère). L’argument ties.method permet d’indiquer comment gérer les égalités (ici les naissances multiples, e.g. les jumeaux). Comme nous voulons comparer le rang de l’enfant au nombre idéal d’enfants, nous allons retenir la méthode "max" pour obtenir, dans le cas présent, le nombre total d’enfants déjà nés1. Avant de calculer un rang, il est impératif de trier préalablement le tableau (voir le chapitre Tris).
+
setorder(enfants, id_femme, date_naissance)
+enfants[, rang :=rank(date_naissance, ties.method ="max"), by = id_femme]
+enfants[, rang_apres_ideal :="non"]
+# note: unclass() requis en raison d'un bug non corrigé dans haven empéchant de comparer haven_labelled_spss et integer
+enfants[rang >unclass(nb_enf_ideal), rang_apres_ideal :="oui"]
+enfants[, rang_apres_ideal :=factor(rang_apres_ideal)]
+enfants[, rang_apres_ideal :=relevel(rang_apres_ideal, "non")]
+
+
+
Préparation des données avec dplyr
+
Tout d’abord, regardons sous quel format elles sont stockées.
+
data(fecondite)
+class(menages)
+
[1] "tbl_df" "tbl" "data.frame"
+
describe(menages)
+
[1814 obs. x 5 variables] tbl_df tbl data.frame
+
+$id_menage: Identifiant du ménage
+numeric: 1 2 3 4 5 6 7 8 9 10 ...
+min: 1 - max: 1814 - NAs: 0 (0%) - 1814 unique values
+
+$taille: Taille du ménage (nombre de membres)
+numeric: 7 3 6 5 7 6 15 6 5 19 ...
+min: 1 - max: 31 - NAs: 0 (0%) - 30 unique values
+
+$sexe_chef: Sexe du chef de ménage
+labelled double: 2 1 1 1 1 2 2 2 1 1 ...
+min: 1 - max: 2 - NAs: 0 (0%) - 2 unique values
+2 value labels: [1] homme [2] femme
+
+$structure: Structure démographique du ménage
+labelled double: 4 2 5 4 4 4 5 2 5 5 ...
+min: 1 - max: 5 - NAs: 0 (0%) - 5 unique values
+6 value labels: [0] pas d'adulte [1] un adulte [2] deux adultes de sexe opposé [3] deux adultes de même sexe [4] trois adultes ou plus avec lien de parenté [5] adultes sans lien de parenté
+
+$richesse: Niveau de vie (quintiles)
+labelled double: 1 2 2 1 1 3 2 5 4 3 ...
+min: 1 - max: 5 - NAs: 0 (0%) - 5 unique values
+5 value labels: [1] très pauvre [2] pauvre [3] moyen [4] riche [5] très riche
+
Les tableaux de données sont déjà au format tibble (c’est-à-dire sont de la classe tbl_df)2 et les variables catégorielles sont du type labelled (voir le chapitre sur les vecteurs labellisés). Ce format correspond au format de données si on les avait importées depuis SPSS avec l’extension haven (voir le chapitre sur l’import de données).
+
Nous allons charger en mémoire l’extension labelled pour la gestion des vecteurs labellisés en plus de dplyr.
+
library(dplyr)
+library(labelled)
+
En premier lieu, il nous faut calculer la durée d’observation des enfants, à savoir le temps passé entre la date de naissance (variable du fichier enfants) et la date de passation de l’entretien (fournie par le tableau de données femmes). Pour récupérer des variables du fichier femmes dans le fichier enfants, nous allons procéder à une fusion de table (voir le chapitre dédié). Pour le calcul de la durée d’observation, nous allons utiliser le package lubridate (voir le chapitre calculer un âge et celui sur la gestion des dates). Nous effectuerons l’analyse en mois (puisque l’âge au décès est connu en mois). Dès lors, la durée d’observation sera calculée en mois.
+
ATTENTION : les étiquettes de valeurs sont le plus souvent perdus lors des fusions de tables avec dplyr. Pour les récupérer après fusion, nous allons conserver une version originale du tableau de données enfants et utiliser la fonction copy_labels_from de l’extension labelled.
ATTENTION : il y 11 enfants soi-disant nés après la date d’enquête ! Quelle que soit l’enquête, il est rare de ne pas observer d’incohérences. Dans le cas présent, il est fort possible que la date d’entretien puisse parfois être erronnée (par exemple si l’enquêteur a inscrit une date sur le questionnaire papier le jour du recensement du ménage mais n’ai pu effectué le questionnaire individuel que plus tard). Nous décidons ici de procéder à une correction en ajoutant un mois aux dates d’entretien problématiques. D’autres approches auraient pu être envisagées, comme par exemple exclure ces observations problématiques. Cependant, cela aurait impacté le calcul du range de naissance pour les autres enfants issus de la même mère. Quoiqu’il en soit, il n’y a pas de réponse unique. À vous de vous adapter au contexte particulier de votre analyse.
The following object is masked from 'package:data.table':
+
+ :=
+
Warning: replacing previous import 'ggplot2:::=' by
+'data.table:::=' when loading 'survMisc'
+
ggsurvplot(km_global)
Courbe de survie de Kaplan-Meier
On peut facilement représenter à la place la courbe cumulée des évènements (l’inverse de la courbe de survie) et la table des effectifs en fonction du temps.
-
ggsurvplot(km_global, fun ="event", risk.table =TRUE, surv.scale ="percent")
+
ggsurvplot(km_global, fun ="event", risk.table =TRUE, surv.scale ="percent")
Courbe cumulée des évènements et table des effectifs
Pour comparer deux groupes (ici les filles et les garçons), il suffit d’indiquer la variable de comparaison à survfit.
-
km_sexe <-survfit(Surv(time, deces) ~ sexe, data = enfants)
-km_sexe
+
km_sexe <-survfit(Surv(time, deces) ~ sexe, data = enfants)
+km_sexe
Call: survfit(formula = Surv(time, deces) ~ sexe, data = enfants)
n events median 0.95LCL 0.95UCL
sexe=masculin 762 94 NA NA NA
sexe=féminin 822 48 NA NA NA
La fonction survdiff permets de calculer le test du logrank afin de comparer des courbes de survie. La mortalité infanto-juvénile diffère-t-elle significativement selon le sexe de l’enfant ?
-
survdiff(Surv(time, deces) ~ sexe, data = enfants)
+
survdiff(Surv(time, deces) ~ sexe, data = enfants)
Likelihood ratio test=38.16 on 15 df, p=0.0008546
n= 1584, number of events= 142
De nombreuses variables ne sont pas significatives. Voyons si nous pouvons, avec la fonction step, améliorer notre modèle par minimisation de l’AIC ou Akaike Information Criterion (voir la section Sélection de modèles du chapitre sur la Régression logistique).
On peut obtenir facilement les coefficients du modèle avec l’excellente fonction tidy de l’extension broom. Ne pas oublier de préciser exponentiate = TRUE. En effet, dans le cas d’un modèle de Cox, l’exponentiel des coefficients corresponds au ratio des risques instantannés ou hazard ratio (HR) en anglais.
Pour représenter ces rapports de risque, on peut ici encore avoir recours à la fonction ggcoef_model de l’extension GGally.
-
library(GGally, quietly =TRUE)
+
library(GGally, quietly =TRUE)
Registered S3 method overwritten by 'GGally':
method from
+.gg ggplot2
-
ggcoef_model(mod2, exponentiate =TRUE)
+
ggcoef_model(mod2, exponentiate =TRUE)
Coefficients du modèle avec ggcoef_model
L’extension survminer fournit également une fonction ggforest qui permet de représenter de manière plus esthétique et complète les coefficients d’un modèle de Cox.
-
ggforest(mod2)
+
ggforest(mod2)
Coefficients du modèle avec ggforest
@@ -28499,8 +29828,8 @@
Modèle de Cox
Vérification de la validité du modèle
Un modèle de Cox n’est valable que sous l’hypothèse de la proportionnalité des risques relatifs. Selon cette hypothèse les résidus de Schoenfeld ne dépendent pas du temps. Cette hypothèse peut être testée avec la fonction cox.zph.
-
test <-cox.zph(mod2)
-test
+
test <-cox.zph(mod2)
+test
chisq df p
sexe 0.505 1 0.48
milieu 0.136 1 0.71
@@ -28508,7 +29837,7 @@
Vérification de la validité du modèle
GLOBAL 0.805 3 0.85
Une valeur de p inférieure à 5 % indique que l’hypothèse n’est pas vérifiée. Il apparaît que p est supérieur à 5 % globalement et pour chaque variable prise individuellement. Notre modèle est donc valide.
Il est possible de représenter la distribution des résidus de Schoenfeld à l’aide de ggcoxzph de l’extension survminer, afin de voir si leur répartition change au cours du temps.
-
ggcoxzph(test)
+
ggcoxzph(test)
Warning in regularize.values(x, y, ties, missing(ties),
na.rm = na.rm): collapsing to unique 'x' values
@@ -28536,7 +29865,7 @@
Données pondérées
-
+
@@ -28744,58 +30073,48 @@
Représentations graphiques
Tapis des séquences triés par multidimensional scaling
On voit mieux apparaître ainsi l’hétérogénéité de certaines classes. Les classes 1, 3 et 4, par exemple, semblent regrouper des carrières relativement stables (respectivement de professions intermédiaires, d’employés et de cadres) et des carrières plus « mobiles » commencées comme ouvrier (classes 1 et 3, en orange) ou comme profession intermédiaire (classe 4, en rouge). De même, la majorité des membres de la dernière classe commencent leur carrière dans un groupe professionnel distinct de celui qu’ils occuperont par la suite (indépendants). Ces distinctions apparaissent d’ailleurs si on relance le programme avec un nombre plus élevé de classes (en remplaçant le 5 de la ligne nbcl <- 5 par 7, seconde inflexion de la courbe des sauts d’inertie, et en exécutant de nouveau le programme à partir de cette ligne) : les stables et les mobiles se trouvent alors dans des classes distinctes.
-
Le package JLutils, disponible sur GitHub, propose une fonction seq_heatmap permettant de représenter le tapis de l’ensemble des séquences selon l’ordre du dendrogramme.
-
Pour installer JLutils, on aura recours au package devtools et à sa fonction install_github :
Registered S3 method overwritten by 'GGally':
- method from
- +.gg ggplot2
-
seq_heatmap(seq, seq.dist, labCol =14:50)
-
+
Le package seqhandbook propose une fonction seq_heatmap permettant de représenter le tapis de l’ensemble des séquences selon l’ordre du dendrogramme.
+
library(seqhandbook)
+
+Attaching package: 'seqhandbook'
+
The following object is masked from 'package:JLutils':
+
+ seq_heatmap
+
seq_heatmap(seq, seq.dist, labCol =14:50)
+
Tapis des séquences trié selon le dendrogramme
Il est possible de reproduire un tapis de séquence avec ggplot2. Outre le fait que cela fournit plus d’options de personnalisation du graphique, cela permets également à ce que la hauteur de chaque classe sur le graphique soit proportionnelle aux nombre d’invidus.
En premier lieu, on a va ajouter à notre fichier de données des identifiants individuels, la typologie crée et l’ordre obtenu par multidimensional scaling.
Ensuite, il est impératif que nos données soient dans un format long et tidy, c’est-à-dire avec une ligne par individu et par pas de temps. Pour cela on aura recours à la fonction gather (voir le chapitre dédié).
-
library(tidyr)
-
-Attaching package: 'tidyr'
-
The following object is masked from 'package:magrittr':
-
- extract
-
long <- donnees %>%gather(csp1:csp37, key = annee, value = csp)
On va mettre en forme la variable csp sous forme de facteur, récupérer l’année grace à la fonction str_sub de l’extension stringr (voir le chapitre sur la manipulation de texte) et recalculer l’âge.
Il n’y a plus qu’à faire notre graphique grace à geom_raster qui permet de colorier chaque pixel. Techniquement, pour un tapis de séquence, il s’agit de représenter le temps sur l’axe horizontal et les individus sur l’axe vertical. Petite astuce : plutôt que d’utiliser id pour l’axe vertical, nous utilisons ordre afin de trier les observations. Par ailleurs, il est impératif de transformer au passage ordre en facteur afin que ggplot2 puisse recalculer proprement et séparément les axes pour chaque facette15, à condition de ne pas oublier l’option scales = "free_y" dans l’appel à facet_grid. Les autres commandes ont surtout pour vocation d’améliorer le rendu du graphique (voir le chapitre dédié à ggplot2).
La distance des séquences d’une classe au centre de cette classe, obtenue avec disscenter, permet de mesurer plus précisément l’homogénéité des classes. Nous utilisons ici aggregate pour calculer la moyenne par classe :
-
aggregate(disscenter(as.dist(seq.om), group = seq.part), list(seq.part), mean)
+
aggregate(disscenter(as.dist(seq.om), group = seq.part), list(seq.part), mean)
On poursuit ensuite la description des classes en croisant la typologie avec la variable generation :
Pearson's Chi-squared test
@@ -28916,7 +30235,7 @@
Bibliographie
-
+
@@ -28941,88 +30260,203 @@
Trajectoires de soins : un exemple de données lon
+
+
Une version de ce chapitre ustilisant l’extension data.table plutôt que dplyrest également disponible.
+
Dans ce chapitre, nous allons aborder plusieurs méthodes d’analyse à partir d’un jeu de données longitudinales. Tout d’abord, importons les données dans R avec la commande suivante :
Nous obtenons un objet appelé care_trajectories. La fonction class nous montre qu’il s’agit d’un tableau de données au format data.table (voir le chapitre dédié). Chargeons donc cette extension ainsi que le tidyverse.
+
Nous obtenons un objet appelé care_trajectories. La fonction class nous montre qu’il s’agit d’un tableau de données au format data.table (voir le chapitre dédié). Chargeons le tidyverse et utilisons as_tibble pour convertir le tableau de données.
[1] less than 10 km 26804 54.3
[2] 10 km or more 22561 45.7
Total 49365 100.0
+
Nous pouvons également utiliser look_for de labelled.
+
care_trajectories %>%look_for()
+
+
+
+pos
+
+
+variable
+
+
+label
+
+
+col_type
+
+
+levels
+
+
+value_labels
+
+
+
+
+
+1
+
+
+id
+
+
+patient identifier
+
+
+int
+
+
+NULL
+
+
+NULL
+
+
+
+
+2
+
+
+month
+
+
+month(s) since diagnosis
+
+
+dbl
+
+
+NULL
+
+
+NULL
+
+
+
+
+3
+
+
+care_status
+
+
+care status
+
+
+chr+lbl
+
+
+NULL
+
+
+D, C, T, S
+
+
+
+
+4
+
+
+sex
+
+
+sex
+
+
+dbl+lbl
+
+
+NULL
+
+
+0, 1
+
+
+
+
+5
+
+
+age
+
+
+age group
+
+
+dbl+lbl
+
+
+NULL
+
+
+1, 2, 3
+
+
+
+
+6
+
+
+education
+
+
+education level
+
+
+dbl+lbl
+
+
+NULL
+
+
+1, 2, 3
+
+
+
+
+7
+
+
+wealth
+
+
+wealth group (assets score)
+
+
+dbl+lbl
+
+
+NULL
+
+
+1, 2, 3
+
+
+
+
+8
+
+
+distance_clinic
+
+
+distance to nearest clinic
+
+
+dbl+lbl
+
+
+NULL
+
+
+1, 2
+
+
+
+
Dans cette étude, on a suivi des patients à partir du moment où ils ont été diagnostiqués pour une pathologie grave et chronique et on a suivi leurs parcours de soins chaque mois à partir du diagnostic. La variable status contient le statut dans les soins de chaque individu pour chaque mois de suivi :
@@ -29119,66 +30739,56 @@
Première description des données
Il est important de noter que nous avons ici des statuts hiérarchiquement ordonnés (D < C < T < S), ce qui aura son importance pour les choix méthodologiques que nous aurons à faire.
Nous disposons également d’autres variables (âge, sexe, niveau d’éducation…) qui sont ici dépendantes du temps, c’est-à-dire que le cas échéant, elles peuvent varier d’un mois à l’autre en cas de changement.
Avant de démarrer les analyses, françisons certaines de ces variables.
Le fichier contient 49 365 lignes, ce qui ne veut pas dire qu’il y a ce nombre d’invidus suivis au cours du temps, puisque plusieurs lignes correspondent à un même individu. On peut obtenir le nombre d’individus différents assez facilement avec la commande :
Le fichier contient 49365 lignes, ce qui ne veut pas dire qu’il y a ce nombre d’invidus suivis au cours du temps, puisque plusieurs lignes correspondent à un même individu. On peut obtenir le nombre d’individus différents assez facilement avec la commande :
+
length(unique(care_trajectories$id))
[1] 2929
Précision : dans ce fichier, tous les individus ne sont pas suivis pendant la même durée, car ils n’ont pas tous été diagnostiqués au même moment. Cependant, il n’y a pas de trous dans le suivi (ce qui serait le cas si certains individus sortaient de l’observation pendant quelques mois puis re-rentraient dans la cohorte de suivi).
Avant d’aller plus avant, il nous faut avoir une idée du nombre d’individus observé au cours du temps, ce que l’on peut obtenir avec :
Améliorons ce graphique en y ajoutant la distribution selon le statut dans les soins chaque mois, en améliorant l’axe du temps (tous les 6 mois est plus facile à lire) et en y ajoutant un titre et des étiquettes appropriées. Afin de disposer d’une palette de couleurs à fort contraste, nous allons utiliser l’extension viridis. Enfin, nous allons utiliser une petite astuce pour indiquer les effectifs sur l’axe horizontal. Au passage, nous allons également franciser les étiquettes de la variable care_status avec val_labels (notez aussi le recours à to_factor dans aes qui nous permet de transformer à la volée la variable en facteur, format attendu par ggplot2 pour les variables catégorielles). On se référera au chapitre dédié à ggplot2 pour plus de détails sur les différentes fonctions de cette extension graphique.
-
library(viridis)
-n <- care_trajectories[month %in% (0:8*6), .(n = .N), by = month]$n
-etiquettes <-paste0("M", 0:8*6, "\n(n=", n, ")")
-val_labels(care_trajectories$care_status) <-c(
-"diagnostiqué, mais pas suivi"="D",
-"suivi, mais pas sous traitement"="C",
-"sous traitement, mais infection non contrôlée"="T",
-"sous traitement et infection contrôlée"="S"
-)
-ggplot(care_trajectories) +
-aes(x = month, fill =to_factor(care_status)) +
-geom_bar(color ="gray50", width =1) +
-scale_x_continuous(breaks =0:8*6, labels = etiquettes) +
-ggtitle("Distribution du statut dans les soins chaque mois") +
-xlab("") +ylab("") +
-theme_light() +
-theme(legend.position ="bottom") +
-labs(fill ="Statut dans les soins") +
-scale_fill_viridis(discrete =TRUE, direction =-1) +
-guides(fill =guide_legend(nrow =2))
Améliorons ce graphique en y ajoutant la distribution selon le statut dans les soins chaque mois, en améliorant l’axe du temps (tous les 6 mois est plus facile à lire) et en y ajoutant un titre et des étiquettes appropriées. Afin de disposer d’une palette de couleurs à fort contraste, nous allons utiliser l’extension viridis. Enfin, nous allons utiliser une petite astuce pour indiquer les effectifs sur l’axe horizontal. Au passage, nous allons également franciser les étiquettes de la variable care_status avec val_labels (notez aussi le recours à to_factor dans aes qui nous permet de transformer à la volée la variable en facteur, format attendu par ggplot2 pour les variables catégorielles). On se référera au chapitre dédié à ggplot2 pour plus de détails sur les différentes fonctions de cette extension graphique.
+
library(viridis)
+n <- care_trajectories %>%
+filter(month %in% (0:8*6)) %>%
+group_by(month) %>%
+count() %>%
+pluck("n")
+etiquettes <-paste0("M", 0:8*6, "\n(n=", n, ")")
+val_labels(care_trajectories$care_status) <-c(
+"diagnostiqué, mais pas suivi"="D",
+"suivi, mais pas sous traitement"="C",
+"sous traitement, mais infection non contrôlée"="T",
+"sous traitement et infection contrôlée"="S"
+)
+ggplot(care_trajectories) +
+aes(x = month, fill =to_factor(care_status)) +
+geom_bar(color ="gray50", width =1) +
+scale_x_continuous(breaks =0:8*6, labels = etiquettes) +
+ggtitle("Distribution du statut dans les soins chaque mois") +
+xlab("") +ylab("") +
+theme_light() +
+theme(legend.position ="bottom") +
+labs(fill ="Statut dans les soins") +
+scale_fill_viridis(discrete =TRUE, direction =-1) +
+guides(fill =guide_legend(nrow =2))
+
On s’aperçoit qu’une majorité des personnes suivies ne l’ont été que peu de temps, avec une décroissance rapide des effectifs.
Évolution de la cascade de soins au cours du temps
On nomme communément cascade de soins la proportion d’individus dans chaque statut à un moment du temps donné. On peut facilement obtenir celle-ci à partir du code du graphique précédent en ajoutant l’option position = fill à geom_bar.
-
ggplot(care_trajectories) +
-aes(x = month, fill =to_factor(care_status)) +
-geom_bar(color ="gray50", width =1, position ="fill") +
-scale_x_continuous(breaks =0:8*6, labels = etiquettes) +
-scale_y_continuous(labels = scales::percent) +
-ggtitle("Cascade des soins observée, selon le temps depuis le diagnostic") +
-xlab("") +ylab("") +
-theme_light() +
-theme(legend.position ="bottom") +
-labs(fill ="Statut dans les soins") +
-scale_fill_viridis(discrete =TRUE, direction =-1) +
-guides(fill =guide_legend(nrow =2))
-
-
Les effectifs sont très faibles au-delà de 36 mois et il serait préférable de couper la cascade au-delà de M36, ce que l’on peut faire aisément ne gardant que les lignes correspondantes de care_trajectories.
Évolution de la cascade de soins au cours du temps
theme(legend.position ="bottom") +labs(fill ="Statut dans les soins") +scale_fill_viridis(discrete =TRUE, direction =-1) +
-guides(fill =guide_legend(nrow =2))
-casc_obs
+guides(fill =guide_legend(nrow =2))
+
Les effectifs sont très faibles au-delà de 36 mois et il serait préférable de couper la cascade au-delà de M36, ce que l’on peut faire aisément ne gardant que les lignes correspondantes de care_trajectories.
+
casc_obs <-ggplot(care_trajectories %>%filter(month <=36)) +
+aes(x = month, fill =to_factor(care_status)) +
+geom_bar(color ="gray50", width =1, position ="fill") +
+scale_x_continuous(breaks =0:8*6, labels = etiquettes) +
+scale_y_continuous(labels = scales::percent) +
+ggtitle("Cascade des soins observée, selon le temps depuis le diagnostic") +
+xlab("") +ylab("") +
+theme_light() +
+theme(legend.position ="bottom") +
+labs(fill ="Statut dans les soins") +
+scale_fill_viridis(discrete =TRUE, direction =-1) +
+guides(fill =guide_legend(nrow =2))
+casc_obs
+
Analyse de survie classique
@@ -29204,106 +30828,70 @@
Analyse de survie classique
En toute rigueur, il faudrait également considérer les transitions inverses (sortie de soins, arrêt du traitement, échec virologique). De même, il est possible que certains aient vécus plusieurs transitions successives (entrée en soin, initiation du traitement, sortie de soins et arrêt du traitement, nouvelle entrée en soins…).
Pour le moment, contentons-nous de regarder la première entrée en soins, la première initiation du traitement et la première atteinte du contrôle de l’infection et de calculer la date de ces trois événements, dans un fichier ind contenant une ligne par individu.
Il nous faut également la durée de suivi par individu.
-
ind <-merge(
- ind,
- care_trajectories[, .(suivi =max(month)), by = id],
-by ="id",
-all.x =TRUE
-)
+
ind <- ind %>%
+left_join(
+ care_trajectories %>%
+group_by(id) %>%
+ dplyr::summarise(suivi =max(month)),
+by ="id"
+ )
Pour faciliter la suite des analyses, nous allons nous créer une petite fonction qui, en fonction de la date d’origine et la date d’événement retenues, calculera la courbe de Kaplan-Meier correspondante (voir le chapitre sur l’analyse de suivie pour le calcul de la courbe de survie et celui dédié à l’écriture de fonctions).
-
km <-function(date_origine, date_evenement, nom) {
-library(survival)
-
-# ne garder que les observations avec date d'origine
- tmp <- ind[!is.na(ind[[date_origine]]), ]
-
-# pre-remplir la variable time avec duree de suivi
-# depuis date d'origine
- tmp$time <- tmp$suivi - tmp[[date_origine]]
-# et considérer que l'événement n'a pas été vécu
- tmp$event <-FALSE
-
-# si date_evement documentée, événement vécu
- tmp[!is.na(tmp[[date_evenement]]), ]$event <-TRUE
- tmp[tmp$event ==TRUE, ]$time <-
- tmp[tmp$event ==TRUE, ][[date_evenement]] -
- tmp[tmp$event ==TRUE, ][[date_origine]]
-
- kaplan <-survfit(Surv(time, event) ~1, data = tmp)
- res <- broom::tidy(kaplan, conf.int =TRUE)
- res$nom <- nom
- res
-}
+
km <-function(date_origine, date_evenement, nom) {
+library(survival)
+
+# ne garder que les observations avec date d'origine
+ tmp <- ind[!is.na(ind[[date_origine]]), ]
+
+# pre-remplir la variable time avec duree de suivi
+# depuis date d'origine
+ tmp$time <- tmp$suivi - tmp[[date_origine]]
+# et considérer que l'événement n'a pas été vécu
+ tmp$event <-FALSE
+
+# si date_evement documentée, événement vécu
+ tmp[!is.na(tmp[[date_evenement]]), ]$event <-TRUE
+ tmp[tmp$event ==TRUE, ]$time <-
+ tmp[tmp$event ==TRUE, ][[date_evenement]] -
+ tmp[tmp$event ==TRUE, ][[date_origine]]
+
+ kaplan <-survfit(Surv(time, event) ~1, data = tmp)
+ res <- broom::tidy(kaplan, conf.int =TRUE)
+ res$nom <- nom
+ res
+}
Une première approche consiste à regarder la survenue de chacun des trois événements mentions plus haut en fonction du temps depuis le diagnostic.
-
depuis_diag <- dplyr::bind_rows(
-km("diagnostic", "entree_soins", "entrée en soins"),
-km("diagnostic", "initiation_tt", "initiation du traitement"),
-km("diagnostic", "controle", "contrôle de l'infection")
-)
-
-g_diag <-ggplot(data = depuis_diag) +
-aes(x = time, y =1- estimate,
-color =as_factor(nom), fill =as_factor(nom),
-ymin =1- conf.high, ymax =1- conf.low) +
-geom_ribbon(alpha = .25, mapping =aes(color =NULL)) +
-geom_line(size =1) +
-theme_classic() +
-theme(
-legend.position ="bottom",
-panel.grid.major.y =element_line(colour ="grey")
- ) +
-scale_x_continuous(breaks =0:6*6, limits =c(0, 36)) +
-scale_y_continuous(labels = scales::percent, limits =c(0, 1)) +
-xlab("mois depuis le diagnostic") +
-ylab("") +labs(color ="", fill ="")
-g_diag
-
-
Ce graphique ressemble à la cascade des soins observée que nous avions calculée plus haut, à quelques différences près :
-
-
avec la méthode de Kaplan-Meier, la censure à droite, i.e. le fait que tous les individus n’ont pas la même durée de suivi, est correctement prise en compte et la courbe est corrigée en conséquence ;
-
par contre, les transitions inverses ne sont pas considérées : lorsqu’un individu a atteint une étape, on ne regarde pas s’il en ressort.
-
-
Une autre manière d’appréhender nos trajectoires est de considérer le temps requis pour atteindre une étape une fois la précédente étape atteinte. Ce qu’on obtient facilement en adaptant légèrement notre code précédent.
-
depuis_prec <- dplyr::bind_rows(
+
depuis_diag <- dplyr::bind_rows(km("diagnostic", "entree_soins", "entrée en soins"),
-km("entree_soins", "initiation_tt", "initiation du traitement"),
-km("initiation_tt", "controle", "contrôle de l'infection")
+km("diagnostic", "initiation_tt", "initiation du traitement"),
+km("diagnostic", "controle", "contrôle de l'infection"))
-g_prec <-ggplot(data = depuis_prec) +
+g_diag <-ggplot(data = depuis_diag) +aes(x = time, y =1- estimate, color =as_factor(nom), fill =as_factor(nom),ymin =1- conf.high, ymax =1- conf.low) +
@@ -29316,14 +30904,43 @@
Analyse de survie classique
) +scale_x_continuous(breaks =0:6*6, limits =c(0, 36)) +scale_y_continuous(labels = scales::percent, limits =c(0, 1)) +
-xlab("mois depuis l'étape précédente") +
+xlab("mois depuis le diagnostic") +ylab("") +labs(color ="", fill ="")
-
-g_prec
+g_diag
-
Attention : cette représentation graphique peut éventuellement prêter à confusion dans sa lecture car l’échelle de temps n’est pas tout à fait la même pour chaque courbe, dans la mesure où la date d’origine diffère pour chacune. Dès lors, il peut être plus pertinent de présenter chaque courbe l’une à côté de l’autre.
-
g_prec +facet_grid(~as_factor(nom))
+
Ce graphique ressemble à la cascade des soins observée que nous avions calculée plus haut, à quelques différences près :
+
+
avec la méthode de Kaplan-Meier, la censure à droite, i.e. le fait que tous les individus n’ont pas la même durée de suivi, est correctement prise en compte et la courbe est corrigée en conséquence ;
+
par contre, les transitions inverses ne sont pas considérées : lorsqu’un individu a atteint une étape, on ne regarde pas s’il en ressort.
+
+
Une autre manière d’appréhender nos trajectoires est de considérer le temps requis pour atteindre une étape une fois la précédente étape atteinte. Ce qu’on obtient facilement en adaptant légèrement notre code précédent.
Attention : cette représentation graphique peut éventuellement prêter à confusion dans sa lecture car l’échelle de temps n’est pas tout à fait la même pour chaque courbe, dans la mesure où la date d’origine diffère pour chacune. Dès lors, il peut être plus pertinent de présenter chaque courbe l’une à côté de l’autre.
+
g_prec +facet_grid(~as_factor(nom))
+
Ici, on voit plus clairement que l’étape où il y a le plus de perdition est celle de l’entrée en soins, moins des trois quarts des personnes diagnostiquées étant venu en clinique au mois une fois dans les trois ans suivant le diagnostic. Par contre, l’initiation du traitement une fois entré en clinique et le contrôle de l’infection une fois le traitement initié sont beaucoup plus rapide.
Pour aller plus loin avec les outils de l’analyse de survie classique, il serait possible de faire des analyses bivariées (Kaplan-Meier) ou multivariées (Cox) pour chacune de ces étapes. Cependant, il serait plus intéressant de trouver une approache statistique permettant de considérer dans un même modèle l’ensemble des transitions possibles.
@@ -29331,407 +30948,1130 @@
Analyse de survie classique
Une première analyse de séquences sur l’ensemble du fichier
L’analyse de séquences permet d’appréhender l’ensemble de la trajectoire de soins à travers la succession des états dans lesquels se trouvent les patients observés.
Nous allons donc réaliser une analyse de séquences (voir le chapitre dédié) sur l’ensemble de notre fichier. Pour cela, il va falloir préalable que nous transformions nos donnée actuellement dans un format long en un tableau large, c’est-à-dire avec une ligne par individu et une variable différentes par pas de temps. On peut réaliser cela facilement avec pivot_wider de tidyr (voir le chapitre dédié à tidyr). Il faut noter que le résultat ne sera pas un tableau data.table, d’où le recours à setDT pour convertir le résultat.
On utilise seqdef de TraMineR pour créer nos séquences, avec les arguments alphabet pour forcer l’ordre de l’alphabet, states pour spécifier des étiquettes courtes à chaque état et cpal pour indiquer le code couleur de chaque état (et être raccord avec nos graphiques précédents).
-
library(TraMineR, quietly =TRUE)
-
-TraMineR stable version 2.2-1 (Built: 2020-11-02)
-
Website: http://traminer.unige.ch
-
Please type 'citation("TraMineR")' for citation information.
library(TraMineR, quietly =TRUE)
+seq_all <-seqdef(
+ large %>% dplyr::select(m0:m50),
+id = large$id,
+alphabet =c("D", "C", "T", "S"),
+states =c("diagnostiqué", "en soins", "sous traitement", "inf. contrôlée"),
+cpal =viridis(4, direction =-1)
+)
[>] found missing values ('NA') in sequence data
[>] preparing 2929 sequences
[>] coding void elements with '%' and missing values with '*'
@@ -29744,207 +32084,556 @@
Une première analyse de séquences sur l’ensemble du fichier
[>] 2929 sequences in the data set
[>] min/max sequence length: 1/51
On peut retrouver la cascade de soins avec seqdplot.
-
seqdplot(seq_all, legend.prop = .25)
-
+
seqdplot(seq_all, legend.prop = .25)
+
Nous allons maintenant calculer une matrice des distances entre individus par optimal matching. Dans le cas présent, nos différents status sont hiérarchiquement ordonnés. Il n’est donc pas raisonnable de penser que les coûts sont constants entre les différents statuts, puisqu’en un sens, passer directement de D à T peut être considéré comme être passé d’abord de D à C puis de C à D. Nous allons donc faire une matrice de coûts hiérarchisée. seqcost nous permets de produire une matrice de coûts constants, que nous allons ensuite modifier manuellement. Pour le coût indel, le plus simple est de considérer la moitié du coût de substitution maximum.
-
couts <-seqcost(seq_all, method ="CONSTANT")
+
couts <-seqcost(seq_all, method ="CONSTANT")
[>] creating 4x4 substitution-cost matrix using 2 as constant value
The following object is masked from 'package:JLutils':
+
+ seq_heatmap
+
seq_heatmap(seq_all, arbre_all)
+
+
Il apparaît que les différentes séquences sont principalement regroupées en fonction de leur longueur. En effet, pour passer d’une séquence courte à une séquence longue il faut ajouter des statuts pour compléter la séquence ce qui induit de facto une distance élevée (en raison du coût indel). Dès lors, lorsque l’on travaille avec des séquences aux longueurs très disparates, une classification ascendante hiérarchique va produire une typologie de séquences courtes et de séquences longues, ce qui n’est pas forcément ce que l’on recherche.
+
Dans notre exemple, nous pouvons considérer que les séquences courtes ne sont pas pertinentes à retenir dans l’analyse car l’observation n’est pas assez longue pour voir le parcours de soins des patients. Une solution consiste à ne retenir que les individus observées au moins n mois et analyser leur trajectoire sur seulement n mois, ce qui permet de n’avoir que des séquences de même longueur. Dès lors, la distance entre deux séquences ne dépendra plus que des différences de parcours. On serait tenté de prendre un n élévé pour avoir ainsi des parcours de soins longs. Mais dans ce cas là, l’analyse ne se fera que sur un tout petit nombre d’individus et on manquera de puissance. Si, à l’inverse, on prends un n petit, nous aurons des effectifs élevés mais les séquences seront peut-être trop courtes pour mettre en évidence la variété des trajectoires. Il faut dès lors trouver un compromis entre ces deux contraintes.
+
Si l’on regarde notre premier graphique montrant le nombre d’observations au cours du temps, il apparaît une sorte de point d’inflexion au niveau de M18 avec un brusque décrochage. D’un autre côté, 18 mois offre un minimum de durée d’observations pour espérer voir émerger des trajectoires plus complexes.
+
+
+
Une seconde analyse de séquences limitées aux 18 premiers mois
+
Reprenons notre analyse en se limitant aux individus observés au moins 18 mois (soit 19 status entre M0 et M18) et en se limitant aux 18 premiers mois pour modéliser les séquences. La fonction seqlength permets de récupérer la longueur de chaque séquence.
Reste maintenant à décider du nombre de classes à retenir. Encore une fois, c’est un équilibre à trouver entre le niveau de détails voulus et le niveau de simplification requis pour permettre l’analyse.
+
Pour faciliter ce choix, on peut avoir recours à la fonction as.seqtree de l’extension WeightedCluster, couplée à la fonction seqtreedisplay. ATTENTION : pour que le graphique puisse être produit, il faut que le logiciel libre GraphViz (https://graphviz.gitlab.io/) soit installé sur votre PC. On peut également installer GraphViz avec le code ci-dessous :
+
if (!requireNamespace("BiocManager", quietly =TRUE))
+install.packages("BiocManager")
+BiocManager::install("Rgraphviz")
+
La combinaison des deux fonctions va permettre de représenter l’évolution des catégories au fur-et-à-mesure que l’on coupe le dendrogramme plus bas. On peut choisir le type de graphique utilisé avec l’argument type (voir l’aide de seqplot) et le nombre maximum de clusters avec nclust.
Représentation graphique produite par seqtreedisplay
+
+
Afin d’éviter de multiplier les sous-groupes, nous n’allons conserver que 4 catégories.
+
large_m18$typo_cah <-cutree(arbre_m18, 4)
+
+
Il est aussi possible de se baser sur divers indicateurs statistiques sur la qualité de chaque partition. Pour cela, on pourra par exemple avoir recours à la fonction as.clustrange de l’extension WeightedCluster. Essayez par exemple les commandes ci-après :
On peut représenter le découpage du dendrogramme avec fviz_dend fournie par factoextra
+
library(factoextra)
+
Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
+
fviz_dend(arbre_m18, k =4, show_labels =FALSE, rect =TRUE)
+
+
Comme expliqué par Matthias Studer dans le manuel de la librairie WeightedCluster, plusieurs critiques peuvent être adressées aux procédures hiérarchiques, en particulier le fait que la fusion de deux groupes se fait en maximisant un critère local.
+
L’algorithme PAM pour Partitioning Around Medoids suit une autre logique que les algorithmes hiérarchiques et vise à obtenir la meilleure partition d’un ensemble de données en un nombre prédéfini de groupes. Il a l’avantage de maximiser un critère global et non uniquement un critère local. Par contre, le nombre de classes doit être fixé à l’avance.
+
Ayant décidé de retenir 4 classes au regard de notre classification ascendante hiérarchique, nous pouvons voir si l’algorithme PAM permets d’améliorer nos 4 classes. Nous allons utiliser la fonction wcKMedoids de l’extension WeightedCluster en lui indiquant comme partition initiale celle obtenue avec la classigication hiérarchique.
+
pam_m18 <-wcKMedoids(dist_m18, k =4, initialclust = arbre_m18)
+large_m18$typo_pam <- pam_m18$clustering
+
Un tableau croisé nous permets de voir que les deux typologies restent proches.
+
library(gtsummary)
+large_m18 %>%
+tbl_cross(row = typo_cah, col = typo_pam)
+
+
+
+
+
Characteristic
+
+ typo_pam
+
+
Total
+
+
+
6
+
23
+
85
+
410
+
+
+
+
+
typo_cah
+
+
+
+
+
+
+
+
1
+
522
+
4
+
0
+
39
+
565
+
+
+
2
+
0
+
3
+
330
+
3
+
336
+
+
+
3
+
1
+
203
+
59
+
13
+
276
+
+
+
4
+
23
+
0
+
5
+
112
+
140
+
+
+
Total
+
546
+
210
+
394
+
167
+
1,317
+
-
+
Regardons les tapis de séquence des deux typologies.
-
large_m18$ordre_cmd <-cmdscale(as.dist(dist_m18), k =1)
-seqIplot(seq_m18, group = large_m18$typo_cah, sortv = large_m18$ordre_cmd)
-
-
seqIplot(seq_m18, group = large_m18$typo_pam, sortv = large_m18$ordre_cmd)
-
+
large_m18$ordre_cmd <-cmdscale(as.dist(dist_m18), k =1)
+seqIplot(seq_m18, group = large_m18$typo_cah, sortv = large_m18$ordre_cmd)
+
+
seqIplot(seq_m18, group = large_m18$typo_pam, sortv = large_m18$ordre_cmd)
+
Comme on le voit les deux typologies obtenues sont très proches. Suivant le cas, à vous de choisir celle qui semble la plus pertinente d’un point de vue sociologique. Il existe également divers indicateurs statisques pour mesurer la qualité d’une partition (voir le manuel de la librairie WeightedCluster de Matthias Studer). Ils peuvent être calculés avec la fonction wcClusterQuality. Comparons les deux typologies obtenues.
Une seconde analyse de séquences limitées aux 18 premiers mois
Selon ces indicateurs calculés, l’approche PAM obtiendrait une partition légèrement de meilleure qualité que celle obtenuepar CAH.
L’extension WeightedCluster fournie aussi une fonction wcSilhouetteObs permettant de mesurer la silhouette de chaque séquence. Plus cette métrique est élevée et proche de 1, plus la séquence est proche du centre de classe et caractéristique de la classe. On peut utiliser cette métrique pour classer les séquences sur le tapis de séquences.
-
large_m18$sil <-wcSilhouetteObs(dist_m18, large_m18$typo_pam)
-seqIplot(seq_m18, group = large_m18$typo_pam, sortv = large_m18$sil)
-
+
large_m18$sil <-wcSilhouetteObs(dist_m18, large_m18$typo_pam)
+seqIplot(seq_m18, group = large_m18$typo_pam, sortv = large_m18$sil)
+
Nous voyons émerger quatre groupes distincts :
les rapides, qui entrent en soins et initient le traitement dès les premiers mois suivant le diagnostic ;
@@ -30348,75 +33037,78 @@
Une seconde analyse de séquences limitées aux 18 premiers mois
les hors soins, qui ne sont pas entrés en soins où n’y sont pas restés.
Le graphique obtenu avec seqIplot affiche visuellement chaque groupe avec la même hauteur, pouvant laisser accroire que chaque groupe a le même poids dans l’échantillon. Pour produire une représentation graphique des tapis de séquences plus correcte, où chaque la hauteur de chaque groupe correspondrait à son poids dans l’échantillon, nous allons passer par ggplot2. Un tapis de séquences peut-être vu comme un raster et dès lors représenté avec geom_raster. Pour travailler avec ggplot2, nos données doivent être au format tidy, c’est-à-dire avec une ligne par point d’observation, soit une ligne par personne et par jour. Nous allons donc repartir du fichier care_trajectories. Le mois d’observation indiquera la position en abscisse. Quant à la position en ordonnée, il faudra que nous la calculions, séparément pour chaque groupe, afin d’éviter des lignes vides dans le graphique.
-
# nommer les groupes
-large_m18$groupe <-factor(
- large_m18$typo_pam,
-c(85, 23, 410, 6),
-c("Rapides", "Lents", "Inaboutis", "Hors soins")
-)
-
-# calculer le rang des individus dans chaque groupe
-setorder(large_m18, "ordre_cmd")
-large_m18[, rang_cmd :=1:.N, by = groupe]
-
-# créer un fichier long
-long_m18 <- care_trajectories[id %in% large_m18$id & month <=18]
-long_m18 <-merge(
- long_m18,
- large_m18[, .(id, groupe, rang_cmd)],
-by ="id",
-all.x =TRUE
-)
-long_m18$care_statusF <-to_factor(long_m18$care_status)
-
-# calculer les effectifs par groupe
-tmp <- large_m18[, .(n = .N), by = groupe]
-tmp[, groupe_n :=paste0(groupe, "\n(n=", n, ")")]
-long_m18 <-merge(
- long_m18,
- tmp[, .(groupe, groupe_n)],
-by ="groupe",
-all.x =TRUE
-)
-
# graphique des tapis de séquences
-ggplot(long_m18) +
-aes(x = month, y = rang_cmd, fill = care_statusF) +
-geom_raster() +
-facet_grid(groupe_n ~ ., space ="free", scales ="free") +
-scale_x_continuous(breaks =0:6*3, labels =paste0("M", 0:6*3)) +
-scale_y_continuous(breaks =0:5*100, minor_breaks =NULL) +
-xlab("") +ylab("") +
-theme_light() +
-theme(legend.position ="bottom") +
-labs(fill ="Statut dans les soins") +
-scale_fill_viridis(discrete =TRUE, direction =-1) +
-guides(fill =guide_legend(nrow =2))
-
+
# nommer les groupes
+large_m18$groupe <-factor(
+ large_m18$typo_pam,
+c(85, 23, 410, 6),
+c("Rapides", "Lents", "Inaboutis", "Hors soins")
+)
+
+# calculer le rang des individus dans chaque groupe
+large_m18 <- large_m18 %>%
+group_by(groupe) %>%
+arrange(ordre_cmd) %>%
+mutate(rang_cmd =rank(ordre_cmd, ties.method ="first"))
+
+# créer un fichier long
+long_m18 <- care_trajectories %>%
+filter(id %in% large_m18$id & month <=18) %>%
+left_join(
+ large_m18 %>% dplyr::select(id, groupe, rang_cmd),
+by ="id"
+ )
+
+long_m18$care_statusF <-to_factor(long_m18$care_status)
+
+# calculer les effectifs par groupe
+tmp <- large_m18 %>%
+ dplyr::count(groupe) %>%
+mutate(groupe_n =paste0(groupe, "\n(n=", n, ")"))
+
+long_m18 <- long_m18 %>%
+left_join(
+ tmp %>% dplyr::select(groupe, groupe_n),
+by ="groupe"
+ )
+
# graphique des tapis de séquences
+ggplot(long_m18) +
+aes(x = month, y = rang_cmd, fill = care_statusF) +
+geom_raster() +
+facet_grid(groupe_n ~ ., space ="free", scales ="free") +
+scale_x_continuous(breaks =0:6*3, labels =paste0("M", 0:6*3)) +
+scale_y_continuous(breaks =0:5*100, minor_breaks =NULL) +
+xlab("") +ylab("") +
+theme_light() +
+theme(legend.position ="bottom") +
+labs(fill ="Statut dans les soins") +
+scale_fill_viridis(discrete =TRUE, direction =-1) +
+guides(fill =guide_legend(nrow =2))
+
Facteurs associés à l’appartenance à chaque groupe
Une fois les différents groupes de trajectoires identifiés, il est courant de vouloir regarder si certains facteurs influencent l’appartenance à un groupe plutôt qu’un autre. Dans nos données d’exemple, nous nous intéresserons aux variables suivantes : sexe, groupe d’âges et niveau d’éducation.
Ces différentes variables sont renseignées pour chaque mois dans le tableau de données care_trajectories, en tenant compte des éventuels changements au cours du temps. Ici, notre analyse est menée au niveau individuel. Nous allons donc récupérer la valeur de ces différentes variables au moment du diagnostic, à savoir à M0.
La fonction tbl_summary de l’extension gtsummary permets de produire aisément une série de tableaux croisés. À noter le recours à add_p() pour ajouter les p-valeurs au test du Chi².1 La fonction add_overall permets d’ajouter une colonne avec l’ensemble de l’échantillon. Notez que nous avons utiliser unlabelled de l’extension labelled pour convertir en amont les vecteurs labellisés en facteurs.
Facteurs associés à l’appartenance à chaque groupe
p-value2
-Total, N = 1 3171
+Overall, N = 1,3171
@@ -30777,7 +33469,7 @@
Facteurs associés à l’appartenance à chaque groupe
-
<0,001
+
<0.001
@@ -30804,7 +33496,7 @@
Facteurs associés à l’appartenance à chaque groupe
-
<0,001
+
<0.001
@@ -30827,12 +33519,12 @@
Facteurs associés à l’appartenance à chaque groupe
60+
-
29 (7,4%)
-
11 (5,2%)
-
14 (8,4%)
-
30 (5,5%)
+
29 (7.4%)
+
11 (5.2%)
+
14 (8.4%)
+
30 (5.5%)
-
84 (6,4%)
+
84 (6.4%)
Education
@@ -30840,7 +33532,7 @@
Facteurs associés à l’appartenance à chaque groupe
-
0,002
+
0.002
@@ -30877,114 +33569,118 @@
Facteurs associés à l’appartenance à chaque groupe
1
- Statistique présentée: n (%)
+ n (%)
2
- Test statistique réalisé: test du khi-deux d'indépendance
+ Pearson's Chi-squared test
Une manière de présenter ces mêmes données de manière plus visuelle consiste à réaliser un diagramme en barres cumulées (voir le chapitre sur les graphiques bivariés). Ici, nous utilisons une boucle for (voir le chapitre sur les structures conditionnelles) pour calculer les différents tableaux croisés et les fusionner avec bind_rows (voir la section concaténation de tables du chapitre dédié à dplyr). Nous en profitions également pour calculer le test du Chi² et l’afficher avec le nom de la variable sur le graphique.
Note : pour afficher les proportions sur le graphique, aurons recours recours à la statistique stat_prop de l’extension GGally.
Pour mieux visualiser les relations entre les variables, on peut avoir recours à la fonction ggtable de l’extension GGally qui permets de représenter les résidus du Chi².
Pour mieux visualiser les relations entre les variables, on peut avoir recours à la fonction ggtable de l’extension GGally qui permets de représenter les résidus du Chi².
On considère qu’une cellule est surreprésentée si son résidu standardisé du Chi² est supérieur à 2 ou 3, et qu’elle est sous-représentée si le résidu est inférieur à -2 ou -3.
On peut ainsi noter que les hommes, les jeunes et ceux vivant à plus de 10 kilomètres d’une clinique sont plus souvent dans le groupe Hors soins. Inversement, les femmes et les plus éduqués sont plus souvent dans le groupe des Rapides. Résultat inattendu, les ménages les plus riches sont moins souvent dans les groupes ayant initiés un traitement (Rapides et Lents), ce résultat pouvant s’expliquer en partie par le fait que les plus aisés peuvent plus facilement accéder à des soins dans le secteur privé, et donc faussement apparaître Hors soins, car seuls les soins reçus dans le secteur public ont été mesurés dans cette étude.
Pour affiner les résultats, on aura recours à un modèle multivarié en exécutant une régression logistique multinomiale avec multinom de l’extension nnet (pour plus de détails, voir le chapitre dédié).
-
library(nnet)
-large_m18$groupe2 <-relevel(large_m18$groupe, "Hors soins")
-regm <-multinom(
- groupe2 ~ sex + age + education,
-data =to_factor(large_m18)
-)
Un autre type d’approche envisageable pour identifier des classes de trajectoires est celle des modèles mixtes à classes latentes. Ce type de modèles peut prendre en compte une grande variété d’indicateurs, continus, binaires ou ordinaux. On peut y intégrer des co-variables et il n’est pas nécessaire de disposer du même nombre d’observations par individu.
Nous n’aborderons que brièvement ici ce type de modèles complexes. Sous R, ils peuvent être réalisés via l’extension lcmm et sa fonction homonyme lcmm.
Ici, nous allons modéliser le statut dans les soins en fonction du temps. Il faut indiquer au modèle, via le paramètre ng le nombre de groupes ou classes latentes souhaité. Ici, nous avons retenu 4 en lien avec les résultats de notre analyse de séquences. L’argument link = "thresholds" permets d’indiquer que notre variable d’intérêt est ordinale. Les modèles lcmm peuvent également prendre en compte des variables continues.
Attention : le temps de calcul de ce type de modèle peut être long (plusieurs heures dans notre exemple), suivant le nombre de paramètres, le nombre d’observations et la puissance de votre machine.
-
library(lcmm)
-mod4 <-lcmm(
- num_status ~ month, random =~ month, subject ='id',
-mixture =~ month, ng =4, idiag =TRUE, data = care_trajectories,
-link ="thresholds"
-)
+
library(lcmm)
+mod4 <-lcmm(
+ num_status ~ month, random =~ month, subject ='id',
+mixture =~ month, ng =4, idiag =TRUE, data = care_trajectories,
+link ="thresholds"
+)
Voyons comment se présentent les résultats.
-
summary(mod4)
+
summary(mod4)
General latent class mixed model
fitted by maximum likelihood method
@@ -31066,7 +33762,7 @@
Modèle mixte à classes latentes
thresh. parm3 1.08980 0.00727 149.816 0.00000
On dispose d’un AIC et d’un BIC. Ainsi, une stratégie possible pour déterminer le nombre de classes consiste à calculer un modèle différent pour chaque nombre de classes envisagé puis à retenir le modèle ayant le plus faible AIC ou BIC.
Pour chaque observation, le modèle a calculé la probabilité qu’elle appartienne à chacune des 4 classes identifiées. La fonction postprob fournit des statistiques sur cette classification.
Une manière alternative de présenter les classes consiste à représenter chaque mois, non pas la distribution dans chaque état, mais un état moyen en considérant que le statut dans les soins peut être assimilé à un score allant de 1 à 4.
Il faut cependant rester vigilant, dans la mesure où ce type de représentation synthétique peut masquer la diversité des trajectoires sous-jacentes, notamment en termes de durée ou d’enchaînement des événements. Même si cela est peut-être plus difficile à lire, il est toujours bon de regarder les tapis de séquences.
-
care_trajectories[, tmp_rang :=as.integer(fct_infreq(factor(id))), by = mod4_class]
-
ggplot(care_trajectories[month <=36]) +
-aes(x = month, y = tmp_rang, fill = care_statusF) +
-geom_raster() +
-facet_grid(mod4_class2 ~ ., space ="free", scales ="free") +
-scale_x_continuous(breaks =0:6*6, labels =paste0("M", 0:6*6)) +
-scale_y_continuous(breaks =0:5*200, minor_breaks =NULL) +
-xlab("") +ylab("") +
-theme_light() +
-theme(legend.position ="bottom") +
-labs(fill ="Statut dans les soins") +
-scale_fill_viridis(discrete =TRUE, direction =-1) +
-guides(fill =guide_legend(nrow =2))
ggplot(care_trajectories %>%filter(month <=36)) +
+aes(x = month, y = tmp_rang, fill = care_statusF) +
+geom_raster() +
+facet_grid(mod4_class2 ~ ., space ="free", scales ="free") +
+scale_x_continuous(breaks =0:6*6, labels =paste0("M", 0:6*6)) +
+scale_y_continuous(breaks =0:5*200, minor_breaks =NULL) +
+xlab("") +ylab("") +
+theme_light() +
+theme(legend.position ="bottom") +
+labs(fill ="Statut dans les soins") +
+scale_fill_viridis(discrete =TRUE, direction =-1) +
+guides(fill =guide_legend(nrow =2))
@@ -31242,24 +34037,28 @@
Modèle à observations répétées
Il s’agit de modèles classiques sauf qu’au lieu de considérer une ligne par individu, nous allons intégrer dans le modèle une ligne par individu et par pas de temps. Dans la mesure où nous avons plusieurs observations pour une même personne, cela doit être pris en compte par l’ajout d’un effet aléatoire dans le cadre d’un modèle mixte ou en ayant recours à un (voir le chapitre sur les modèles à effets aléatoires).
Vue la nature de notre variable d’intérêt (plusieurs modalités ordonnées), nous aurons recours à une régression logistique ordinale (voir le chapitre dédié). Pour un modèle mixte ordinal on peut utiliser la fonction clmm de l’extension ordinal. Pour un modèle GEE ordinal, citons ordgee de l’extension geepack ou encore ordLORgee de multgee. Il importe également que la dimension temporelle soit inclue dans les variables du modèle.
Ici, nous allons utiliser ordgee. Il nous faut tout d’abord transformer notre variable d’intérêt en un facteur ordonné.
Nous allons transformer nos variables explicatives en facteurs. Pour le temps, dans la mesure où sont effet n’est pas forcément linéaire, nous allons l’intégrer en tant que variable catégorielle. Par contre, comme nous n’avons que très peu d’observations individuelles après 3 ans, nous ne prendrons en compte que les observations des 36 premiers mois. Nous allons aussi retirer les observations à M0 puisqu’à ce moment précis tous les individus sont dans la même situation (diagnostiqués mais pas en soins.)
Les coefficients du modèle s’obtiennent avec summary. Malheureusement, il n’existe pas de tieder pour ce type de modèle. Nous allons donc procéder manuellement.
-
res <-summary(mod_td)$mean
-res$term <-rownames(res)
-head(res)
+
res <-summary(mod_td)$mean
+res$term <-rownames(res)
+head(res)
@@ -31321,26 +34120,26 @@
Modèle à observations répétées
Les intervalles de confiance à 95% ne sont pas déjà calculés. Faisons-le donc nous même.
-
mult <- stats::qnorm((1+ .95) /2)
-res$conf.low <- res$estimate - mult * res$san.se
-res$conf.high <- res$estimate + mult * res$san.se
+
mult <- stats::qnorm((1+ .95) /2)
+res$conf.low <- res$estimate - mult * res$san.se
+res$conf.high <- res$estimate + mult * res$san.se
Enfin, nous souhaitons disposer des odds ratios et non des coefficients bruts. Il faut avoir recours à la fonction exp (exponentielle).
Préparons un tableau avec les résultats. Pour le rendre plus lisible, nous allons mettre en forme les odds ratios avec une seule décimale. Améliorer le rendu des p-values et nous allons utiliser la virgule comme séparateur de décimal, comme il se doit. Nous aurons recours aux fonctions number et pvalue de l’extension scales (voir le chapitre sur la mise en forme des nombres).
On peut facilement représenter tout cela graphiquement. On va supprimer les termes seuils, grâce à str_detect de stringr. On notera le recours à fct_inorder de forcats pour conserver l’ordre des termes selon leur ordre d’apparition.
Sur ce graphique, on visualise bien l’évolution temporelle traduite par les odds ratios associés à chaque mois, ainsi que les effets globaux de nos covariables : les femmes ont une meilleure progression dans la cascade de soins que les hommes, de même que les plus éduqués et les plus âgés.
@@ -31712,13 +34506,15 @@
Modèle à observations répétées
Modèle de survie multi-états
Depuis la fin du XXe siècle, de nombreux développements ont réalisés pour étendre les modèles de survie à des processus multi-états. Ces modèles permettent de considérer une grande variété de processus. Plusieurs implémentations existent dans R3. Ici, nous allons utiliser l’extension msm qui repose sur des modèles de Markov multi-états et peut prendre en compte des co-variables dans le modèle.
En premier lieu, pour cette extension, ils nous faut disposer des données sous une forme longue, c’est-à-dire avec une ligne par individu et point d’observation dans le temps, ce qui est déjà le cas du fichier care_trajectories. Les différents status possibles doivent également être codés sous la forme de nombres entiers croissants (ici 1 correspondra à D, 2 à C, 3 à T et 4 à S).
Par ailleurs, nous n’allons concerver dans l’analyse que les individus avec au moins deux points d’observation, ici ceux observés au moins jusqu’à un mois.
Par ailleurs, nous n’allons conserver dans l’analyse que les individus avec au moins deux points d’observation, ici ceux observés au moins jusqu’à un mois.
La fonction statetable.msm permet de calculer le nombre et le type de transitions observées dans les données.
-
statetable.msm(status, id, data = ct)
+
statetable.msm(status, id, data = ct)
from/to
@@ -31759,49 +34555,49 @@
Modèle de survie multi-états
Il faut ensuite définir les transitions possibles dans le modèle en faisant une matrice carrée. On indiquera 0 si la transition n’est pas possible, une valeur positive sinon.
-
tr <-rbind(
-c(0, 1, 0, 0), # de 1 : vers 2
-c(1, 0, 1, 0), # de 2 : vers 1 ou vers 3
-c(0, 1, 0, 1), # de 3 : vers 2 ou vers 4
-c(0, 0, 1, 0) # de 4 : vers 3
-)
+
tr <-rbind(
+c(0, 1, 0, 0), # de 1 : vers 2
+c(1, 0, 1, 0), # de 2 : vers 1 ou vers 3
+c(0, 1, 0, 1), # de 3 : vers 2 ou vers 4
+c(0, 0, 1, 0) # de 4 : vers 3
+)
Dans notre matrice de transitions, nous n’avons pas défini de transition directe entre le statut 1 et le statut 3, alors que de telles transitions sont pourtant observées dans notre fichier. En fait, nous pouvons considérer qu’une transition de 1 vers 3 correspond en fait à deux transitions successives, de 1 vers 2 puis de 2 vers 3. La fonction msm s’aura identifier d’elle-mêmes ces doubles, voire triples, transitions.
On peut facilement représenter notre matrice de transition sous forme de schéma à l’aide de l’excellente extension DiagrammeR.
-
library(DiagrammeR)
-mermaid("
-graph TD
-1[diagnostiqué, mais pas suivi]
-2[suivi, mais pas sous traitement]
-3[sous traitement, mais infection non contrôlée]
-4[sous traitement et infection contrôlée]
-
-1--entrée<br />en soins-->2
-2--initiation<br />traitement-->3
-2--sortie<br />de soins-->1
-3--contrôle<br />infection-->4
-3--arrêt<br />traitement-->2
-4--échec<br/>virologique-->3
-", height =300)
+
library(DiagrammeR)
+mermaid("
+graph TD
+1[diagnostiqué, mais pas suivi]
+2[suivi, mais pas sous traitement]
+3[sous traitement, mais infection non contrôlée]
+4[sous traitement et infection contrôlée]
+
+1--entrée<br />en soins-->2
+2--initiation<br />traitement-->3
+2--sortie<br />de soins-->1
+3--contrôle<br />infection-->4
+3--arrêt<br />traitement-->2
+4--échec<br/>virologique-->3
+", height =300)
Il ne nous reste plus qu’à spécifier notre modèle. L’option obstype = 1 indique à msm que nos données correspondent à des snapshots à certains moments donnés (ici tous les mois) et donc que les transitions d’un état à un autre ont eu lieu entre nos points d’observation. Les types 2 et 3 correspondent à des dates de transition exactes (voir l’aide la fonction pour plus de détails).
-
ms_mod <-msm(
- status ~ month, subject = id, data = ct, qmatrix = tr, obstype =1
-)
+
ms_mod <-msm(
+ status ~ month, subject = id, data = ct, qmatrix = tr, obstype =1
+)
En exécutant cette commande, vous risquez d’obtenir le message d’erreur suivant :
Error in Ccall.msm(params, do.what = "lik", ...) : numerical overflow in calculating likelihood
Cela est dû à un problème d’échelle dans l’optimisation du modèle qui génère des nombres plus grands que ce que peux gérer l’ordinateur. Il peut être résolu de la manière suivante. Tout d’abord, on reexécute le modèle avec l’option control = list(trace = TRUE).
initial value 74796.800445
Error in Ccall.msm(params, do.what = "lik", ...) : numerical overflow in calculating likelihood
Ce qui importe de retenir, c’est la valeur initiale du paramètre d’optimisation avant que l’erreur ne se produise. On va l’utiliser (ou une valeur proche) comme paramètre d’échelle pour l’optimation avec l’option fnscale :
initial value 0.997291
iter 10 value 0.520302
iter 20 value 0.497598
@@ -31810,19 +34606,19 @@
Modèle de survie multi-états
converged
Used 39 function and 37 gradient evaluations
On peut comparer la prévalence dans chaque état au cours du temps telle que modélisée par le modèle avec les valeurs observées avec la fonction plot.prevalence.msm.
-
plot.prevalence.msm(ms_mod)
+
plot.prevalence.msm(ms_mod)
Par défaut, msm considère que les intensités de transition d’un état à un autre sont constantes au cours du temps. Or, dans notre example, il apparait que les prévalences observées varient différemment pendant les premiers mois après le diagnostic. Nous allons donc recalculer le modèle en spécifiant avec le paramètre pci que nous souhaitons considérer des intensités de transition différentes pour les trois premiers mois, la première année, la seconde année et après la seconde année. Comme il faudra plus d’itérations pour faire converger notre modèle, nous avons également augmenter la valeur du paramètre maxit (100 par défaut).
Ceci dit, le format dans lequel sont renvoyées les prévalences n’est que peu pratique pour les exploiter ensuite, par exemple avec ggplot2. L’extension JLutils fournit une fonction expérimentale tidy.prevalence.msm4 permettant de transformer ce résultat dans un format tidy.
-
library(JLutils)
-prev <-tidy.prevalence.msm(prevalence.msm(ms_mod, times =0:36))
-head(prev)
+
Ceci dit, le format dans lequel sont renvoyées les prévalences n’est que peu pratique pour les exploiter ensuite, par exemple avec ggplot2. L’extension JLutils fournit une fonction expérimentale tidy.prevalence.msm4 permettant de transformer ce résultat dans un format tidy. JLutils est seulement disponible sur GitHub. On l’installera donc (ou on la mettra à jour) avec la commande devtools::install_github("larmarange/JLutils").
+
library(JLutils)
+prev <-tidy.prevalence.msm(prevalence.msm(ms_mod, times =0:36))
+head(prev)
time
@@ -31942,43 +34738,43 @@
Modèle de survie multi-états
-
prev$status <-to_factor(prev$status)
-casc_status <-c(
-"diagnostiqué, mais pas suivi",
-"suivi, mais pas sous traitement",
-"sous traitement, mais infection non contrôlée",
-"sous traitement et infection contrôlée"
-)
-levels(prev$status) <- casc_status
+
prev$status <-to_factor(prev$status)
+casc_status <-c(
+"diagnostiqué, mais pas suivi",
+"suivi, mais pas sous traitement",
+"sous traitement, mais infection non contrôlée",
+"sous traitement et infection contrôlée"
+)
+levels(prev$status) <- casc_status
Il est alors ensuite facile de produire le graphique de la cascade de soins, estimée par le modèle, que l’on pourra mettre en comparaison de la cascade observée que nous avions calculé tout à l’heure.
-
casc_est <-ggplot(prev) +
-aes(x = time, fill = status, weight = expected) +
-geom_bar(color ="gray50", width =1, position ="fill") +
-scale_x_continuous(breaks =0:6*6, labels =paste0("M", 0:6*6)) +
-scale_y_continuous(labels = scales::percent) +
-ggtitle("Cascade des soins estimée, selon le temps depuis le diagnostic") +
-xlab("") +ylab("") +
-theme_light() +
-theme(legend.position ="bottom") +
-labs(fill ="Statut dans les soins") +
-scale_fill_viridis(discrete =TRUE, direction =-1) +
-guides(fill =guide_legend(nrow =2))
-multiplot(casc_obs, casc_est)
+
casc_est <-ggplot(prev) +
+aes(x = time, fill = status, weight = expected) +
+geom_bar(color ="gray50", width =1, position ="fill") +
+scale_x_continuous(breaks =0:6*6, labels =paste0("M", 0:6*6)) +
+scale_y_continuous(labels = scales::percent) +
+ggtitle("Cascade des soins estimée, selon le temps depuis le diagnostic") +
+xlab("") +ylab("") +
+theme_light() +
+theme(legend.position ="bottom") +
+labs(fill ="Statut dans les soins") +
+scale_fill_viridis(discrete =TRUE, direction =-1) +
+guides(fill =guide_legend(nrow =2))
+cowplot::plot_grid(casc_obs, casc_est, ncol =1)
Comme on peut le voir, dans le modèle, l’évolution est plus lissée comparativement aux données brutes observées.
Un des intérêts de msm est la possibilité d’y intégrer des covariables, permettant ainsi de calculer un modèle multivarié. Les variables peuvent être dépendantes du temps, puisqu’il suffit de renseigner leur valeur à chaque point d’observation.
ms_mod_mult <-msm(
+ status ~ month, subject = id, data = ct, qmatrix = tr, obstype =1,
+pci =c(3, 12, 24),
+control =list(fnscale =75000, trace =TRUE, maxit =500),
+covariates =~ sex + age + education
+)
Les risques relatifs, ou hazard ratios en anglais, associés à chaque covariable et à chaque transition s’obtiennent avec hazard.msm. Une fois encore, le format de sortie n’est pas le plus adapté pour un traitement graphique, mais on pourra avoir recours à tidy.hazard.msm de JLutils.
On va recoder certaines étiquettes en vue de faire un graphique des résultats.
-
hr$type <-"Transition ascendante"
-hr[hr$transition ==c("State 2 - State 1"),]$type <-"Transition descendante"
-hr[hr$transition ==c("State 3 - State 2"),]$type <-"Transition descendante"
-hr[hr$transition ==c("State 4 - State 3"),]$type <-"Transition descendante"
-
-hr$term <-fct_recode(
-fct_inorder(hr$term),
-"femme vs. homme"="sexfemme",
-"30-59 vs. 16-29"="age30-59",
-"60+ vs. 16-29"="age60+",
-"éduc. secondaire vs. primaire"="educationsecondaire",
-"éduc. supérieure vs. primaire"="educationsupérieur"
-)
-
-hr$transition <-fct_recode(
-fct_inorder(hr$transition),
-"entrée en soins"="State 1 - State 2",
-"sortie de soins"="State 2 - State 1",
-"initiation traitement"="State 2 - State 3",
-"arrêt traitement"="State 3 - State 2",
-"contrôle infection"="State 3 - State 4",
-"échec virologique"="State 4 - State 3"
-)
-
-head(hr)
+
hr$type <-case_when(
+ hr$transition ==c("State 2 - State 1") ~"Transition descendante",
+ hr$transition ==c("State 3 - State 2") ~"Transition descendante",
+ hr$transition ==c("State 4 - State 3") ~"Transition descendante",
+TRUE~"Transition ascendante"
+)
+
+hr$term <- hr$term %>%
+fct_inorder() %>%
+fct_recode(
+"femme vs. homme"="sexfemme",
+"30-59 vs. 16-29"="age30-59",
+"60+ vs. 16-29"="age60+",
+"éduc. secondaire vs. primaire"="educationsecondaire",
+"éduc. supérieure vs. primaire"="educationsupérieur"
+ )
+
+hr$transition <- hr$transition %>%
+fct_inorder() %>%
+fct_recode(
+"entrée en soins"="State 1 - State 2",
+"sortie de soins"="State 2 - State 1",
+"initiation traitement"="State 2 - State 3",
+"arrêt traitement"="State 3 - State 2",
+"contrôle infection"="State 3 - State 4",
+"échec virologique"="State 4 - State 3"
+ )
+
+head(hr)
-
-
term
-
transition
-
from
-
to
-
estimate
-
conf.low
-
conf.high
-
type
+
+
+term
+
+
+transition
+
+
+from
+
+
+to
+
+
+estimate
+
+
+conf.low
+
+
+conf.high
+
+
+type
+
-
-
femme vs. homme
-
entrée en soins
-
State 1
-
State 2
-
1.8855232
-
1.7059612
-
2.083985
-
Transition ascendante
+
+
+femme vs. homme
+
+
+entrée en soins
+
+
+State 1
+
+
+State 2
+
+
+1.8855232
+
+
+1.7059612
+
+
+2.083985
+
+
+Transition ascendante
+
-
-
femme vs. homme
-
sortie de soins
-
State 2
-
State 1
-
0.9245660
-
0.7655997
-
1.116540
-
Transition descendante
+
+
+femme vs. homme
+
+
+sortie de soins
+
+
+State 2
+
+
+State 1
+
+
+0.9245660
+
+
+0.7655997
+
+
+1.116540
+
+
+Transition descendante
+
-
-
femme vs. homme
-
initiation traitement
-
State 2
-
State 3
-
0.9015913
-
0.8044458
-
1.010468
-
Transition ascendante
+
+
+femme vs. homme
+
+
+initiation traitement
+
+
+State 2
+
+
+State 3
+
+
+0.9015913
+
+
+0.8044458
+
+
+1.010468
+
+
+Transition ascendante
+
-
-
femme vs. homme
-
arrêt traitement
-
State 3
-
State 2
-
1.0959810
-
0.8854966
-
1.356498
-
Transition descendante
+
+
+femme vs. homme
+
+
+arrêt traitement
+
+
+State 3
+
+
+State 2
+
+
+1.0959810
+
+
+0.8854966
+
+
+1.356498
+
+
+Transition descendante
+
-
-
femme vs. homme
-
contrôle infection
-
State 3
-
State 4
-
1.3985443
-
1.2300980
-
1.590057
-
Transition ascendante
+
+
+femme vs. homme
+
+
+contrôle infection
+
+
+State 3
+
+
+State 4
+
+
+1.3985443
+
+
+1.2300980
+
+
+1.590057
+
+
+Transition ascendante
+
-
-
femme vs. homme
-
échec virologique
-
State 4
-
State 3
-
0.8610144
-
0.6678918
-
1.109979
-
Transition descendante
+
+
+femme vs. homme
+
+
+échec virologique
+
+
+State 4
+
+
+State 3
+
+
+0.8610144
+
+
+0.6678918
+
+
+1.109979
+
+
+Transition descendante
+
Vu le nombre de coefficients (un risque relatif par covariable et par transition), on va organiser le graphique en distinguant les transitions acsendantes et les transitions descendantes, d’une part, et en regroupant les risques relatifs d’une même covariable, d’autre part. Pour alléger le graphique, nous allons également retirer la covariable timeperiod créé par l’argument pci, en ayant recours à str_detect de l’extension stringr pour repérer les lignes en questions (voir le chapitre sur la manipulation de texte).
Ce modèle de survie multi-états permet de mettre en évidence des effets différenciés selon la transition considérée. Par exemple, les femmes sont plus rapides en matière d’entrée en soins et de contrôle de l’infection (une fois le traitement initié) mais plus lentes à démarrer le traitement (une fois entrées en soins). Par contre, le sexe ne semble pas jouer sur les transitions descendantes (échec virologique, arrête du traitement ou sortie de soins).
Ce chapitre est évoqué dans le webin-R #13 (exemples de graphiques avancés) sur YouTube.
+
Ce chapitre est évoqué dans le webin-R #14 (exemples de graphiques avancés - 2) sur YouTube.
+
De nombreuses extensions permettent d’étendre les possibilités graphiques de ggplot2. Certaines ont déjà été abordées dans les différents chapitres d’analyse-R. Le présent chapitre ne se veut pas exhaustif et ne présente qu’une sélection choisie d’extensions.
L’extension ggh4x propose une géométrie geom_pointpath pour afficher des points reliés entre eux par un tiret. La géométrie geom_pointpath est également proposée par l’extension lemon.
+
library(ggh4x)
+
+Attaching package: 'ggh4x'
+
The following objects are masked from 'package:lemon':
+
+ geom_pointpath, GeomPointPath
+
ggplot(pressure) +
+aes(x = temperature, y = pressure) +
+geom_pointpath()
+
+
Accolades de comparaison (bracket)
La géométrie geom_braket de l’extension ggpubr permets d’ajouter sur un graphique des accolades de comparaison entre groupes.
L’extension ggridges fournit une géométrie geom_density_ridges_gradient pour la création de diagramme en crêtes.
-
library(ggridges)
-ggplot(lincoln_weather, aes(x =`Mean Temperature [F]`, y = Month, fill =stat(x))) +
-geom_density_ridges_gradient(scale =3, rel_min_height =0.01) +
-scale_fill_viridis_c(name ="Temp. [F]", option ="C") +
-labs(title ="Temperatures in Lincoln NE in 2016") +
-theme_ridges()
+
library(ggridges)
+ggplot(lincoln_weather, aes(x =`Mean Temperature [F]`, y = Month, fill =stat(x))) +
+geom_density_ridges_gradient(scale =3, rel_min_height =0.01) +
+scale_fill_viridis_c(name ="Temp. [F]", option ="C") +
+labs(title ="Temperatures in Lincoln NE in 2016") +
+theme_ridges()
Par défaut, des encoches (ticks) sont dessinées sur les axes uniquement pour la grille principale (major breaks). La fonction guide_axis_minor de l’extension ggh4x permet de rajouter des encoches aux points de la grille mineure (minor breaks).
+
p <-ggplot(mtcars) +
+aes(x = cyl, y = mpg) +
+geom_point() +
+theme_classic()
+p
+
+
library(ggh4x)
+p +
+scale_x_continuous(
+minor_breaks =16:32/4,
+guide =guide_axis_minor()
+ ) +
+# to control relative length of minor ticks
+theme(ggh4x.axis.ticks.length.minor =rel(.8))
+
+
+
+
Relations imbriquées
+
La fonction guide_axis_nested de l’extension ggh4x permet d’afficher sur un axe des relations imbriquées.
L’extension ggcharts fournit plusieurs fonctions de haut niveau pour faciliter la réalisation de graphiques divergents en barres (diverging_bar_chart), en sucettes (diverging_lollipop_chart) voire même une pyramide des âges (pyramid_chart).
library(ggplot2)
+library(gganimate)
+library(gapminder)
+
+ggplot(gapminder, aes(gdpPercap, lifeExp, size = pop, colour = country)) +
+geom_point(alpha =0.7, show.legend =FALSE) +
+scale_colour_manual(values = country_colors) +
+scale_size(range =c(2, 12)) +
+scale_x_log10() +
+facet_wrap(~continent) +
+# Here comes the gganimate specific bits
+labs(title ="Year: {frame_time}", x ="GDP per capita", y ="life expectancy") +
+transition_time(year) +
+ease_aes("linear")
+
Voir le site de l’extension (https://gganimate.com/) pour la documentation et des tutoriels. Il est conseillé d’installer également l’extension gifski avec gganimate.
Surligner certaines données
L’extension gghighlight fournit une fonction gghiglight qui permets de surligner les données qui remplissent des conditions spécifiées.
L’extension hrbrthemes fournit plusieurs thèmes graphiques pour ggplot2. Un exemple ci-dessous. Pour plus d’informations, voir https://github.com/hrbrmstr/hrbrthemes.
-
library(ggplot2)
-library(hrbrthemes)
-ggplot(mtcars, aes(mpg, wt)) +
-geom_point(aes(color =factor(carb))) +
-labs(
-x ="Fuel efficiency (mpg)", y ="Weight (tons)",
-title ="Seminal ggplot2 scatterplot example",
-subtitle ="A plot that is only useful for demonstration purposes",
-caption ="Brought to you by the letter 'g'"
- ) +
-scale_color_ipsum() +
-theme_ipsum_rc()
-
+
library(ggplot2)
+library(hrbrthemes)
+ggplot(mtcars, aes(mpg, wt)) +
+geom_point(aes(color =factor(carb))) +
+labs(
+x ="Fuel efficiency (mpg)", y ="Weight (tons)",
+title ="Seminal ggplot2 scatterplot example",
+subtitle ="A plot that is only useful for demonstration purposes",
+caption ="Brought to you by the letter 'g'"
+ ) +
+scale_color_ipsum() +
+theme_ipsum_rc()
Ce chapitre est évoqué dans le webin-R #13 (exemples de graphiques avancés) sur YouTube.
+
Ce chapitre est évoqué dans le webin-R #14 (exemples de graphiques avancés - 2) sur YouTube.
Dans ce chapitre, nous présentons plusieurs graphiques avec une mise en forme avancée et détaillons pas à pas leur création.
Chargeons quelques extensions de base.
@@ -33555,9 +36362,302 @@
Code final du graphique
+
+
Évolutions temporelles
+
Pour cet exemple, nous allons utiliser des données fictives. Dans les années 1950, au sortir de la guerre, la Syldavie et la Bordurie (deux pays fictifs inventés par Hergé) ont décidé de mettre un programme de distribution de bandes dessinées de Tintin dans des écoles et des collèges afin de motiver les élèves à lire plus. Nous disposons des rapports d’activités des différents agents du programme qui indident, par mois et par type d’établissement, le nombre d’élèves rencontrés et le nombre de bandes dessinées distribués.
+
Commençons par importer les données, qui sont au format Excel. Nous utiliserons donc la fonction read_excel de l’extension readxl. Profitons-en pour charger également en mémoire le tidyverse ainsi que l’extension lubridate qui nous servira pour la gestion des dates.
En premier lieu, nous souhaiterions représenter le nombre de bandes dessinées distribuées par mois. Comme on peut le voir, nous n’avons pas directement une variable date, l’année et le mois étant indiqués dans deux variables différentes. Nous allons donc commencer par créer une telle variable, qui correspondra au premier jour du mois.
+
Nous allons utiliser la fonction paste0 qui permets de concaténer du texte pour créer des chaines de texte au format ANNEE-MOIS-JOUR puis transformer cela en objet Date avec ymd de l’extension lubridate.
Min. 1st Qu. Median Mean
+"1949-12-01" "1950-07-01" "1950-12-01" "1950-11-13"
+ 3rd Qu. Max.
+"1951-04-01" "1951-07-01"
+
Nous pouvons maintenant commencer à réaliser notre graphique. Nous allons compter le nombre de BD distribuées par mois et représenter cela sous forme d’un graphique en barres. geom_bar utilise la statistique stat_count pour compter le nombre d’observations. Cependant, le nombre de BD distribuées varie d’une ligne du fichier à l’autre. Nous allons donc utiliser l’esthétique weight pour pondérer le calcul et obtenir le total de bandes dessinées distribuées. Nous allons également utiliser l’esthétique fill pour distinguer les écoles primaires et les collèges.
+
p <-ggplot(bd) +
+aes(x = date, weight = bd_distribuees, fill = type_etablissement) +
+geom_bar()
+p
+
+
Avec facet_grid nous pouvons comparer la Syldavie et la Bordurie.
+
p +facet_grid(cols =vars(pays))
+
+
Essayons de faire un peu mieux en répétant l’axe des y sur le graphique de droite, ce qui peut être fait avec la fonction facet_rep_grid de l’extension lemon.
De même, nous allons utiliser une palette de couleurs adaptées aux personnes daltoniennes, en s’appuyant sur scale_fill_bright disponible avec khroma.
+
library(khroma)
+p <- p +scale_fill_bright()
+p
+
+
Reste à améliorer la présentation de l’axe des x avec scale_x_date{data-pkg="ggplot2}. L’option date_labels permet de personnaliser l’affichage des dates via des raccourcis décrits dans l’aide de la fonction strptime. Ainsi, %b correspond au mois écrit textuellement de manière abrégée et %Y à l’année sur 4 chiffres. \n est utilisé pour indiquer un retour à la ligne. L’option date_breaks permet d’indiquer de manière facilitée la durée attendue entre deux étiquettes, ici trois mois.
+
p <- p +scale_x_date(
+date_labels ="%b\n%Y",
+date_breaks ="3 months"
+)
+p
+
+
Pour améliorer le graphique, nous aimerions rajouter une encoche sous chaque mois. Ceci est rendu possible avec la fonction guide_axis_minor de ggh4x.
Évolution du nombre moyen de bandes dessinées remises par élève
+
Nous disposons dans les rapports d’activités du nombre de bandes dessinées distribuées et du nombre d’élèves rencontrés. Nous pouvons donc regarder comment à évaluer au cours du temps le nombre moyen de bandes dessinées distribuées par élève, à savoir l’évolution du ratio bd_distribuees / eleves_rencontres.
Comme on le voit sur la figure précédente, si l’on représente directement cet indicateur, nous obtenons plusieurs points pour chaque date. En effet, nous avons plusieurs observations par date et nous obtenons donc un point pour chacune de ces observations. Nous souhaiterions avoir la moyenne pour chaque date. Cela est possible en ayant recours à stat_weighted_mean de l’extension GGally. Pour que le calcul soit correct, il faut réaliser une moyenne pondérée par le dénominateur, à savoir ici le nombre d’élèves.
Pour améliorer la comparaison, il serait pertinent d’ajouter les intervalles de confiance sur le graphique. Cela peut être réalisé par exemple avec des rubans de couleur. Mais pour cela, il va falloir avoir recours à quelques astuces pour permettre le calcul de ces intervalles de confiance à la volée.
+
Tout d’abord, pour calculer l’intervalle de confiance d’un ratio de nombres entiers, nous pouvons utiliser un test de Poisson avec poisson.test auquel nous devons transmettre le numérateur et le dénominateur. Par exemple, pour 48 bandes dessinées distribuées à 27 élèves :
+
test <-poisson.test(48, 27)
+test
+
+ Exact Poisson test
+
+data: 48 time base: 27
+number of events = 48, time base = 27, p-value =
+0.0002246
+alternative hypothesis: true event rate is not equal to 1
+95 percent confidence interval:
+ 1.310793 2.357075
+sample estimates:
+event rate
+ 1.777778
+
str(test)
+
List of 9
+ $ statistic : Named num 48
+ ..- attr(*, "names")= chr "number of events"
+ $ parameter : Named num 27
+ ..- attr(*, "names")= chr "time base"
+ $ p.value : num 0.000225
+ $ conf.int : num [1:2] 1.31 2.36
+ ..- attr(*, "conf.level")= num 0.95
+ $ estimate : Named num 1.78
+ ..- attr(*, "names")= chr "event rate"
+ $ null.value : Named num 1
+ ..- attr(*, "names")= chr "event rate"
+ $ alternative: chr "two.sided"
+ $ method : chr "Exact Poisson test"
+ $ data.name : chr "48 time base: 27"
+ - attr(*, "class")= chr "htest"
+
Comme nous le montre str, poisson.test renvoie une liste avec différents éléments et l’intervalle de confiance est disponible dans le sous-objet "conf.int". Il est possible d’y accéder avec l’opérateur $ ou bien avec la fonction pluck de l’extension purrr chargée par défaut avec le tidyverse.
Pour obtenir la borne inférieure de l’intervalle de confiance, il nous faut maintenant extraire le premier élément.
+
test$conf.int[1]
+
[1] 1.310793
+
test %>%pluck("conf.int", 1)
+
[1] 1.310793
+
Nous avons trouvé la manière de calculer la borne inférieure de l’intervalle de confiance pour une ligne d’observations. Mais comment procéder pour plusieurs observations. Supposons que nous ayons 4 observations :
+
num <-c(50, 25, 32, 48)
+denom <-c(25, 17, 25, 31)
+
Si nous passons ces vecteurs de données à poisson.test, nous obtenons un message d’erreur car poisson.test ne permets de calculer qu’un seul test à la fois.
+
poisson.test(num, denom)
+
Error in poisson.test(num, denom) : le cas k > 2 n'est pas implémenté
+
Dans notre cas, nous souhaitons exécuter poisson.test ligne à ligne. Cela peut être réalisée avec la famille de fonctions map fournies par l’extension purrr. Plus précisément, dans notre situation, nous allons avoir recours à map2 qui permet de fournir deux vecteurs en entrée et d’appliquer une fonction ligne à ligne. On passera à map2 soit une fonction à utiliser telle quelle ou bien une formule décrivant une fonction personnalisée à la levée, ce que nous allons faire. Dans cette formule décrivant le calcul à effectuer, on utilisera .x pour se référer à l’élément du premier vecteur et à .y pour se référer au deuxième vecteur.
Le résultat obtenu est une liste de 4 éléments, chaque élément contenu la borne inférieur des intervalles. Nous préférerions que la fonction nous retourne un vecteur de valeurs numériques (double) plutôt qu’une liste. Nous allons donc utiliser map2_dbl plutôt que map2.
Nous y sommes presque. Il nous reste plus qu’à avoir recours à after_stat pour avoir accès aux numérateurs et dénominateurs calculés par stat_weighted_mean.
Faisons quelques ajustements pour une meilleur lisibilité : effaçons les limites supérieures et inférieures des rubans avec color = "transparent", augmentons la transparence des rubans avec alpha = .2 et épaississons les courbes principales avec size = 1.
Pour des besoins plus pointus, on pourra aussi utiliser l’extension stringi sur laquelle est elle-même basée stringr.
-
+
@@ -38608,7 +41708,7 @@
Ressources
-
+
@@ -38620,10 +41720,13 @@
Mettre en forme des nombres
Ce chapitre est en cours d’écriture.
+
+
Ce chapitre est évoqué dans le webin-R #09 (Graphiques uni- et bivariés avec ggplot2) sur YouTube.
+
Le plus simple est d’avoir recours à l’extension scales et aux fonctions number, comma, percent, unit, ordinal et pvalue.
-
+
@@ -38651,6 +41754,10 @@
Couleurs et Palettes
+
+
Ce chapitre est évoqué dans le webin-R #08 (ggplot2 et la grammaires des graphiques) sur YouTube.
+
Ce chapitre est évoqué dans le webin-R #13 (exemples de graphiques avancés) sur YouTube.
+
Noms de couleur
Lorsque l’on doit indiquer à R une couleur, notamment dans les fonctions graphiques, on peut mentionner certaines couleurs en toutes lettres (en anglais) comme "red" ou "blue". La liste des couleurs reconnues par R est disponible sur http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf.
Ce chapitre est évoqué dans le webin-R #15 (analyses de survie) sur YouTube.
+
Ce chapitre est évoqué dans le webin-R #15 (analyse de survie) sur YouTube.
Ressources en ligne
@@ -904,10 +905,10 @@
Kaplan-Meier
survdiff(formula = Surv(time, deces) ~ sexe, data = enfants)
N Observed Expected (O-E)^2/E (O-E)^2/V
-sexe=masculin 762 94 66.2 11.7 21.9
-sexe=féminin 822 48 75.8 10.2 21.9
+sexe=masculin 762 94 66.2 11.6 21.8
+sexe=féminin 822 48 75.8 10.2 21.8
- Chisq= 21.9 on 1 degrees of freedom, p= 3e-06
+ Chisq= 21.8 on 1 degrees of freedom, p= 3e-06
Une fois encore, on aura recours à ggsurvplot pour représenter les courbes de survie.
ggsurvplot(km_sexe, conf.int =TRUE, risk.table =TRUE, pval =TRUE, data = enfants)
@@ -932,163 +933,162 @@
Modèle de Cox
data = enfants)
coef
-sexeféminin -0.810964
-milieu 0.657039
-richessepauvre -0.074736
-richessemoyen 0.328896
-richesseriche 0.357901
-richessetrès riche 0.471011
-structureun adulte -0.146177
-structuredeux adultes de même sexe 0.602447
-structuretrois adultes ou plus avec lien de parenté 0.053835
-structureadultes sans lien de parenté -0.126713
-educ2primaire -0.032997
-educ2secondaire ou plus -0.201774
-gpage_mere_naissance19 ou moins -0.310449
-gpage_mere_naissance30 et plus -0.004781
-rang_apres_idealoui 1.414030
+sexeféminin -0.809568
+milieu 0.656241
+richessepauvre -0.082223
+richessemoyen 0.318645
+richesseriche 0.353483
+richessetrès riche 0.464590
+structureun adulte -0.150231
+structuredeux adultes de même sexe 0.604592
+structuretrois adultes ou plus avec lien de parenté 0.049430
+structureadultes sans lien de parenté -0.131369
+educ2primaire -0.030251
+educ2secondaire ou plus -0.203903
+gpage_mere_naissance19 ou moins -0.310248
+gpage_mere_naissance30 et plus -0.002586
+rang_apres_idealoui 1.355106
exp(coef)
-sexeféminin 0.444429
-milieu 1.929072
-richessepauvre 0.927989
-richessemoyen 1.389433
-richesseriche 1.430324
-richessetrès riche 1.601612
-structureun adulte 0.864005
-structuredeux adultes de même sexe 1.826583
-structuretrois adultes ou plus avec lien de parenté 1.055311
-structureadultes sans lien de parenté 0.880986
-educ2primaire 0.967542
-educ2secondaire ou plus 0.817279
-gpage_mere_naissance19 ou moins 0.733117
-gpage_mere_naissance30 et plus 0.995230
-rang_apres_idealoui 4.112495
+sexeféminin 0.445050
+milieu 1.927533
+richessepauvre 0.921066
+richessemoyen 1.375263
+richesseriche 1.424019
+richessetrès riche 1.591361
+structureun adulte 0.860509
+structuredeux adultes de même sexe 1.830506
+structuretrois adultes ou plus avec lien de parenté 1.050672
+structureadultes sans lien de parenté 0.876895
+educ2primaire 0.970202
+educ2secondaire ou plus 0.815541
+gpage_mere_naissance19 ou moins 0.733265
+gpage_mere_naissance30 et plus 0.997417
+rang_apres_idealoui 3.877173
se(coef)
-sexeféminin 0.177804
-milieu 0.269893
-richessepauvre 0.250820
-richessemoyen 0.248013
-richesseriche 0.298660
-richessetrès riche 0.428707
-structureun adulte 0.600331
-structuredeux adultes de même sexe 0.376554
-structuretrois adultes ou plus avec lien de parenté 0.196517
-structureadultes sans lien de parenté 0.305436
-educ2primaire 0.205792
-educ2secondaire ou plus 0.367133
-gpage_mere_naissance19 ou moins 0.268083
-gpage_mere_naissance30 et plus 0.191732
-rang_apres_idealoui 0.603549
+sexeféminin 0.177806
+milieu 0.269928
+richessepauvre 0.250417
+richessemoyen 0.247868
+richesseriche 0.298425
+richessetrès riche 0.428583
+structureun adulte 0.600329
+structuredeux adultes de même sexe 0.376445
+structuretrois adultes ou plus avec lien de parenté 0.196666
+structureadultes sans lien de parenté 0.305478
+educ2primaire 0.205751
+educ2secondaire ou plus 0.366889
+gpage_mere_naissance19 ou moins 0.268062
+gpage_mere_naissance30 et plus 0.191557
+rang_apres_idealoui 0.602401
z
-sexeféminin -4.561
-milieu 2.434
-richessepauvre -0.298
-richessemoyen 1.326
-richesseriche 1.198
-richessetrès riche 1.099
-structureun adulte -0.243
-structuredeux adultes de même sexe 1.600
-structuretrois adultes ou plus avec lien de parenté 0.274
-structureadultes sans lien de parenté -0.415
-educ2primaire -0.160
-educ2secondaire ou plus -0.550
-gpage_mere_naissance19 ou moins -1.158
-gpage_mere_naissance30 et plus -0.025
-rang_apres_idealoui 2.343
+sexeféminin -4.553
+milieu 2.431
+richessepauvre -0.328
+richessemoyen 1.286
+richesseriche 1.184
+richessetrès riche 1.084
+structureun adulte -0.250
+structuredeux adultes de même sexe 1.606
+structuretrois adultes ou plus avec lien de parenté 0.251
+structureadultes sans lien de parenté -0.430
+educ2primaire -0.147
+educ2secondaire ou plus -0.556
+gpage_mere_naissance19 ou moins -1.157
+gpage_mere_naissance30 et plus -0.014
+rang_apres_idealoui 2.250
p
-sexeféminin 5.09e-06
-milieu 0.0149
-richessepauvre 0.7657
-richessemoyen 0.1848
-richesseriche 0.2308
-richessetrès riche 0.2719
-structureun adulte 0.8076
-structuredeux adultes de même sexe 0.1096
-structuretrois adultes ou plus avec lien de parenté 0.7841
-structureadultes sans lien de parenté 0.6782
-educ2primaire 0.8726
-educ2secondaire ou plus 0.5826
-gpage_mere_naissance19 ou moins 0.2469
-gpage_mere_naissance30 et plus 0.9801
-rang_apres_idealoui 0.0191
-
-Likelihood ratio test=38.44 on 15 df, p=0.0007773
+sexeféminin 5.29e-06
+milieu 0.0151
+richessepauvre 0.7427
+richessemoyen 0.1986
+richesseriche 0.2362
+richessetrès riche 0.2784
+structureun adulte 0.8024
+structuredeux adultes de même sexe 0.1083
+structuretrois adultes ou plus avec lien de parenté 0.8016
+structureadultes sans lien de parenté 0.6672
+educ2primaire 0.8831
+educ2secondaire ou plus 0.5784
+gpage_mere_naissance19 ou moins 0.2471
+gpage_mere_naissance30 et plus 0.9892
+rang_apres_idealoui 0.0245
+
+Likelihood ratio test=38.16 on 15 df, p=0.0008546
n= 1584, number of events= 142
De nombreuses variables ne sont pas significatives. Voyons si nous pouvons, avec la fonction step, améliorer notre modèle par minimisation de l’AIC ou Akaike Information Criterion (voir la section Sélection de modèles du chapitre sur la Régression logistique).
On peut obtenir facilement les coefficients du modèle avec l’excellente fonction tidy de l’extension broom. Ne pas oublier de préciser exponentiate = TRUE. En effet, dans le cas d’un modèle de Cox, l’exponentiel des coefficients corresponds au ratio des risques instantannés ou hazard ratio (HR) en anglais.
-
library(broom, quietly =TRUE)
-
Warning: package 'broom' was built under R version 4.0.5
Pour représenter ces rapports de risque, on peut ici encore avoir recours à la fonction ggcoef_model de l’extension GGally.
-
library(GGally, quietly =TRUE)
+
library(GGally, quietly =TRUE)
Registered S3 method overwritten by 'GGally':
method from
+.gg ggplot2
-
ggcoef_model(mod2, exponentiate =TRUE)
+
ggcoef_model(mod2, exponentiate =TRUE)
Coefficients du modèle avec ggcoef_model
@@ -1096,7 +1096,7 @@
Modèle de Cox
L’extension survminer fournit également une fonction ggforest qui permet de représenter de manière plus esthétique et complète les coefficients d’un modèle de Cox.
-
ggforest(mod2)
+
ggforest(mod2)
Coefficients du modèle avec ggforest
@@ -1106,17 +1106,17 @@
Modèle de Cox
Vérification de la validité du modèle
Un modèle de Cox n’est valable que sous l’hypothèse de la proportionnalité des risques relatifs. Selon cette hypothèse les résidus de Schoenfeld ne dépendent pas du temps. Cette hypothèse peut être testée avec la fonction cox.zph.
Une valeur de p inférieure à 5 % indique que l’hypothèse n’est pas vérifiée. Il apparaît que p est supérieur à 5 % globalement et pour chaque variable prise individuellement. Notre modèle est donc valide.
Il est possible de représenter la distribution des résidus de Schoenfeld à l’aide de ggcoxzph de l’extension survminer, afin de voir si leur répartition change au cours du temps.
-
ggcoxzph(test)
+
ggcoxzph(test)
Warning in regularize.values(x, y, ties, missing(ties),
na.rm = na.rm): collapsing to unique 'x' values
Exemple 2 : une régression logistique avec deux variables catégorielles
odds.ratio(0.6122, 0.7114)
[1] 0.6404246
De la même manière, les différents coefficients associés à groupe_age correspondent à la différence entre chaque groupe d’âges et sa modalité de référence (ici 16-24 ans), quand les autres variables (ici le sexe) sont à leur référence (ici les hommes).
-
Pour prédire la probabilité de faire du sport pour un profil particulier, il faut prendre en compte toutes les termes qui s’appliquent et qui s’ajoutent à l’intercept. Par exemple, pour une femme de 50 ans il faut considérer l’intercept (0.9021), le coefficient sexeFemme (-0.4455) et le coefficient groupe_ageè45-64 (-1.6535). Sa probabilité de faire du sport est donc de 23%.
+
Pour prédire la probabilité de faire du sport pour un profil particulier, il faut prendre en compte toutes les termes qui s’appliquent et qui s’ajoutent à l’intercept. Par exemple, pour une femme de 50 ans il faut considérer l’intercept (0.9021), le coefficient sexeFemme (-0.4455) et le coefficient groupe_age45-64 (-1.6535). Sa probabilité de faire du sport est donc de 23%.
Un graphique similaire peut s’obtenir très facilement en ayant recours à la fonction 
