La version originale de ce chapitre a été écrite par Joseph Larmarange dans le cadre du support de cours Introduction à l’analyse d’enquêtes avec R.

La régression logistique est fréquemment utilisée en sciences sociales car elle permet d’effectuer un raisonnement dit toutes choses étant égales par ailleurs. Plus précisément, la régression logistique a pour but d’isoler les effets de chaque variable, c’est-à-dire d’identifier les effets résiduels d’une variable explicative sur une variable d’intérêt, une fois pris en compte les autres variables explicatives introduites dans le modèle. La régression logistique est ainsi prisée en épidémiologie pour identifier les facteurs associés à telle ou telle pathologie.

La régression logistique ordinaire ou régression logistique binaire vise à expliquer une variable d’intérêt binaire (c’est-à-dire de type « oui / non » ou « vrai / faux »). Les variables explicatives qui seront introduites dans le modèle peuvent être quantitatives ou qualitatives.

La régression logistique multinomiale est une extension de la régression logistique aux variables qualitatives à trois modalités ou plus.

Préparation des données

Dans ce chapite, nous allons encore une fois utiliser les données de l’enquête Histoire de vie, fournies avec l’extension questionr.

library(questionr)
data(hdv2003)
d <- hdv2003

À titre d’exemple, nous allons étudier l’effet de l’âge, du sexe, du niveau d’étude, de la pratique religieuse et du nombre moyen d’heures passées à regarder la télévision par jour.

En premier lieu, il importe de vérifier que notre variable d’intérêt (ici sport) est correctement codée. Une possibilité consiste à créer une variable booléenne (vrai / faux) selon que l’individu a pratiqué du sport ou non :

d$sport2 <- FALSE
d$sport2[d$sport == "Oui"] <- TRUE

Dans le cas présent, cette variable n’a pas de valeur manquante. Mais, le cas échéant, il importe de bien coder les valeurs manquantes en NA, les individus en question étant alors exclu de l’analyse.

Il n’est pas forcément nécessaire de transformer notre variable d’intérêt en variable booléenne. En effet, R accepte sans problème une variable de type facteur. Cependant, l’ordre des valeurs d’un facteur a de l’importance. En effet, R considère toujours la première modalité comme étant la modalité de référence. Dans le cas de la variable d’intérêt, la modalité de référence correspond au fait de ne pas remplir le critère étudié, dans notre exemple au fait de ne pas avoir eu d’activité sportive au cours des douze derniers mois.

Pour connaître l’ordre des modalités d’une variable de type facteur, on peut utiliser la fonction levels ou bien encore tout simplement la fonction freq de l’extension questionr :

levels(d$sport)
[1] "Non" "Oui"
freq(d$sport)
       n    % val%
Non 1277 63.8 63.8
Oui  723 36.1 36.1

Dans notre exemple, la modalité « Non » est déjà la première modalité. Il n’y a donc pas besoin de modifier notre variable. Si ce n’est pas le cas, il faudra modifier la modalité de référence avec la fonction relevel comme nous allons le voir un peu plus loin.

Il est possible d’indiquer un facteur à plus de deux modalités. Dans une telle situation, R considérera que tous les modalités, sauf la modalité de référence, est une réalisation de la variable d’intérêt. Cela serait correct, par exemple, si notre variable sport était codée ainsi : « Non », « Oui, toutes les semaines », « Oui, au moins une fois par mois », « Oui, moins d’une fois par mois ». Cependant, afin d’éviter tout risque d’erreur ou de mauvaise interprétation, il est vivement conseillé de recoder au préalable sa variable d’intérêt en un facteur à deux modalités.

La notion de modalité de référence s’applique également aux variables explicatives qualitatives. En effet, dans un modèle, tous les coefficients sont calculés par rapport à la modalité de référence. Il importe de choisir une modalité de référence qui fasse sens afin de faciliter l’interprétation. Par ailleurs, ce choix peut également dépendre de la manière dont on souhaite présenter les résultats. De manière générale on évitera de choisir comme référence une modalité peu représentée dans l’échantillon ou bien une modalité correspondant à une situation atypique.

Prenons l’exemple de la variable sexe. Souhaite-t-on connaitre l’effet d’être une femme par rapport au fait d’être un homme ou bien l’effet d’être un homme par rapport au fait d’être une femme ? Si l’on opte pour le second, alors notre modalité de référence sera le sexe féminin. Comme est codée cette variable ?

freq(d$sexe)
         n  % val%
Homme  899 45   45
Femme 1101 55   55

La modalité « Femme » s’avère ne pas être la première modalité. Nous devons appliquer la fonction relevel :

d$sexe <- relevel(d$sexe, "Femme")
freq(d$sexe)
         n  % val%
Femme 1101 55   55
Homme  899 45   45

Données labellisées

Si l’on utilise des données labellisées (voir le chapitre dédié), nos variables catégorielles seront stockées sous la forme d’un vecteur numérique avec des étiquettes. Il sera donc nécessaire de convertir ces variables en facteurs, tout simplement avec la fonction to_factor de l’extension labelled qui pourra utiliser les étiquettes de valeurs comme modalités du facteur.

Les variables age et heures.tv sont des variables quantitatives. Il importe de vérifier qu’elles sont bien enregistrées en tant que variables numériques. En effet, il arrive parfois que dans le fichier source les variables quantitatives soient renseignées sous forme de valeur textuelle et non sous forme numérique.

str(d$age)
 int [1:2000] 28 23 59 34 71 35 60 47 20 28 ...
str(d$heures.tv)
 num [1:2000] 0 1 0 2 3 2 2.9 1 2 2 ...

Nos deux variables sont bien renseignées sous forme numérique.

Cependant, l’effet de l’âge est rarement linéaire. Un exemple trivial est par exemple le fait d’occuper un emploi qui sera moins fréquent aux jeunes âges et aux âges élevés. Dès lors, on pourra transformer la variable age en groupe d’âges avec la fonction cut (voir le chapitre Manipulation de données) :

d$grpage <- cut(d$age, c(16, 25, 45, 65, 93), right = FALSE, 
  include.lowest = TRUE)
freq(d$grpage)
          n    % val%
[16,25) 169  8.5  8.5
[25,45) 706 35.3 35.3
[45,65) 745 37.2 37.3
[65,93] 378 18.9 18.9
NA        2  0.1   NA

Jetons maintenant un oeil à la variable nivetud :

freq(d$nivetud)
                                                                  n
N'a jamais fait d'etudes                                         39
A arrete ses etudes, avant la derniere annee d'etudes primaires  86
Derniere annee d'etudes primaires                               341
1er cycle                                                       204
2eme cycle                                                      183
Enseignement technique ou professionnel court                   463
Enseignement technique ou professionnel long                    131
Enseignement superieur y compris technique superieur            441
NA                                                              112
                                                                   %
N'a jamais fait d'etudes                                         2.0
A arrete ses etudes, avant la derniere annee d'etudes primaires  4.3
Derniere annee d'etudes primaires                               17.1
1er cycle                                                       10.2
2eme cycle                                                       9.2
Enseignement technique ou professionnel court                   23.2
Enseignement technique ou professionnel long                     6.6
Enseignement superieur y compris technique superieur            22.1
NA                                                               5.6
                                                                val%
N'a jamais fait d'etudes                                         2.1
A arrete ses etudes, avant la derniere annee d'etudes primaires  4.6
Derniere annee d'etudes primaires                               18.1
1er cycle                                                       10.8
2eme cycle                                                       9.7
Enseignement technique ou professionnel court                   24.5
Enseignement technique ou professionnel long                     6.9
Enseignement superieur y compris technique superieur            23.4
NA                                                                NA

En premier lieu, cette variable est détaillée en pas moins de huit modalités dont certaines sont peu représentées (seulement 39 individus soit 2 % n’ont jamais fait d’études par exemple). Afin d’améliorier notre modèle logistique, il peut être pertinent de regrouper certaines modalités (voir le chapitre Manipulation de données) :

d$etud <- d$nivetud
levels(d$etud) <- c("Primaire", "Primaire", "Primaire", 
  "Secondaire", "Secondaire", "Technique/Professionnel", 
  "Technique/Professionnel", "Supérieur")
freq(d$etud)
                          n    % val%
Primaire                466 23.3 24.7
Secondaire              387 19.4 20.5
Technique/Professionnel 594 29.7 31.5
Supérieur               441 22.1 23.4
NA                      112  5.6   NA

Notre variable comporte également 112 individus avec une valeur manquante. Si nous conservons cette valeur manquante, ces 112 individus seront, par défaut, exclus de l’analyse. Ces valeurs manquantes n’étant pas négligeable (5,6 %), nous pouvons également faire le choix de considérer ces valeurs manquantes comme une modalité supplémentaire. Auquel cas, nous utiliserons la fonction add.NA :

levels(d$etud)
[1] "Primaire"               
[2] "Secondaire"             
[3] "Technique/Professionnel"
[4] "Supérieur"              
d$etud <- addNA(d$etud)
levels(d$etud)
[1] "Primaire"               
[2] "Secondaire"             
[3] "Technique/Professionnel"
[4] "Supérieur"              
[5] NA                       

Régression logistique binaire

La fonction glm (pour generalized linear models soit modèle linéaire généralisé en français) permet de calculer une grande variété de modèles statistiques. La régression logistique ordinaire correspond au modèle logit de la famille des modèles binomiaux, ce que l’on indique à glm avec l’argument family=binomial(logit).

Le modèle proprement dit sera renseigné sous la forme d’une formule (que nous avons déjà rencontrée dans le chapitre sur la statistique bivariée et présentée plus en détails dans un chapitre dédié). On indiquera d’abord la variable d’intérêt, suivie du signe ~ (que l’on obtient en appuyant sur les touches Alt Gr et 3 sur un clavier de type PC) puis de la liste des variables explicatives séparées par un signe +. Enfin, l’argument data permettra d’indiquer notre tableau de données.

reg <- glm(sport ~ sexe + grpage + etud + relig + heures.tv, 
  data = d, family = binomial(logit))
reg

Call:  glm(formula = sport ~ sexe + grpage + etud + relig + heures.tv, 
    family = binomial(logit), data = d)

Coefficients:
                     (Intercept)  
                       -0.797364  
                       sexeHomme  
                        0.439002  
                   grpage[25,45)  
                       -0.420306  
                   grpage[45,65)  
                       -1.085463  
                   grpage[65,93]  
                       -1.379274  
                  etudSecondaire  
                        0.948312  
     etudTechnique/Professionnel  
                        1.047156  
                   etudSupérieur  
                        1.889621  
                          etudNA  
                        2.148545  
     religPratiquant occasionnel  
                       -0.020597  
 religAppartenance sans pratique  
                       -0.006176  
religNi croyance ni appartenance  
                       -0.214067  
                      religRejet  
                       -0.382744  
                religNSP ou NVPR  
                       -0.083361  
                       heures.tv  
                       -0.120724  

Degrees of Freedom: 1992 Total (i.e. Null);  1978 Residual
  (7 observations deleted due to missingness)
Null Deviance:      2607 
Residual Deviance: 2206     AIC: 2236

Il est possible de spécifier des modèles plus complexes. Par exemple, x:y permet d’indiquer l’interaction entre les variables x et y. x * y sera équivalent à x + y + x:y. Pour aller plus loin, voir http://ww2.coastal.edu/kingw/statistics/R-tutorials/formulae.html.

Une présentation plus complète des résultats est obtenue avec la méthode summary :

summary(reg)

Call:
glm(formula = sport ~ sexe + grpage + etud + relig + heures.tv, 
    family = binomial(logit), data = d)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8783  -0.8864  -0.4814   1.0033   2.4202  

Coefficients:
                                  Estimate
(Intercept)                      -0.797364
sexeHomme                         0.439002
grpage[25,45)                    -0.420306
grpage[45,65)                    -1.085463
grpage[65,93]                    -1.379274
etudSecondaire                    0.948312
etudTechnique/Professionnel       1.047156
etudSupérieur                     1.889621
etudNA                            2.148545
religPratiquant occasionnel      -0.020597
religAppartenance sans pratique  -0.006176
religNi croyance ni appartenance -0.214067
religRejet                       -0.382744
religNSP ou NVPR                 -0.083361
heures.tv                        -0.120724
                                 Std. Error
(Intercept)                        0.323833
sexeHomme                          0.106063
grpage[25,45)                      0.228042
grpage[45,65)                      0.237704
grpage[65,93]                      0.273794
etudSecondaire                     0.197427
etudTechnique/Professionnel        0.189777
etudSupérieur                      0.195190
etudNA                             0.330198
religPratiquant occasionnel        0.189203
religAppartenance sans pratique    0.174709
religNi croyance ni appartenance   0.193096
religRejet                         0.285878
religNSP ou NVPR                   0.410968
heures.tv                          0.033589
                                 z value Pr(>|z|)
(Intercept)                       -2.462 0.013806
sexeHomme                          4.139 3.49e-05
grpage[25,45)                     -1.843 0.065313
grpage[45,65)                     -4.566 4.96e-06
grpage[65,93]                     -5.038 4.71e-07
etudSecondaire                     4.803 1.56e-06
etudTechnique/Professionnel        5.518 3.43e-08
etudSupérieur                      9.681  < 2e-16
etudNA                             6.507 7.67e-11
religPratiquant occasionnel       -0.109 0.913311
religAppartenance sans pratique   -0.035 0.971801
religNi croyance ni appartenance  -1.109 0.267600
religRejet                        -1.339 0.180624
religNSP ou NVPR                  -0.203 0.839260
heures.tv                         -3.594 0.000325
                                    
(Intercept)                      *  
sexeHomme                        ***
grpage[25,45)                    .  
grpage[45,65)                    ***
grpage[65,93]                    ***
etudSecondaire                   ***
etudTechnique/Professionnel      ***
etudSupérieur                    ***
etudNA                           ***
religPratiquant occasionnel         
religAppartenance sans pratique     
religNi croyance ni appartenance    
religRejet                          
religNSP ou NVPR                    
heures.tv                        ***
---
Signif. codes:  
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2607.4  on 1992  degrees of freedom
Residual deviance: 2205.9  on 1978  degrees of freedom
  (7 observations deleted due to missingness)
AIC: 2235.9

Number of Fisher Scoring iterations: 4

Dans le cadre d’un modèle logistique, généralement on ne présente pas les coefficients du modèle mais leur valeur exponentielle, cette dernière correspondant en effet à des odds ratio, également appelés rapports des cotes. L’odds ratio diffère du risque relatif. Cependent son interprétation est similaire. Un odds ratio de 1 signifie l’absence d’effet. Un odds ratio largement supérieur à 1 correspond à une augmentation du phénomène étudié et un odds ratio largement inféieur à 1 correspond à une diminution du phénomène étudié1.

La fonction coef permet d’obtenir les coefficients d’un modèle, confint leurs intervalles de confiance et exp de calculer l’exponentiel. Les odds ratio et leurs intervalles de confiance s’obtiennent ainsi :

exp(coef(reg))
                     (Intercept) 
                       0.4505151 
                       sexeHomme 
                       1.5511577 
                   grpage[25,45) 
                       0.6568456 
                   grpage[45,65) 
                       0.3377455 
                   grpage[65,93] 
                       0.2517614 
                  etudSecondaire 
                       2.5813476 
     etudTechnique/Professionnel 
                       2.8495358 
                   etudSupérieur 
                       6.6168632 
                          etudNA 
                       8.5723745 
     religPratiquant occasionnel 
                       0.9796133 
 religAppartenance sans pratique 
                       0.9938431 
religNi croyance ni appartenance 
                       0.8072940 
                      religRejet 
                       0.6819874 
                religNSP ou NVPR 
                       0.9200190 
                       heures.tv 
                       0.8862785 
exp(confint(reg))
Waiting for profiling to be done...
                                     2.5 %
(Intercept)                      0.2379725
sexeHomme                        1.2605519
grpage[25,45)                    0.4195082
grpage[45,65)                    0.2115297
grpage[65,93]                    0.1466915
etudSecondaire                   1.7613423
etudTechnique/Professionnel      1.9764492
etudSupérieur                    4.5427464
etudNA                           4.5265421
religPratiquant occasionnel      0.6767575
religAppartenance sans pratique  0.7067055
religNi croyance ni appartenance 0.5531366
religRejet                       0.3872078
religNSP ou NVPR                 0.3998878
heures.tv                        0.8291800
                                     97.5 %
(Intercept)                       0.8480538
sexeHomme                         1.9106855
grpage[25,45)                     1.0275786
grpage[45,65)                     0.5380619
grpage[65,93]                     0.4296919
etudSecondaire                    3.8240276
etudTechnique/Professionnel       4.1639745
etudSupérieur                     9.7745023
etudNA                           16.5563346
religPratiquant occasionnel       1.4217179
religAppartenance sans pratique   1.4026763
religNi croyance ni appartenance  1.1798422
religRejet                        1.1898605
religNSP ou NVPR                  2.0221544
heures.tv                         0.9459484

On pourra faciliter la lecture en combinant les deux :

exp(cbind(coef(reg), confint(reg)))
Waiting for profiling to be done...
                                          
(Intercept)                      0.4505151
sexeHomme                        1.5511577
grpage[25,45)                    0.6568456
grpage[45,65)                    0.3377455
grpage[65,93]                    0.2517614
etudSecondaire                   2.5813476
etudTechnique/Professionnel      2.8495358
etudSupérieur                    6.6168632
etudNA                           8.5723745
religPratiquant occasionnel      0.9796133
religAppartenance sans pratique  0.9938431
religNi croyance ni appartenance 0.8072940
religRejet                       0.6819874
religNSP ou NVPR                 0.9200190
heures.tv                        0.8862785
                                     2.5 %
(Intercept)                      0.2379725
sexeHomme                        1.2605519
grpage[25,45)                    0.4195082
grpage[45,65)                    0.2115297
grpage[65,93]                    0.1466915
etudSecondaire                   1.7613423
etudTechnique/Professionnel      1.9764492
etudSupérieur                    4.5427464
etudNA                           4.5265421
religPratiquant occasionnel      0.6767575
religAppartenance sans pratique  0.7067055
religNi croyance ni appartenance 0.5531366
religRejet                       0.3872078
religNSP ou NVPR                 0.3998878
heures.tv                        0.8291800
                                     97.5 %
(Intercept)                       0.8480538
sexeHomme                         1.9106855
grpage[25,45)                     1.0275786
grpage[45,65)                     0.5380619
grpage[65,93]                     0.4296919
etudSecondaire                    3.8240276
etudTechnique/Professionnel       4.1639745
etudSupérieur                     9.7745023
etudNA                           16.5563346
religPratiquant occasionnel       1.4217179
religAppartenance sans pratique   1.4026763
religNi croyance ni appartenance  1.1798422
religRejet                        1.1898605
religNSP ou NVPR                  2.0221544
heures.tv                         0.9459484

Pour savoir si un odds ratio diffère significativement de 1 (ce qui est identique au fait que le coefficient soit différent de 0), on pourra se référer à la colonne Pr(>|z|) obtenue avec summary.

Si vous disposez de l’extension questionr, la fonction odds.ratio permet de calculer directement les odds ratio, leur intervalles de confiance et les p-value :

library(questionr)
odds.ratio(reg)
Waiting for profiling to be done...
                                      OR   2.5 %
(Intercept)                      0.45052 0.23797
sexeHomme                        1.55116 1.26055
grpage[25,45)                    0.65685 0.41951
grpage[45,65)                    0.33775 0.21153
grpage[65,93]                    0.25176 0.14669
etudSecondaire                   2.58135 1.76134
etudTechnique/Professionnel      2.84954 1.97645
etudSupérieur                    6.61686 4.54275
etudNA                           8.57237 4.52654
religPratiquant occasionnel      0.97961 0.67676
religAppartenance sans pratique  0.99384 0.70671
religNi croyance ni appartenance 0.80729 0.55314
religRejet                       0.68199 0.38721
religNSP ou NVPR                 0.92002 0.39989
heures.tv                        0.88628 0.82918
                                  97.5 %
(Intercept)                       0.8481
sexeHomme                         1.9107
grpage[25,45)                     1.0276
grpage[45,65)                     0.5381
grpage[65,93]                     0.4297
etudSecondaire                    3.8240
etudTechnique/Professionnel       4.1640
etudSupérieur                     9.7745
etudNA                           16.5563
religPratiquant occasionnel       1.4217
religAppartenance sans pratique   1.4027
religNi croyance ni appartenance  1.1798
religRejet                        1.1899
religNSP ou NVPR                  2.0222
heures.tv                         0.9459
                                         p    
(Intercept)                      0.0138062 *  
sexeHomme                        3.488e-05 ***
grpage[25,45)                    0.0653132 .  
grpage[45,65)                    4.961e-06 ***
grpage[65,93]                    4.713e-07 ***
etudSecondaire                   1.560e-06 ***
etudTechnique/Professionnel      3.432e-08 ***
etudSupérieur                    < 2.2e-16 ***
etudNA                           7.674e-11 ***
religPratiquant occasionnel      0.9133105    
religAppartenance sans pratique  0.9718010    
religNi croyance ni appartenance 0.2675998    
religRejet                       0.1806239    
religNSP ou NVPR                 0.8392596    
heures.tv                        0.0003255 ***
---
Signif. codes:  
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Il est possible de représenter graphiquement les différents odds ratios. Pour cela, on va utiliser la fonction tidy de l’extension broom pour récupérer les coefficients du modèle sous la forme d’un tableau de données exploitable avec ggplot2. On précisera conf.int = TRUE pour obtenir les intervalles de confiance et exponentiate = TRUE pour avoir les odds ratio plutôt que les coefficients bruts. geom_errorbarh permets de représenter les intervalles de confiance sous forme de barres d’erreurs, geom_vline une ligne verticale au niveau x = 1, scale_x_log10 pour afficher l’axe des x de manière logarithmique, les odds ratios étant de nature multiplicative et non additive.

library(broom)
tmp <- tidy(reg, conf.int = TRUE, exponentiate = TRUE)
str(tmp)
'data.frame':   15 obs. of  7 variables:
 $ term     : chr  "(Intercept)" "sexeHomme" "grpage[25,45)" "grpage[45,65)" ...
 $ estimate : num  0.451 1.551 0.657 0.338 0.252 ...
 $ std.error: num  0.324 0.106 0.228 0.238 0.274 ...
 $ statistic: num  -2.46 4.14 -1.84 -4.57 -5.04 ...
 $ p.value  : num  1.38e-02 3.49e-05 6.53e-02 4.96e-06 4.71e-07 ...
 $ conf.low : num  0.238 1.261 0.42 0.212 0.147 ...
 $ conf.high: num  0.848 1.911 1.028 0.538 0.43 ...
library(ggplot2)
ggplot(tmp) + aes(x = estimate, y = term, xmin = conf.low, 
  xmax = conf.high) + geom_vline(xintercept = 1) + 
  geom_errorbarh() + geom_point() + scale_x_log10()
Représentation graphique des odds ratios

La fonction ggcoef de l’extension GGally permet d’effectuer le graphique précédent directement à partir de notre modèle. Voir l’aide de cette fonction pour la liste complète des paramètres personnalisables.

library(GGally)
ggcoef(reg, exponentiate = TRUE)
La fonction ggcoef

L’extension JLutils, disponible uniquement sur GitHub, propose une fonction intéressante dans ce contexte. Pour l’installer ou la mettre à jour, si ce n’est déjà fait, on aura recours à la commande ci-après.

devtools::install_github("larmarange/JLutils")

La fonction tidy_detailed est une version élargie de tidy qui rajoute des colonnes supplémentaires avec le nom des variables et des modalités.

str(tidy(reg))
'data.frame':   15 obs. of  5 variables:
 $ term     : chr  "(Intercept)" "sexeHomme" "grpage[25,45)" "grpage[45,65)" ...
 $ estimate : num  -0.797 0.439 -0.42 -1.085 -1.379 ...
 $ std.error: num  0.324 0.106 0.228 0.238 0.274 ...
 $ statistic: num  -2.46 4.14 -1.84 -4.57 -5.04 ...
 $ p.value  : num  1.38e-02 3.49e-05 6.53e-02 4.96e-06 4.71e-07 ...
library(JLutils)
Loading required package: plyr
str(tidy_detailed(reg))
'data.frame':   14 obs. of  10 variables:
 $ term          : chr  "etudNA" "etudSecondaire" "etudSupérieur" "etudTechnique/Professionnel" ...
 $ estimate      : num  2.149 0.948 1.89 1.047 -0.42 ...
 $ std.error     : num  0.33 0.197 0.195 0.19 0.228 ...
 $ statistic     : num  6.51 4.8 9.68 5.52 -1.84 ...
 $ p.value       : num  7.67e-11 1.56e-06 3.63e-22 3.43e-08 6.53e-02 ...
 $ variable      : chr  "etud" "etud" "etud" "etud" ...
 $ variable_label: chr  "etud" "etud" "etud" "etud" ...
 $ level         : chr  NA "Secondaire" "Supérieur" "Technique/Professionnel" ...
 $ level_detail  : chr  "NA vs. Primaire" "Secondaire vs. Primaire" "Supérieur vs. Primaire" "Technique/Professionnel vs. Primaire" ...
 $ label         : chr  "NA vs. Primaire" "Secondaire vs. Primaire" "Supérieur vs. Primaire" "Technique/Professionnel vs. Primaire" ...

Il est possible de combiner tidy_detailed avec ggcoef pour personnaliser un peu plus le résultat.

td <- tidy_detailed(reg, exponentiate = TRUE, conf.int = TRUE)
td$label <- factor(td$label, rev(td$label))  # Pour fixer l'ordre pour ggplot2
ggcoef(td, mapping = aes(y = level_detail, x = estimate, 
  colour = variable_label), exponentiate = TRUE)

L’extension effects propose une représentation graphique résumant les effets de chaque variable du modèle. Pour cela, il suffit d’appliquer la méthode plot au résultat de la fonction allEffects. Nous obtenons alors la figure ci-dessous.

library(effects)
plot(allEffects(reg))
Représentation graphique de l’effet de chaque variable du modèle logistique

Une manière de tester la qualité d’un modèle est le calcul d’une matrice de confusion, c’est-à-dire le tableau croisé des valeurs observées et celles des valeurs prédites en appliquant le modèle aux données d’origine.

La méthode predict avec l’argument type="response" permet d’appliquer notre modèle logistique à un tableau de données et renvoie pour chaque individu la probabilité qu’il ait vécu le phénomène étudié.

sport.pred <- predict(reg, type = "response", newdata = d)
head(sport.pred)
         1          2          3          4 
0.61251214 0.73426878 0.16004519 0.70335574 
         5          6 
0.07318189 0.34841147 

Or notre variable étudiée est de type binaire. Nous devons donc transformer nos probabilités prédites en une variable du type « oui / non ». Usuellement, les probabilités prédites seront réunies en deux groupes selon qu’elles soient supérieures ou inférieures à la moitié. La matrice de confusion est alors égale à :

table(sport.pred > 0.5, d$sport)
       
         Non  Oui
  FALSE 1074  384
  TRUE   199  336

Nous avons donc 583 (384+199) prédictions incorrectes sur un total de 1993, soit un taux de mauvais classement de 29,3 %.

Sélection de modèles

Il est toujours tentant lorsque l’on recherche les facteurs associés à un phénomène d’inclure un nombre important de variables explicatives potentielles dans un mmodèle logistique. Cependant, un tel modèle n’est pas forcément le plus efficace et certaines variables n’auront probablement pas d’effet significatif sur la variable d’intérêt.

La technique de sélection descendante pas à pas est une approche visant à améliorer son modèle explicatif2. On réalise un premier modèle avec toutes les variables spécifiées, puis on regarde s’il est possible d’améliorer le modèle en supprimant une des variables du modèle. Si plusieurs variables permettent d’améliorer le modèle, on supprimera la variable dont la suppression améliorera le plus le modèle. Puis on recommence le même procédé pour voir si la suppression d’une seconde variable peut encore améliorer le modèle et ainsi de suite. Lorsque le modèle ne peut plus être améliorer par la suppresion d’une variable, on s’arrête.

Il faut également définir un critère pour déterminer la qualité d’un modèle. L’un des plus utilisés est le Akaike Information Criterion ou AIC. Plus l’AIC sera faible, meilleure sera le modèle.

La fonction step permet justement de sélectionner le meilleur modèle par une procédure pas à pas descendante basée sur la minimisation de l’AIC. La fonction affiche à l’écran les différentes étapes de la sélection et renvoie le modèle final.

reg2 <- step(reg)
Start:  AIC=2235.94
sport ~ sexe + grpage + etud + relig + heures.tv

            Df Deviance    AIC
- relig      5   2210.2 2230.2
<none>           2205.9 2235.9
- heures.tv  1   2219.3 2247.3
- sexe       1   2223.2 2251.2
- grpage     3   2258.7 2282.7
- etud       4   2329.6 2351.6

Step:  AIC=2230.15
sport ~ sexe + grpage + etud + heures.tv

            Df Deviance    AIC
<none>           2210.2 2230.2
- heures.tv  1   2223.7 2241.7
- sexe       1   2226.1 2244.1
- grpage     3   2260.3 2274.3
- etud       4   2333.8 2345.8

Le modèle initial a un AIC de 2235,9. À la première étape, il apparait que la suppression de la variable religion permet diminuer l’AIC à 2230,2. Lors de la seconde étape, toute suppression d’une autre variable ferait augmenter l’AIC. La procédure s’arrête donc.

Régression logistique multinomiale

La régression logistique multinomiale est une extension de la régression logistique aux variables qualitatives à trois modalités ou plus. Dans ce cas de figure, chaque modalité de la variable d’intérêt sera comparée à la modalité de réference. Les odds ratio seront donc exprimés par rapport à cette dernière.

Nous allons prendre pour exemple la variable trav.satisf, à savoir la satisfaction ou l’insatisfaction au travail.

freq(d$trav.satisf)
                 n    % val%
Satisfaction   480 24.0 45.8
Insatisfaction 117  5.9 11.2
Equilibre      451 22.6 43.0
NA             952 47.6   NA

Nous allons choisir comme modalité de référence la position intermédiaire, à savoir l’« équilibre ».

d$trav.satisf <- relevel(d$trav.satisf, "Equilibre")

Enfin, nous allons aussi en profiter pour raccourcir les étiquettes de la variable trav.imp :

levels(d$trav.imp) <- c("Le plus", "Aussi", "Moins", 
  "Peu")

Pour calculer un modèle logistique multinomial, nous allons utiliser la fonction multinom de l’extension nnet3. La syntaxe de multinom est similaire à celle de glm, le paramètre family en moins.

library(nnet)
regm <- multinom(trav.satisf ~ sexe + etud + grpage + 
  trav.imp, data = d)
# weights:  39 (24 variable)
initial  value 1151.345679 
iter  10 value 977.348901
iter  20 value 969.849189
iter  30 value 969.522965
final  value 969.521855 
converged

Comme pour la régression logistique, il est possible de réaliser une sélection pas à pas descendante :

regm2 <- step(regm)
Start:  AIC=1987.04
trav.satisf ~ sexe + etud + grpage + trav.imp

trying - sexe 
# weights:  36 (22 variable)
initial  value 1151.345679 
iter  10 value 978.538886
iter  20 value 970.453555
iter  30 value 970.294459
final  value 970.293988 
converged
trying - etud 
# weights:  27 (16 variable)
initial  value 1151.345679 
iter  10 value 987.907714
iter  20 value 981.785467
iter  30 value 981.762800
final  value 981.762781 
converged
trying - grpage 
# weights:  30 (18 variable)
initial  value 1151.345679 
iter  10 value 979.485430
iter  20 value 973.175923
final  value 973.172389 
converged
trying - trav.imp 
# weights:  30 (18 variable)
initial  value 1151.345679 
iter  10 value 998.803976
iter  20 value 994.417973
iter  30 value 994.378914
final  value 994.378869 
converged
           Df      AIC
- grpage   18 1982.345
- sexe     22 1984.588
<none>     24 1987.044
- etud     16 1995.526
- trav.imp 18 2024.758
# weights:  30 (18 variable)
initial  value 1151.345679 
iter  10 value 979.485430
iter  20 value 973.175923
final  value 973.172389 
converged

Step:  AIC=1982.34
trav.satisf ~ sexe + etud + trav.imp

trying - sexe 
# weights:  27 (16 variable)
initial  value 1151.345679 
iter  10 value 976.669670
iter  20 value 973.928385
iter  20 value 973.928377
iter  20 value 973.928377
final  value 973.928377 
converged
trying - etud 
# weights:  18 (10 variable)
initial  value 1151.345679 
iter  10 value 988.413720
final  value 985.085797 
converged
trying - trav.imp 
# weights:  21 (12 variable)
initial  value 1151.345679 
iter  10 value 1001.517287
final  value 998.204280 
converged
           Df      AIC
- sexe     16 1979.857
<none>     18 1982.345
- etud     10 1990.172
- trav.imp 12 2020.409
# weights:  27 (16 variable)
initial  value 1151.345679 
iter  10 value 976.669670
iter  20 value 973.928385
iter  20 value 973.928377
iter  20 value 973.928377
final  value 973.928377 
converged

Step:  AIC=1979.86
trav.satisf ~ etud + trav.imp

trying - etud 
# weights:  15 (8 variable)
initial  value 1151.345679 
iter  10 value 986.124104
final  value 986.034023 
converged
trying - trav.imp 
# weights:  18 (10 variable)
initial  value 1151.345679 
iter  10 value 1000.225356
final  value 998.395273 
converged
           Df      AIC
<none>     16 1979.857
- etud      8 1988.068
- trav.imp 10 2016.791

La plupart des fonctions vues précédemment fonctionnent :

summary(regm2)
Call:
multinom(formula = trav.satisf ~ etud + trav.imp, data = d)

Coefficients:
               (Intercept) etudSecondaire
Satisfaction    -0.1110996     0.04916210
Insatisfaction  -1.1213760    -0.09737523
               etudTechnique/Professionnel
Satisfaction                    0.07793241
Insatisfaction                  0.08392603
               etudSupérieur      etudNA
Satisfaction      0.69950061 -0.53841577
Insatisfaction    0.07755307 -0.04364055
               trav.impAussi trav.impMoins
Satisfaction       0.2578973    -0.1756206
Insatisfaction    -0.2279774    -0.5330349
               trav.impPeu
Satisfaction    -0.5995051
Insatisfaction   1.3401509

Std. Errors:
               (Intercept) etudSecondaire
Satisfaction     0.4520902      0.2635573
Insatisfaction   0.6516992      0.3999875
               etudTechnique/Professionnel
Satisfaction                     0.2408483
Insatisfaction                   0.3579684
               etudSupérieur    etudNA
Satisfaction       0.2472571 0.5910993
Insatisfaction     0.3831110 0.8407592
               trav.impAussi trav.impMoins
Satisfaction       0.4260623     0.4115818
Insatisfaction     0.6213781     0.5941721
               trav.impPeu
Satisfaction     0.5580115
Insatisfaction   0.6587383

Residual Deviance: 1947.857 
AIC: 1979.857 
odds.ratio(regm2)
                                                 OR
Satisfaction/(Intercept)                   0.894850
Satisfaction/etudSecondaire                1.050391
Satisfaction/etudTechnique/Professionnel   1.081050
Satisfaction/etudSupérieur                 2.012747
Satisfaction/etudNA                        0.583672
Satisfaction/trav.impAussi                 1.294206
Satisfaction/trav.impMoins                 0.838936
Satisfaction/trav.impPeu                   0.549083
Insatisfaction/(Intercept)                 0.325831
Insatisfaction/etudSecondaire              0.907216
Insatisfaction/etudTechnique/Professionnel 1.087548
Insatisfaction/etudSupérieur               1.080640
Insatisfaction/etudNA                      0.957298
Insatisfaction/trav.impAussi               0.796142
Insatisfaction/trav.impMoins               0.586821
Insatisfaction/trav.impPeu                 3.819620
                                              2.5 %
Satisfaction/(Intercept)                   0.368918
Satisfaction/etudSecondaire                0.626629
Satisfaction/etudTechnique/Professionnel   0.674272
Satisfaction/etudSupérieur                 1.239720
Satisfaction/etudNA                        0.183242
Satisfaction/trav.impAussi                 0.561485
Satisfaction/trav.impMoins                 0.374447
Satisfaction/trav.impPeu                   0.183932
Insatisfaction/(Intercept)                 0.090838
Insatisfaction/etudSecondaire              0.414229
Insatisfaction/etudTechnique/Professionnel 0.539194
Insatisfaction/etudSupérieur               0.510007
Insatisfaction/etudNA                      0.184243
Insatisfaction/trav.impAussi               0.235544
Insatisfaction/trav.impMoins               0.183124
Insatisfaction/trav.impPeu                 1.050270
                                            97.5 %
Satisfaction/(Intercept)                    2.1706
Satisfaction/etudSecondaire                 1.7607
Satisfaction/etudTechnique/Professionnel    1.7332
Satisfaction/etudSupérieur                  3.2678
Satisfaction/etudNA                         1.8591
Satisfaction/trav.impAussi                  2.9831
Satisfaction/trav.impMoins                  1.8796
Satisfaction/trav.impPeu                    1.6391
Insatisfaction/(Intercept)                  1.1687
Insatisfaction/etudSecondaire               1.9869
Insatisfaction/etudTechnique/Professionnel  2.1936
Insatisfaction/etudSupérieur                2.2897
Insatisfaction/etudNA                       4.9740
Insatisfaction/trav.impAussi                2.6910
Insatisfaction/trav.impMoins                1.8805
Insatisfaction/trav.impPeu                 13.8912
                                                  p
Satisfaction/(Intercept)                   0.805878
Satisfaction/etudSecondaire                0.852027
Satisfaction/etudTechnique/Professionnel   0.746260
Satisfaction/etudSupérieur                 0.004669
Satisfaction/etudNA                        0.362363
Satisfaction/trav.impAussi                 0.544977
Satisfaction/trav.impMoins                 0.669600
Satisfaction/trav.impPeu                   0.282661
Insatisfaction/(Intercept)                 0.085306
Insatisfaction/etudSecondaire              0.807660
Insatisfaction/etudTechnique/Professionnel 0.814635
Insatisfaction/etudSupérieur               0.839581
Insatisfaction/etudNA                      0.958603
Insatisfaction/trav.impAussi               0.713701
Insatisfaction/trav.impMoins               0.369663
Insatisfaction/trav.impPeu                 0.041909
                                             
Satisfaction/(Intercept)                     
Satisfaction/etudSecondaire                  
Satisfaction/etudTechnique/Professionnel     
Satisfaction/etudSupérieur                 **
Satisfaction/etudNA                          
Satisfaction/trav.impAussi                   
Satisfaction/trav.impMoins                   
Satisfaction/trav.impPeu                     
Insatisfaction/(Intercept)                 . 
Insatisfaction/etudSecondaire                
Insatisfaction/etudTechnique/Professionnel   
Insatisfaction/etudSupérieur                 
Insatisfaction/etudNA                        
Insatisfaction/trav.impAussi                 
Insatisfaction/trav.impMoins                 
Insatisfaction/trav.impPeu                 * 
---
Signif. codes:  
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

De même, il est possible de calculer la matrice de confusion :

table(predict(regm2, newdata = d), d$trav.satisf)
                
                 Equilibre Satisfaction
  Equilibre            262          211
  Satisfaction         171          258
  Insatisfaction        18           11
                
                 Insatisfaction
  Equilibre                  49
  Satisfaction               45
  Insatisfaction             23

Données pondérées et l’extension survey

Lorsque l’on utilise des données pondérées, on aura recours à l’extension survey4.

Préparons des données d’exemple :

library(survey)
dw <- svydesign(ids = ~1, data = d, weights = ~poids)

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

odds.ratio(reg)
                                      OR   2.5 %
(Intercept)                      4.64553 2.57758
sexeHomme                        1.44089 1.12150
age                              0.95957 0.95181
religPratiquant occasionnel      1.05735 0.65902
religAppartenance sans pratique  1.17782 0.75911
religNi croyance ni appartenance 1.04068 0.64638
religRejet                       0.86190 0.42752
religNSP ou NVPR                 0.79707 0.32481
heures.tv                        0.83356 0.77019
                                 97.5 %         p
(Intercept)                      8.3726 3.522e-07
sexeHomme                        1.8513  0.004325
age                              0.9674 < 2.2e-16
religPratiquant occasionnel      1.6965  0.817180
religAppartenance sans pratique  1.8275  0.465321
religNi croyance ni appartenance 1.6755  0.869658
religRejet                       1.7376  0.677858
religNSP ou NVPR                 1.9560  0.620504
heures.tv                        0.9022 6.786e-06
                                    
(Intercept)                      ***
sexeHomme                        ** 
age                              ***
religPratiquant occasionnel         
religAppartenance sans pratique     
religNi croyance ni appartenance    
religRejet                          
religNSP ou NVPR                    
heures.tv                        ***
---
Signif. codes:  
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dans ses dernières versions, survey fournit une méthode AIC.svyglm permettant d’estimer un AIC sur un modèle calculé avec svyglm. Il est dès lors possible d’utiliser la fonction step pour réaliser une sélection descendante pas à pas.

L’extension effects n’est quant à elle pas compatible avec svyglm5.


  1. Pour plus de détails, voir http://www.spc.univ-lyon1.fr/polycop/odds%20ratio.htm.

  2. Il existe également des méthodes de sélection ascendante pas à pas, mais nous les aborderons pas ici.

  3. Une alternative est d’avoir recours à l’extension mlogit que nous n’aborderons pas ici. Voir http://www.ats.ucla.edu/stat/r/dae/mlogit.htm (en anglais) pour plus de détails.

  4. Voir le chapitre dédié aux données pondérées.

  5. Compatibilité qui pourra éventuellement être introduite dans une future version de l’extension.