-
Notifications
You must be signed in to change notification settings - Fork 40
/
effets-d-interaction.Rmd
328 lines (245 loc) · 14.4 KB
/
effets-d-interaction.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
---
title: "Effets d'interaction dans un modèle"
---
```{r options_communes, include=FALSE, cache=FALSE}
source("options_communes.R")
```
<div class="guide-R">
Une version actualisée de ce chapitre est disponible sur **guide-R** : [Interactions](https://larmarange.github.io/guide-R/analyses/interactions.html)
</div>
<div class="webin-R">
Ce chapitre est évoqué dans le webin-R #07 (régression logistique partie 2) sur [YouTube](https://youtu.be/BUo9i7XTLYQ).
</div>
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 <dfn data-index="effet d'interaction">effets d'interaction</dfn><dfn data-index="interaction"></dfn>^[Pour une présentation plus statistique et mathématique des effets d'interaction, on pourra se référer au [cours de Jean-François Bickel disponible en ligne](http://commonweb.unifr.ch/artsdean/pub/gestens/f/as/files/4665/9547_131825.pdf).].
## Exemple d'interaction
Reprenons le modèle que nous avons utilisé dans le chapitre sur la [régression logistique](regression-logistique.html).
```{r}
library(questionr)
data(hdv2003)
d <- hdv2003
d$sexe <- relevel(d$sexe, "Femme")
d$grpage <- cut(d$age, c(16, 25, 45, 65, 99), right = FALSE, include.lowest = TRUE)
d$etud <- d$nivetud
levels(d$etud) <- c(
"Primaire", "Primaire", "Primaire",
"Secondaire", "Secondaire", "Technique/Professionnel",
"Technique/Professionnel", "Supérieur"
)
d$etud <- addNAstr(d$etud, "Manquant")
library(labelled)
var_label(d$sport) <- "Pratique du sport ?"
var_label(d$sexe) <- "Sexe"
var_label(d$grpage) <- "Groupe d'âges"
var_label(d$etud) <- "Niveau d'étude"
var_label(d$relig) <- "Pratique religieuse"
var_label(d$heures.tv) <- "Nombre d'heures passées devant la télévision par jour"
```
Nous avions alors exploré les facteurs associés au fait de pratiquer du sport.
```{r}
mod <- glm(sport ~ sexe + grpage + etud + heures.tv + relig, data = d, family = binomial())
library(gtsummary)
tbl_regression(mod, 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`{data-pkg="effects" data-rdoc="effects"} de l'extension `effects`{.pkg}.
<figure>
```{r}
library(effects)
plot(allEffects(mod))
```
<figcaption>Représentation graphique des effets du modèle</figcaption>
</figure>
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.
```{r}
mod2 <- glm(sport ~ sexe * grpage + etud + heures.tv + relig, data = d, family = binomial())
tbl_regression(mod2, exponentiate = TRUE)
```
Commençons par regarder les effets du modèle.
<figure>
```{r}
plot(allEffects(mod2))
```
<figcaption>Représentation graphique des effets du modèle avec interaction entre le sexe et le groupe d'âge</figcaption>
</figure>
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`{data-pkg="stats"}.
```{r}
anova(mod2, test = "Chisq")
```
Jetons maintenant un oeil aux coefficients du modèle. Pour rendre les choses plus visuelles, nous aurons recours à `ggcoef_model`{data-pkg="GGally"} de l'extension `GGally`{.pkg}.
<figure>
```{r}
library(GGally)
ggcoef_model(mod2, exponentiate = TRUE)
```
<figcaption>Représentation graphique des coefficients du modèle avec interaction entre le sexe et le groupe d'âge</figcaption>
</figure>
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
<figure>
```{r}
library(breakDown)
library(ggplot2)
logit <- function(x) exp(x)/(1+exp(x))
nouvelle_observation <- d[1, ]
nouvelle_observation$sexe[1] = "Femme"
nouvelle_observation$grpage[1] = "[45,65)"
nouvelle_observation$etud[1] = "Primaire"
nouvelle_observation$relig[1] = "Pratiquant regulier"
nouvelle_observation$heures.tv[1] = 0
plot(
broken(mod2, nouvelle_observation, predict.function = betas),
trans = logit
) + ylim(0, 1) + ylab("Probabilité de faire du sport")
```
<figcaption>Représentation graphique de l'estimation de la probabilité de faire du sport pour une femme de 60 ans</figcaption>
</figure>
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). <q>Femme</q> étant la modalité de référence pour la variable <var>sexe</var>, cela ne modifie pas le calcul de la probabilité de faire du sport. Par contre, il y a une modification induite par la modalité <q>45-65</q> de la variable <var>grpage</var>.
Regardons maintenant la situation d'un homme de 20 ans.
<figure>
```{r}
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")
```
<figcaption>Représentation graphique de l'estimation de la probabilité de faire du sport pour un homme de 20 ans</figcaption>
</figure>
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 <var>sexe</var> 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.
<figure>
```{r}
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")
```
<figcaption>Représentation graphique de l'estimation de la probabilité de faire du sport pour un homme de 60 ans</figcaption>
</figure>
Cette fois-ci, il y a plusieurs modifications d'effet. On applique en effet à la fois le coefficient <q>sexe = Homme</q> (effet du sexe pour les 15-24 ans), le coefficient <q>grpage = [45-65)</q> qui est l'effet de l'âge pour les femmes de 45-64 ans et le coefficient <q>sexe:grpage = Homme:[45-65)</q> 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`.
```{r}
mod3 <- glm(sport ~ sexe : grpage + etud + heures.tv + relig, data = d, family = binomial())
tbl_regression(mod3, exponentiate = TRUE)
```
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`{data-pkg="stats"}.
```{r}
anova(mod2, mod3, test = "Chisq")
```
De même, les effets modélisés sont les mêmes.
<figure>
```{r}
plot(allEffects(mod3))
```
<figcaption>Représentation graphique des effets du modèle avec interaction simple entre le sexe et le groupe d'âge</figcaption>
</figure>
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.
<figure>
```{r}
ggcoef_model(mod3, exponentiate = TRUE)
```
<figcaption>Représentation graphique des coefficients du modèle avec interaction simple entre le sexe et le groupe d'âge</figcaption>
</figure>
Cette fois-ci, il n'y a plus de coefficients globaux pour la variable <var>sexe</var> ni pour <var>grpage</var> mais des coefficients pour chaque combinaison de ces deux variables.
<figure>
```{r}
plot(
broken(mod3, nouvelle_observation, predict.function = betas),
trans = logit
) + ylim(0, 1) + ylab("Probabilité de faire du sport")
```
<figcaption>Représentation graphique de l'estimation de la probabilité de faire du sport pour un homme de 40 ans</figcaption>
</figure>
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 ?
```{r}
mod4 <- glm(sport ~ sexe * etud + grpage + heures.tv + relig, data = d, family = binomial())
```
Regardons d'abord les effets.
<figure>
```{r}
plot(allEffects(mod4))
```
<figcaption>Représentation graphique des effets du modèle avec interaction entre le sexe et le niveau d'étude</figcaption>
</figure>
À 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`{data-pkg="ggeffects"} de l'extension `ggeffects`{.pkg} qui permets de récupérer les effets calculés avec `effect`{data-pkg="effects"} dans un format utilisable avec `ggplot2`{.pkg}.
<figure>
```{r}
library(ggeffects)
plot(ggeffect(mod4, c("etud", "sexe")))
```
<figcaption>Effets du niveau d'étude selon le sexe</figcaption>
</figure>
Cela confirme ce que l'on suppose. Regardons les coefficients du modèle.
<figure>
```{r}
ggcoef_model(mod4, exponentiate = TRUE)
```
<figcaption>Représentation graphique des coefficients du modèle avec interaction simple entre le sexe et le niveau d'étude</figcaption>
</figure>
```{r}
tbl_regression(mod4, exponentiate = TRUE)
```
Si les coefficients associés au niveau d'étude sont significatifs, ceux de l'interaction ne le sont pas (sauf <q>sexeHomme:etudManquant</q>) et celui associé au sexe, précédemment significatif ne l'est plus. Testons avec `anova`{data-pkg="stats"} si l'interaction est belle et bien significative.
```{r}
anova(mod4, 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.
## 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`{data-pkg="ggmulti"} de l'extension du même nom. `glmulti`{data-pkg="ggmulti"} 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'<df>AIC</dfn>). ATTENTION : le temps de calcul de `glmulti`{data-pkg="ggmulti"} peut-être long.
```{r, eval=FALSE}
library(glmulti)
glmulti(sport ~ sexe + grpage + etud + heures.tv + relig, data = d, family = binomial())
```
```
Initialization...
TASK: Exhaustive screening of candidate set.
Fitting...
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
After 100 models:
Best model: sport~1+etud+heures.tv+grpage:heures.tv
Crit= 2267.79462883348
Mean crit= 2360.46497457747
After 150 models:
Best model: sport~1+grpage+etud+heures.tv+sexe:heures.tv
Crit= 2228.88574082404
Mean crit= 2286.60589884071
After 200 models:
Best model: sport~1+grpage+etud+heures.tv+sexe:heures.tv
Crit= 2228.88574082404
Mean crit= 2254.99359340075
After 250 models:
Best model: sport~1+sexe+grpage+etud+heures.tv+etud:sexe+sexe:heures.tv
Crit= 2226.00088609349
Mean crit= 2241.76611580481
After 300 models:
Best model: sport~1+sexe+grpage+etud+heures.tv+grpage:sexe+sexe:heures.tv
Crit= 2222.67161519005
Mean crit= 2234.95020358944
```
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.
```{r}
best <- glm(sport~1+sexe+grpage+etud+heures.tv+grpage:sexe+sexe:heures.tv, data = d, family = binomial())
odds.ratio(best)
```
<figure>
```{r}
ggcoef_model(best, exponentiate = TRUE)
```
<figcaption>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</figcaption>
</figure>
<figure>
```{r}
plot(allEffects(best))
```
<figcaption>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</figcaption>
</figure>
## 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.
On pourra en particulier se référer à la vignette inclue avec l'extension `phia`{.pkg} : <https://cran.r-project.org/web/packages/phia/vignettes/phia.pdf>.