guide-R - 26 Interactions

26.1 Données d’illustration

Reprenons le modèle que nous avons utilisé dans le chapitre sur la régression logistique binaire (cf. Chapitre 22).

library(tidyverse)
library(labelled)

data(hdv2003, package = "questionr")

d <-
  hdv2003 |> 
  mutate(
    sexe = sexe |> fct_relevel("Femme"),
    groupe_ages = age |>
      cut(
        c(18, 25, 45, 65, 99),
        right = FALSE,
        include.lowest = TRUE,
        labels = c("18-24 ans", "25-44 ans",
                   "45-64 ans", "65 ans et plus")
      ),
    etudes = nivetud |> 
      fct_recode(
        "Primaire" = "N'a jamais fait d'etudes",
        "Primaire" = "A arrete ses etudes, avant la derniere annee d'etudes primaires",
        "Primaire" = "Derniere annee d'etudes primaires",
        "Secondaire" = "1er cycle",
        "Secondaire" = "2eme cycle",
        "Technique / Professionnel" = "Enseignement technique ou professionnel court",
        "Technique / Professionnel" = "Enseignement technique ou professionnel long",
        "Supérieur" = "Enseignement superieur y compris technique superieur"
    ) |> 
    fct_na_value_to_level("Non documenté")  
  ) |> 
  set_variable_labels(
    sport = "Pratique un sport ?",
    sexe = "Sexe",
    groupe_ages = "Groupe d'âges",
    etudes = "Niveau d'études",
    heures.tv = "Heures de télévision / jour"
  )

26.2 Modèle sans interaction

Nous avions alors exploré les facteurs associés au fait de pratiquer du sport.

mod <- glm(
  sport ~ sexe + groupe_ages + etudes + heures.tv,
  family = binomial,
  data = d
)
library(gtsummary)
theme_gtsummary_language(
  "fr",
  decimal.mark = ",",
  big.mark = " "
)

mod |> 
  tbl_regression(exponentiate = TRUE) |> 
  bold_labels()

Caractéristique	OR ¹	95% IC ¹	p-valeur
Sexe
Femme	—	—
Homme	1,52	1,24 – 1,87	<0,001
Groupe d'âges
18-24 ans	—	—
25-44 ans	0,68	0,43 – 1,06	0,084
45-64 ans	0,36	0,23 – 0,57	<0,001
65 ans et plus	0,27	0,16 – 0,46	<0,001
Niveau d'études
Primaire	—	—
Secondaire	2,54	1,73 – 3,75	<0,001
Technique / Professionnel	2,81	1,95 – 4,10	<0,001
Supérieur	6,55	4,50 – 9,66	<0,001
Non documenté	8,54	4,51 – 16,5	<0,001
Heures de télévision / jour	0,89	0,83 – 0,95	<0,001
¹ OR = rapport de cotes, IC = intervalle de confiance

Table 26.1: Odds Ratios du modèle logistique simple

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.

Dans le chapitre sur les estimations marginales, cf. Chapitre 24, nous avons présenté la fonction broom.helpers::plot_marginal_predictions() qui permet de représenter les prédictions marginales moyennes du modèle.

mod |> 
  broom.helpers::plot_marginal_predictions(type = "response") |> 
  patchwork::wrap_plots() &
  scale_y_continuous(
    limits = c(0, .8),
    labels = scales::label_percent()
  )

Figure 26.1: Prédictions marginales moyennes du modèle simple

26.3 Définition d’une interaction

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 * groupe_agesdans l’équation du modèle.

mod2 <- glm(
  sport ~ sexe * groupe_ages + etudes + heures.tv,
  family = binomial,
  data = d
)

Commençons par regarder les prédictions marginales du modèle avec interaction.

mod2 |> 
  broom.helpers::plot_marginal_predictions(type = "response") |> 
  patchwork::wrap_plots(ncol = 1) &
  scale_y_continuous(
    labels = scales::label_percent()
  )

Figure 26.2: Prédictions marginales moyennes du modèle avec interaction

Sur ce graphique, on voit que la pratique d’un sport diminue fortement avec l’âge chez les hommes, tandis que cette diminution est bien plus modérée chez les femmes.

Astuce

Par défaut, broom.helpers::plot_marginal_predictions() détecte la présence d’interactions dans le modèle et calcule les prédictions marginales pour chaque combinaison de variables incluent dans une interaction. Il reste possible de calculer des prédictions marginales individuellement pour chaque variable du modèle. Pour cela, il suffit d’indiquer variables_list = "no_interaction".

mod2 |> 
  broom.helpers::plot_marginal_predictions(
    variables_list = "no_interaction",
    type = "response"
  ) |> 
  patchwork::wrap_plots() &
  scale_y_continuous(
    labels = scales::label_percent()
  )

26.4 Significativité de l’interaction

L’ajout d’une interaction au modèle augmente la capacité prédictive du modèle mais, dans le même temps, augmente le nombre de coefficients (et donc de degrés de liberté). La question se pose donc de savoir si l’ajout d’un terme d’interaction améliore notre modèle.

En premier lieu, nous pouvons comparer les AIC des modèles avec et sans interaction.

AIC(mod)

[1] 2230.404

AIC(mod2)

[1] 2223.382

L’AIC du modèle avec interaction est plus faible que celui sans interaction, nous indiquant un gain : notre modèle avec interaction est donc meilleur.

On peut tester avec car::Anova() si l’interaction est statistiquement significative¹.

¹ Lorsqu’il y a une interaction, il est préférable d’utiliser le type III, cf. Section 22.8. En toute rigueur, il serait préférable de coder nos variables catégorielles avec un contraste de type somme (cf. Chapitre 25). En pratique, nous pouvons nous en passer ici.

car::Anova(mod2, type = "III")

Analysis of Deviance Table (Type III tests)

Response: sport
                 LR Chisq Df Pr(>Chisq)    
sexe               19.349  1  1.089e-05 ***
groupe_ages        15.125  3  0.0017131 ** 
etudes            125.575  4  < 2.2e-16 ***
heures.tv          12.847  1  0.0003381 ***
sexe:groupe_ages   13.023  3  0.0045881 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La p-valeur associée au terme d’interaction (sexe:groupe_ages) est inférieure à 1% : l’interaction a donc bien un effet significatif.

Nous pouvons également utiliser gtsummary::add_global_p().

mod2 |> 
  tbl_regression(exponentiate = TRUE) |> 
  add_global_p() |> 
  bold_labels()

Caractéristique	OR ¹	95% IC ¹	p-valeur
Sexe			<0,001
Femme	—	—
Homme	5,10	2,41 – 11,4
Groupe d'âges			0,002
18-24 ans	—	—
25-44 ans	1,06	0,61 – 1,85
45-64 ans	0,64	0,36 – 1,15
65 ans et plus	0,49	0,25 – 0,97
Niveau d'études			<0,001
Primaire	—	—
Secondaire	2,55	1,74 – 3,78
Technique / Professionnel	2,84	1,97 – 4,14
Supérieur	6,69	4,60 – 9,89
Non documenté	8,94	4,64 – 17,6
Heures de télévision / jour	0,89	0,83 – 0,95	<0,001
Sexe * Groupe d'âges			0,005
Homme * 25-44 ans	0,31	0,13 – 0,70
Homme * 45-64 ans	0,24	0,10 – 0,54
Homme * 65 ans et plus	0,23	0,09 – 0,60
¹ OR = rapport de cotes, IC = intervalle de confiance

Table 26.2: Odds Ratios du modèle logistique avec interaction

26.5 Interprétation des coefficients

Jetons maintenant un œil aux coefficients du modèle. Pour rendre les choses plus visuelles, nous aurons recours à ggtstats::ggcoef_model().

mod2 |> 
  ggstats::ggcoef_model(exponentiate = TRUE)

Figure 26.3: Coefficients (odds ratio) du modèle avec interaction

Concernant les variables sexe et groupe_ages, nous avons trois séries de coefficients : une série pour le sexe, une pour le groupe d’âges et enfin des coefficients pour l’interaction entre le sexe et le groupe d’âges.

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 les autres variables correspondent aux modalités de référence (i.e. 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, grâce au package breakDown.

library(breakDown)
logit <- function(x) exp(x)/(1+exp(x))
nouvelle_observation <- d[1, ]
nouvelle_observation$sexe[1] = "Femme"
nouvelle_observation$groupe_ages[1] = "45-64 ans"
nouvelle_observation$etud[1] = "Primaire"
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")

Figure 26.4: 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 à la 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-64 ans de la variable groupe_ages.

Regardons maintenant la situation d’un homme de 20 ans.

nouvelle_observation$sexe[1] = "Homme"
nouvelle_observation$groupe_ages[1] = "18-24 ans"
plot(
  broken(mod2, nouvelle_observation, predict.function = betas),
  trans = logit
) +
  ylim(0, 1.2) +
  ylab("Probabilité de faire du sport")

Figure 26.5: 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.

Regardons enfin la situation d’un homme de 60 ans.

nouvelle_observation$groupe_ages[1] = "45-64 ans"
plot(
  broken(mod2, nouvelle_observation, predict.function = betas),
  trans = logit
) +
  ylim(0, 1.2) +
  ylab("Probabilité de faire du sport")

Figure 26.6: Représentation graphique de l’estimation de la probabilité de faire du sport pour un homme de 60 ans

Cette fois, plusieurs coefficients s’appliquent : à la fois le coefficient sexe = Homme (effet du sexe pour les 18-24 ans), le coefficient groupe_ages = 45-64 ans qui est l’effet de l’âge pour les femmes de 45-64 ans par rapport aux 18-24 ans et le coefficient sexe:groupe_ages = Homme:45-64 ans qui indique l’effet spécifique qui s’applique aux hommes de 45-64 ans, d’une part par rapport aux femmes du même âge et d’autre part par rapport aux hommes de 18-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.

26.6 Définition alternative de l’interaction

Il est cependant possible d’écrire le même modèle différemment. En effet, sexe * groupe_ages dans la formule du modèle est équivalent à l’écriture sexe + groupe_ages + sexe:groupe_ages, c’est-à-dire que l’on demande des coefficients pour la variable sexe à la référence de groupe_ages, des coefficients pour groupe_ages à la référence de sexe et enfin des coefficients pour tenir compte de l’interaction.

On peut se contenter d’une série de coefficients uniques pour l’interaction en indiquant seulement sexe : groupe_ages.

mod3 <- glm(
  sport ~ sexe : groupe_ages + etudes + heures.tv,
  family = binomial,
  data = d
)

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 stats::anova().

anova(mod2, mod3, test = "Chisq")

Analysis of Deviance Table

Model 1: sport ~ sexe * groupe_ages + etudes + heures.tv
Model 2: sport ~ sexe:groupe_ages + etudes + heures.tv
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1      1982     2197.4                     
2      1982     2197.4  0        0

De même, les prédictions marginales sont les mêmes, comme nous pouvons le constater avec ggstats::ggcoef_compare().

ggstats::ggcoef_compare(
  list("sexe * groupe_ages" = mod2, "sexe : groupe_ages" = mod3),
  tidy_fun = broom.helpers::tidy_marginal_predictions,
  significance = NULL,
  vline = FALSE
) +
  scale_x_continuous(labels = scales::label_percent())

Warning: Model matrix is rank deficient. Some variance-covariance parameters are
  missing.
Warning: Model matrix is rank deficient. Some variance-covariance parameters are
  missing.
Warning: Model matrix is rank deficient. Some variance-covariance parameters are
  missing.

Figure 26.7: Comparaison des prédictions marginales moyennes des deux modèles avec interaction

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.

mod3 |> 
  ggstats::ggcoef_model(exponentiate = TRUE)

Figure 26.8: Coefficients (odds ratio) du modèle avec interaction simple entre le sexe et le groupe d’âges

Cette fois-ci, il n’y a plus de coefficients globaux pour la variable sexe ni pour groupe_ages mais des coefficients pour chaque combinaison de ces deux variables. Reprenons l’exemple de notre homme de 60 ans.

plot(
  broken(mod3, nouvelle_observation, predict.function = betas),
  trans = logit
) +
  ylim(0, 1.2) +
  ylab("Probabilité de faire du sport")

Figure 26.9: Représentation graphique de l’estimation de la probabilité de faire du sport pour un homme de 60 ans (interaction simple)

Cette fois-ci, le coefficient d’interaction fournit indique l’effet combiné du sexe et du groupe d’âges par rapport à la situation de référence (femme de 18-24 ans).

Que l’on définisse une interaction simple (sexe:groupe_ages) ou complète (sexe*groupe_ages), les deux modèles calculés sont donc identiques en termes prédictifs et explicatifs, mais l’interprétation de leurs coefficients diffèrent.

26.7 Identifier les interactions pertinentes

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 toujours bon, selon notre connaissance du sujet et de la littérature, d’explorer manuellement les interactions attendues / prévisibles.

Mais, il est tentant de vouloir tester les multiples interactions possibles de manière itératives afin d’identifier celles à retenir.

Une possibilité² est d’avoir recours à une sélection de modèle pas à pas ascendante (voir Chapitre 23). Nous allons partir de notre modèle sans interaction, indiquer à step() l’ensemble des interactions possibles et voir si nous pouvons minimiser l’AIC.

² On pourra également regarder du côté de glmulti::glmulti() pour des approches alternatives.

mod4 <- mod |> 
  step(scope = list(upper = ~ sexe * groupe_ages * etudes * heures.tv))

Start:  AIC=2230.4
sport ~ sexe + groupe_ages + etudes + heures.tv

                        Df Deviance    AIC
+ sexe:groupe_ages       3   2197.4 2223.4
+ sexe:etudes            4   2199.6 2227.6
+ sexe:heures.tv         1   2207.6 2229.6
<none>                       2210.4 2230.4
+ groupe_ages:heures.tv  3   2207.0 2233.0
+ etudes:heures.tv       4   2207.4 2235.4
+ groupe_ages:etudes    11   2194.6 2236.6
- heures.tv              1   2224.0 2242.0
- sexe                   1   2226.4 2244.4
- groupe_ages            3   2260.6 2274.6
- etudes                 4   2334.3 2346.3

Step:  AIC=2223.38
sport ~ sexe + groupe_ages + etudes + heures.tv + sexe:groupe_ages

                        Df Deviance    AIC
+ sexe:heures.tv         1   2194.7 2222.7
<none>                       2197.4 2223.4
+ groupe_ages:heures.tv  3   2193.5 2225.5
+ sexe:etudes            4   2192.1 2226.1
+ etudes:heures.tv       4   2194.6 2228.6
- sexe:groupe_ages       3   2210.4 2230.4
+ groupe_ages:etudes    11   2183.1 2231.1
- heures.tv              1   2210.2 2234.2
- etudes                 4   2323.0 2341.0

Step:  AIC=2222.67
sport ~ sexe + groupe_ages + etudes + heures.tv + sexe:groupe_ages + 
    sexe:heures.tv

                        Df Deviance    AIC
<none>                       2194.7 2222.7
- sexe:heures.tv         1   2197.4 2223.4
+ groupe_ages:heures.tv  3   2190.4 2224.4
+ sexe:etudes            4   2189.0 2225.0
+ etudes:heures.tv       4   2191.6 2227.6
- sexe:groupe_ages       3   2207.6 2229.6
+ groupe_ages:etudes    11   2180.7 2230.7
- etudes                 4   2319.9 2339.9

mod4$formula

sport ~ sexe + groupe_ages + etudes + heures.tv + sexe:groupe_ages + 
    sexe:heures.tv

Le modèle final suggéré comprends une interaction entre le sexe et le groupe d’âges et une interaction entre le sexe et le nombre quotidien d’heures de télévision. Nous pouvons utiliser broom.helpers::plot_marginal_predictions() pour visualiser l’effet de ces deux interactions.

mod4 |> 
  broom.helpers::plot_marginal_predictions(type = "response") |> 
  patchwork::wrap_plots(ncol = 1) &
  scale_y_continuous(
    labels = scales::label_percent()
  )

Figure 26.10: Prédictions marginales moyennes du modèle avec deux interactions

26.8 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.

Les effets d’interaction par Jean-François Bickel
Analysing interactions of fitted models par Helios De Rosario Martínez

26.9 webin-R

Les interactions sont abordées dans le webin-R #07 (régression logistique partie 2) sur YouTube.