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.

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

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())
odds.ratio(mod)
Waiting for profiling to be done...
                                      OR   2.5 %  97.5 %         p    
(Intercept)                      0.45006 0.23770  0.8473 0.0137076 *  
sexeHomme                        1.55223 1.26143  1.9120 3.389e-05 ***
grpage[25,45)                    0.65675 0.41944  1.0275 0.0652358 .  
grpage[45,65)                    0.33776 0.21153  0.5381 4.969e-06 ***
grpage[65,99]                    0.25124 0.14639  0.4288 4.531e-07 ***
etudSecondaire                   2.58719 1.76527  3.8328 1.476e-06 ***
etudTechnique/Professionnel      2.85552 1.98049  4.1729 3.237e-08 ***
etudSupérieur                    6.63042 4.55179  9.7951 < 2.2e-16 ***
etudManquant                     8.58853 4.53478 16.5885 7.419e-11 ***
heures.tv                        0.88611 0.82902  0.9458 0.0003188 ***
religPratiquant occasionnel      0.97833 0.67588  1.4199 0.9078331    
religAppartenance sans pratique  0.99333 0.70630  1.4020 0.9694337    
religNi croyance ni appartenance 0.80623 0.55242  1.1782 0.2646171    
religRejet                       0.68144 0.38688  1.1890 0.1797562    
religNSP ou NVPR                 0.91963 0.39967  2.0216 0.8384696    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

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")
Analysis of Deviance Table

Model: binomial, link: logit

Response: sport

Terms added sequentially (first to last)

            Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                         1994     2609.2              
sexe         1   16.211      1993     2593.0 5.668e-05 ***
grpage       3  216.268      1990     2376.7 < 2.2e-16 ***
etud         4  152.732      1986     2224.0 < 2.2e-16 ***
heures.tv    1   13.588      1985     2210.4 0.0002276 ***
relig        5    4.232      1980     2206.2 0.5165401    
sexe:grpage  3   13.030      1977     2193.1 0.0045714 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Jetons maintenant un oeil aux coefficients du modèle. Pour rendre les choses plus visuelles, nous aurons recours à ggcoef de l’extension GGally.

library(GGally)
ggcoef(mod2, exponentiate = TRUE)
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

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

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")
Analysis of Deviance Table

Model 1: sport ~ sexe * grpage + etud + heures.tv + relig
Model 2: sport ~ sexe:grpage + etud + heures.tv + relig
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1      1977     2193.1                     
2      1977     2193.1  0        0         

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

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 l’âge ?

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)
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(mod4, exponentiate = TRUE)
Représentation graphique des coefficients du modèle avec interaction simple entre le sexe et le niveau d’étude
odds.ratio(mod4)
Waiting for profiling to be done...
                                           OR   2.5 %  97.5 %         p    
(Intercept)                           0.45928 0.22781  0.9126 0.0276831 *  
sexeHomme                             1.41659 0.77234  2.5843 0.2558796    
etudSecondaire                        2.24794 1.35296  3.8167 0.0021242 ** 
etudTechnique/Professionnel           3.09244 1.89953  5.1709 9.260e-06 ***
etudSupérieur                         6.59073 4.01753 11.1205 3.307e-13 ***
etudManquant                          5.21288 2.42496 11.4294 2.847e-05 ***
grpage[25,45)                         0.65985 0.41925  1.0364 0.0710897 .  
grpage[45,65)                         0.34249 0.21353  0.5477 7.877e-06 ***
grpage[65,99]                         0.25402 0.14745  0.4350 6.606e-07 ***
heures.tv                             0.88707 0.82985  0.9468 0.0003646 ***
religPratiquant occasionnel           0.97621 0.67318  1.4196 0.8992576    
religAppartenance sans pratique       1.00874 0.71604  1.4266 0.9605105    
religNi croyance ni appartenance      0.82632 0.56489  1.2108 0.3262194    
religRejet                            0.70306 0.39797  1.2297 0.2200263    
religNSP ou NVPR                      0.89563 0.38588  1.9820 0.7904669    
sexeHomme:etudSecondaire              1.37693 0.65003  2.9290 0.4037975    
sexeHomme:etudTechnique/Professionnel 0.87434 0.43790  1.7538 0.7035417    
sexeHomme:etudSupérieur               1.00776 0.49244  2.0725 0.9831280    
sexeHomme:etudManquant                4.28305 1.32153 15.8357 0.0200885 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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")
Analysis of Deviance Table

Model: binomial, link: logit

Response: sport

Terms added sequentially (first to last)

          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                       1994     2609.2              
sexe       1    16.21      1993     2593.0 5.668e-05 ***
etud       4   319.19      1989     2273.8 < 2.2e-16 ***
grpage     3    49.81      1986     2224.0 8.770e-11 ***
heures.tv  1    13.59      1985     2210.4 0.0002276 ***
relig      5     4.23      1980     2206.2 0.5165401    
sexe:etud  4    10.22      1976     2195.9 0.0368395 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

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.

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...
                                  OR    2.5 %  97.5 %         p    
(Intercept)                 0.285395 0.144811  0.5557 0.0002508 ***
sexeHomme                   4.018732 1.808573  9.4145 0.0008949 ***
grpage[25,45)               1.030205 0.591312  1.8184 0.9170428    
grpage[45,65)               0.627599 0.351868  1.1323 0.1171763    
grpage[65,99]               0.493477 0.249106  0.9762 0.0421203 *  
etudSecondaire              2.541989 1.735573  3.7634 2.211e-06 ***
etudTechnique/Professionnel 2.817570 1.956251  4.1137 4.443e-08 ***
etudSupérieur               6.685085 4.588538  9.8796 < 2.2e-16 ***
etudManquant                8.817567 4.575016 17.3668 1.447e-10 ***
heures.tv                   0.845750 0.771519  0.9237 0.0002612 ***
sexeHomme:grpage[25,45)     0.318004 0.134434  0.7155 0.0069626 ** 
sexeHomme:grpage[45,65)     0.242528 0.101362  0.5524 0.0010034 ** 
sexeHomme:grpage[65,99]     0.224635 0.083753  0.5805 0.0024190 ** 
sexeHomme:heures.tv         1.115351 0.979397  1.2705 0.0997266 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggcoef(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.

On pourra en particulier se référer à la vignette inclue avec l’extension phia : https://cran.r-project.org/web/packages/phia/vignettes/phia.pdf.


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