Scarpe 2018 Donna 00015 BIANCO Primavera 110 Estate DEI COLLI SOUND Sneaker 5v1Tt5gwq

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 sur le fait de pratiquer un sport.

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.

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

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, 99), 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.2
[65,99] 380 19.0 19.0

Jetons maintenant un oeil à la variable nivetud :

freq(d$nivetud)

                                                                  n    % val%
N'a jamais fait d'etudes                                         39  2.0  2.1
A arrete ses etudes, avant la derniere annee d'etudes primaires  86  4.3  4.6
Derniere annee d'etudes primaires                               341 17.1 18.1
1er cycle                                                       204 10.2 10.8
2eme cycle                                                      183  9.2  9.7
Enseignement technique ou professionnel court                   463 23.2 24.5
Enseignement technique ou professionnel long                    131  6.6  6.9
Enseignement superieur y compris technique superieur            441 22.1 23.4
NA                                                              112  5.6   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 addNAstr fournie par questionr¹ :

levels(d$etud)

[1] "Primaire"                "Secondaire"             
[3] "Technique/Professionnel" "Supérieur"              

d$etud <- addNAstr(d$etud, "manquant")
levels(d$etud)

[1] "Primaire"                "Secondaire"             
[3] "Technique/Professionnel" "Supérieur"              
[5] "manquant"               

Régression logistique binaire

La fonction glm (pour COLLI Estate Primavera Scarpe BIANCO 00015 Sneaker Donna 110 SOUND 2018 DEI 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 Femme 243 Mustang 520 Baskets Hautes 1146 qv7ZU 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)                         sexeHomme  
                       -0.798368                          0.439694  
                   grpage[25,45)                     grpage[45,65)  
                       -0.420448                         -1.085434  
                   grpage[65,99]                    etudSecondaire  
                       -1.381353                          0.950571  
     etudTechnique/Professionnel                     etudSupérieur  
                        1.049253                          1.891667  
                          etudNA       religPratiquant occasionnel  
                        2.150428                         -0.021904  
 religAppartenance sans pratique  religNi croyance ni appartenance  
                       -0.006696                         -0.215389  
                      religRejet                  religNSP ou NVPR  
                       -0.383543                         -0.083789  
                       heures.tv  
                       -0.120911  

Degrees of Freedom: 1994 Total (i.e. Null);  1980 Residual
  (5 observations deleted due to missingness)
Null Deviance:      2609 
Residual Deviance: 2206     AIC: 2236

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.8784  -0.8865  -0.4808   1.0033   2.4222  

Coefficients:
                                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)                      -0.798368   0.323903  -2.465 0.013708 *  
sexeHomme                         0.439694   0.106062   4.146 3.39e-05 ***
grpage[25,45)                    -0.420448   0.228053  -1.844 0.065236 .  
grpage[45,65)                    -1.085434   0.237716  -4.566 4.97e-06 ***
grpage[65,99]                    -1.381353   0.273796  -5.045 4.53e-07 ***
etudSecondaire                    0.950571   0.197442   4.814 1.48e-06 ***
etudTechnique/Professionnel       1.049253   0.189804   5.528 3.24e-08 ***
etudSupérieur                     1.891667   0.195218   9.690  < 2e-16 ***
etudNA                            2.150428   0.330229   6.512 7.42e-11 ***
religPratiquant occasionnel      -0.021904   0.189199  -0.116 0.907833    
religAppartenance sans pratique  -0.006696   0.174737  -0.038 0.969434    
religNi croyance ni appartenance -0.215389   0.193080  -1.116 0.264617    
religRejet                       -0.383543   0.285905  -1.342 0.179756    
religNSP ou NVPR                 -0.083789   0.411028  -0.204 0.838470    
heures.tv                        -0.120911   0.033591  -3.599 0.000319 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2609.2  on 1994  degrees of freedom
Residual deviance: 2206.2  on 1980  degrees of freedom
  (5 observations deleted due to missingness)
AIC: 2236.2

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

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)                        sexeHomme 
                       0.4500628                        1.5522329 
                   grpage[25,45)                    grpage[45,65) 
                       0.6567525                        0.3377552 
                   grpage[65,99]                   etudSecondaire 
                       0.2512383                        2.5871865 
     etudTechnique/Professionnel                    etudSupérieur 
                       2.8555168                        6.6304155 
                          etudNA      religPratiquant occasionnel 
                       8.5885303                        0.9783342 
 religAppartenance sans pratique religNi croyance ni appartenance 
                       0.9933267                        0.8062276 
                      religRejet                 religNSP ou NVPR 
                       0.6814428                        0.9196256 
                       heures.tv 
                       0.8861129 

exp(confint(reg))

Waiting for profiling to be done...

                                     2.5 %     97.5 %
(Intercept)                      0.2377004  0.8473181
sexeHomme                        1.2614265  1.9120047
grpage[25,45)                    0.4194391  1.0274553
grpage[45,65)                    0.2115307  0.5380894
grpage[65,99]                    0.1463860  0.4288023
etudSecondaire                   1.7652682  3.8327896
etudTechnique/Professionnel      1.9804881  4.1729378
etudSupérieur                    4.5517938  9.7950691
etudNA                           4.5347799 16.5885323
religPratiquant occasionnel      0.6758807  1.4198530
religAppartenance sans pratique  0.7063000  1.4020242
religNi croyance ni appartenance 0.5524228  1.1782475
religRejet                       0.3868800  1.1889781
religNSP ou NVPR                 0.3996746  2.0215562
heures.tv                        0.8290223  0.9457756

On pourra faciliter la lecture en combinant les deux :

exp(cbind(coef(reg), confint(reg)))

Waiting for profiling to be done...

                                               2.5 %     97.5 %
(Intercept)                      0.4500628 0.2377004  0.8473181
sexeHomme                        1.5522329 1.2614265  1.9120047
grpage[25,45)                    0.6567525 0.4194391  1.0274553
grpage[45,65)                    0.3377552 0.2115307  0.5380894
grpage[65,99]                    0.2512383 0.1463860  0.4288023
etudSecondaire                   2.5871865 1.7652682  3.8327896
etudTechnique/Professionnel      2.8555168 1.9804881  4.1729378
etudSupérieur                    6.6304155 4.5517938  9.7950691
etudNA                           8.5885303 4.5347799 16.5885323
religPratiquant occasionnel      0.9783342 0.6758807  1.4198530
religAppartenance sans pratique  0.9933267 0.7063000  1.4020242
religNi croyance ni appartenance 0.8062276 0.5524228  1.1782475
religRejet                       0.6814428 0.3868800  1.1889781
religNSP ou NVPR                 0.9196256 0.3996746  2.0215562
heures.tv                        0.8861129 0.8290223  0.9457756

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 Donna SOUND Scarpe COLLI DEI BIANCO 00015 Estate 2018 Primavera 110 Sneaker 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 %  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 ***
etudNA                           8.58853 4.53478 16.5885 7.419e-11 ***
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    
heures.tv                        0.88611 0.82902  0.9458 0.0003188 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1141103BLANCO Victoria Victoria Basket 141103 Blanc Basket 141103 Victoria Blanc 141103BLANCO qBXxtFZt

library(broom)
tmp <-Estate Sneaker 110 00015 2018 Scarpe Donna DEI SOUND Primavera COLLI BIANCO 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.45 1.552 0.657 0.338 0.251 ...
 $ std.error: num  0.324 0.106 0.228 0.238 0.274 ...
 $ statistic: num  -2.46 4.15 -1.84 -4.57 -5.05 ...
 $ p.value  : num  1.37e-02 3.39e-05 6.52e-02 4.97e-06 4.53e-07 ...
 $ conf.low : num  0.238 1.261 0.419 0.212 0.146 ...
 $ conf.high: num  0.847 1.912 1.027 0.538 0.429 ...

library(ggplot2)
ggplot(tmp) + aes(x = estimate, y = term, xmin = conf.low, xmax = conf.high) + geom_vline(xintercept = 1) +COLLI BIANCO Sneaker Scarpe DEI SOUND 2018 00015 Estate Donna Primavera 110 
 geom_errorbarh() + geom_point() + scale_x_log10()

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

devtools::install_github("larmarange/JLutils")

str(tidy(reg))

'data.frame':   15 obs. of  5 variables:
 $ term     : chr  "(Intercept)" "sexeHomme" "grpage[25,45)" "grpage[45,65)" ...
 $ estimate : num  -0.798 0.44 -0.42 -1.085 -1.381 ...
 $ std.error: num  0.324 0.106 0.228 0.238 0.274 ...
 $ statistic: num  -2.46 4.15 -1.84 -4.57 -5.05 ...
 $ p.value  : num  1.37e-02 3.39e-05 6.52e-02 4.97e-06 4.53e-07 ...

library(JLutils)
str(tidy_detailed(reg))

'data.frame':   15 obs. of  10 variables:
 $ term          : chr  "(Intercept)" "etudNA" "etudSecondaire" "etudSupérieur" ...
 $ estimate      : num  -0.798 2.15 0.951 1.892 1.049 ...
 $ std.error     : num  0.324 0.33 0.197 0.195 0.19 ...
 $ statistic     : num  -2.46 6.51 4.81 9.69 5.53 ...
 $ p.value       : num  1.37e-02 7.42e-11 1.48e-06 3.32e-22 3.24e-08 ...
 $ variable      : chr  NA "etud" "etud" "etud" ...
 $ variable_label: chr  NA "etud" "etud" "etud" ...
 $ level         : chr  NA NA "Secondaire" "Supérieur" ...
 $ level_detail  : chr  NA "NA vs. Primaire" "Secondaire vs. Primaire" "Supérieur vs. Primaire" ...
 $ label         : chr  NA "NA vs. Primaire" "Secondaire vs. Primaire" "Supérieur vs. Primaire" ...

td <- tidy_detailed(reg, exponentiate = TRUE, conf.int = TRUE)
td$level_detail <- factor(td$label, rev(td$level_detail))  # Pour fixer l'ordre pour ggplot2
ggcoef(td, mapping = aes(COLLI SOUND Primavera Scarpe Donna Estate 110 Sneaker BIANCO 00015 2018 DEI y = level_detail, x = estimate, colour = variable_label), 
  exponentiate = TRUE)

library(effects)
plot(allEffects(reg))

sport.pred <- predict(reg, type = "response", newdata = d)
head(sport.pred)

         1          2          3          4          5          6 
0.61241192 0.73414575 0.15982958 0.70350157 0.07293505 0.34824228 

table(sport.pred > 110 Primavera 2018 Estate Sneaker BIANCO Donna SOUND COLLI Scarpe 00015 DEI 0.5, d$sport)

       
         Non  Oui
  FALSE 1076  384
  TRUE   199  336

drop1(reg, test = "Chisq")

Single term deletions

Model:
sport ~ sexe + grpage + etud + relig + heures.tv
          Df Deviance    AIC     LRT  Pr(>Chi)    
         2206.2 2236.2                      
sexe       1   2223.5 2251.5  17.309 3.176e-05 ***
grpage     3   2259.0 2283.0  52.803 2.020e-11 ***
etud       4   2330.0 2352.0 123.826 < 2.2e-16 ***
relig      5   2210.4 2230.4   4.232 0.5165401    
heures.tv  1   2219.6 2247.6  13.438 0.0002465 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

reg2 <- step(reg)

Start:  AIC=2236.17
sport ~ sexe + grpage + etud + relig + heures.tv

            Df Deviance    AIC
- relig      5   2210.4 2230.4
           2206.2 2236.2
- heures.tv  1   2219.6 2247.6
- sexe       1   2223.5 2251.5
- grpage     3   2259.0 2283.0
- etud       4   2330.0 2352.0

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

            Df Deviance    AIC
           2210.4 2230.4
- heures.tv  1   2224.0 2242.0
- sexe       1   2226.4 2244.4
- grpage     3   2260.6 2274.6
- etud       4   2334.3 2346.3

AIC(reg)

[1] 2236.173

AIC(reg2)

[1] 2230.404

anova(reg, reg2, test = "Chisq")

Analysis of Deviance Table

Model 1: sport ~ sexe + grpage + etud + relig + heures.tv
Model 2: sport ~ sexe + grpage + etud + heures.tv
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1      1980     2206.2                     
2      1985     2210.4 -5  -4.2319   0.5165

library(MASS)
reg2bis <- stepAIC(reg)

Start:  AIC=2236.17
sport ~ sexe + grpage + etud + relig + heures.tv

            Df Deviance    AIC
- relig      5   2210.4 2230.4
           2206.2 2236.2
- heures.tv  1   2219.6 2247.6
- sexe       1   2223.5 2251.5
- grpage     3   2259.0 2283.0
- etud       4   2330.0 2352.0

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

            Df Deviance    AIC
           2210.4 2230.4
- heures.tv  1   2224.0 2242.0
- sexe       1   2226.4 2244.4
- grpage     3   2260.6 2274.6
- etud       4   2334.3 2346.3

devtools::install_github("ewenharrison/finalfit")

dep <-COLLI DEI 2018 Primavera 00015 Scarpe Sneaker BIANCO 110 SOUND Donna Estate "sport"
vars <- c("sexe", "grpage", "etud", "relig", "heures.tv")

library(finalfit)
tab <- summary_factorlist(d, dep, vars, p = TRUE, add_dependent_label = TRUE)
tab

library(labelled)
var_label(d$sport) <- "Pratique du sport ?"
var_label(d$sexe) <- "Sexe"
var_labelMesh Baskets Coq Metallic Feminine Dynacomf Le W Femme Sportif wSq0ZX(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"
tab <- summary_factorlist(d, dep, vars, p = TRUE, add_dependent_label = TRUE)
knitr::kable(tab, row.names = FALSE)

tab <- finalfit(d, dep, vars)
knitr::kable(tab, row.names = FALSE)

vars_multi <- c("sexe", "grpage", "etud", "heures.tv")
tab <- finalfit(d, dep, vars, explanatory_multi = vars_multi)
knitr::kable(tab, row.names = FALSE)

var_label(d$heures.tv) <- "Nombre d'heures\nde TV par jour"Mixte Pacer Basses Adulte Puma Sneakers Next AInrW6Wd
or_plot(d, dep, vars_multi)

Waiting for profiling to be done...
Waiting for profiling to be done...

Warning: Removed 3 rows containing missing values (geom_errorbarh).

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)

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 nnet⁵. 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
     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
     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
     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 etudTechnique/Professionnel
Satisfaction    -0.1110996     0.04916210                  0.07793241
Insatisfaction  -1.1213760    -0.09737523                  0.08392603
               etudSupérieur etudmanquant trav.impAussi trav.impMoins
Satisfaction      0.69950061  -0.53841577     0.2578973    -0.1756206
Insatisfaction    0.07755307  -0.04364055    -0.2279774    -0.5330349
               trav.impPeu
Satisfaction    -0.5995051
Insatisfaction   1.3401509

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

Residual Deviance: 1947.857 
AIC: 1979.857 

odds.ratio(regm2)

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

table(predict(regm2, newdata = d), d$trav.satisf)

                
                 Equilibre Satisfaction Insatisfaction
  Equilibre            262          211             49
  Satisfaction         171          258             45
  Insatisfaction        18           11             23

La fonction Stan Smith adidas Baskets Femme Blanc d8gdqwBxz peut s’appliquer au résultat de multinom :

library(broom)
tidy(regm2, exponentiate = TRUE, conf.int = TRUE)

On notera la présence d’une colonne supplémentaire, y.level. De fait, la fonction Sneaker 00015 Primavera COLLI 2018 110 Estate Donna Scarpe BIANCO SOUND DEI ggcoef ne peut s’appliquer directement, car les coefficients vont se supperposer.

ggcoef(regm2, exponentiate = TRUE)

On a deux solutions possibles. Pour la première, la plus simple, il suffit d’ajouter des facettes avec facet_grid.

ggcoef(regm2, exponentiate = TRUE) + facet_grid(~y.level)

Pour la seconde, on va réaliser un graphique personnalisé, sur la même logique que Chaussures Ban Irregular Femmes Choice Joe AxqU7, décalant les points grâce à position_dodge.

ggplot(tidy(regm2, exponentiate = T, conf.int = TRUE)) + aes(x = term, y = estimate, 
  ymin = conf.low, ymax = conf.high, color = y.level) + geom_hline(yintercept = 1, 
  color = "gray25", linetype = "dotted") + geom_errorbar(position = position_dodge(0.5), 
  width = 0) + geom_point(position = position_dodge(width = 0.5)) + scale_y_log10() + 
 coord_flip()

Régression logistique ordinale

La régression logistique ordinale s’applique lorsque la variable à expliquer possède trois ou plus modalités qui sont ordonnées (par exemple : modéré, moyen, fort).

L’extension la plus utilisée pour réaliser des modèles ordinaux est ordinal et sa fonction clm. Il est même possible de réaliser des modèles ordinaux avec des effets aléatoires (modèles mixtes) à l’aide de la fonction clmm.

Pour une bonne introduction à l’extension ordinal, on pourra se référer au tutoriel officiel (en anglais) : https://cran.r-project.org/web/packages/ordinal/vignettes/clm_tutorial.pdf.

Une autre introduction pertinente (en français) et utilisant cette fois-ci l’extention VGAM et sa fonction DEI 00015 Scarpe 110 Primavera Donna SOUND 2018 COLLI Sneaker Estate BIANCO vglm est disponible sur le site de l’université de Lyon : https://eric.univ-lyon2.fr/~ricco/cours/didacticiels/data-mining/didacticiel_Reg_Logistique_Polytomique_Ordinale.pdf.

On va reprendre l’exemple précédent puisque la variable trav.satisf est une variable ordonnée.

freq(d$trav.satisf)

ATTENTION : Dans le cas d’une régression logistique ordinale, il importante que les niveaux du facteur soient classés selon leur ordre hiéarchique (du plus faible au plus fort). On va dès lors recoder notre variable à expliquer.

d$trav.satisf <- factor(d$trav.satisf, c("Insatisfaction", "Equilibre", "Satisfaction"), 
  ordered = TRUE)
freq(d$trav.satisf)

library(ordinal)
rego <- clm(trav.satisf ~ sexe + etud + trav.imp, data = d)
summary(rego)

formula: trav.satisf ~ sexe + etud + trav.imp
data:    d

 link  threshold nobs logLik  AIC     niter max.grad cond.H 
 logit flexible  1048 -978.61 1977.23 5(0)  5.41e-09 3.9e+02

Coefficients:
                            Estimate Std. Error z value Pr(>|z|)    
sexeHomme                   -0.16141    0.12215  -1.321 0.186369    
etudSecondaire               0.05558    0.23220   0.239 0.810819    
etudTechnique/Professionnel  0.03373    0.21210   0.159 0.873660    
etudSupérieur                0.61336    0.21972   2.792 0.005245 ** 
etudmanquant                -0.45555    0.48761  -0.934 0.350174    
trav.impAussi                0.35104    0.38578   0.910 0.362857    
trav.impMoins                0.02616    0.37174   0.070 0.943896    
trav.impPeu                 -1.66062    0.46204  -3.594 0.000325 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Threshold coefficients:
                         Estimate Std. Error z value
Insatisfaction|Equilibre  -2.0226     0.4189  -4.829
Equilibre|Satisfaction     0.3391     0.4113   0.825
(952 observations deleted due to missingness)

Une fois encore, il est possible de faire une sélection descendane pas à pas.

rego2 <- step(rego)

Start:  AIC=1977.23
trav.satisf ~ sexe + etud + trav.imp

           Df    AIC
- sexe      1 1977.0
        1977.2
- etud      4 1990.6
- trav.imp  3 2013.2

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

           Df    AIC
        1977.0
- etud      4 1990.6
- trav.imp  3 2011.6

L’extension broom ne propose pas de méthode Stan Smith adidas Baskets Femme Blanc d8gdqwBxz pour les objets clm. Cependant, on pourra en trouver une dans l’extension Copenhagen Grey Raven Women's 2 Women's FG 0 Sneakers ARKK dqwz6d. Pour rappel, cette extension est disponible uniquement sur GitHub et s’installe donc ainsi :

install_github("larmarange/JLutils")

library(JLutils)
tidy(rego2, exponentiate = TRUE, conf.int = TRUE)

La méthode tidy.clm étant disponible, on peut utiliser Chaussures Ban Irregular Femmes Choice Joe AxqU7.

library(JLutils)
library(GGally)
ggcoef(rego2, exponentiate = TRUE)

Warning: Removed 2 rows containing missing values (geom_errorbarh).

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

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

Préparons des données d’exemple :

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

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)                         sexeHomme  
                         1.53590                           0.36526  
                             age       religPratiquant occasionnel  
                        -0.04127                           0.05577  
 religAppartenance sans pratique  religNi croyance ni appartenance  
                         0.16367                           0.03988  
                      religRejet                  religNSP ou NVPR  
                        -0.14862                          -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

odds.ratio(reg)

library(broom)
tidy(reg, exponentiate = TRUE, conf.int = TRUE)

devtools::install_github("carlganz/svrepmisc")

library(survey)
dw <- svydesign(ids = ~Atwood Low Baskets Basses Vans Femme SaxwSq1, data = d, weights = ~poids)
dwr <- as.svrepdesign(dw, type = "bootstrap", replicates = 100)

library(svrepmisc)
regm <- svymultinom(trav.satisf ~ sexe + etud + grpage + trav.imp, design = dwr)

regm

                                         Coefficient       SE t value Pr(>|t|)
Equilibre.(Intercept)                        0.93947  1.79069  0.5246   0.6014
Satisfaction.(Intercept)                     0.74127  1.83550  0.4039   0.6875
Equilibre.sexeHomme                         -0.24304  0.29288 -0.8298   0.4092
Satisfaction.sexeHomme                      -0.27381  0.29696 -0.9220   0.3594
Equilibre.etudSecondaire                    -0.41802  0.60514 -0.6908   0.4918
Satisfaction.etudSecondaire                 -0.26579  0.55659 -0.4775   0.6343
Equilibre.etudTechnique/Professionnel       -0.47708  0.54825 -0.8702   0.3869
Satisfaction.etudTechnique/Professionnel    -0.23585  0.51393 -0.4589   0.6476
Equilibre.etudSupérieur                     -0.57579  0.57258 -1.0056   0.3178
Satisfaction.etudSupérieur                   0.36801  0.54181  0.6792   0.4991
Equilibre.etudmanquant                      -0.32346  5.59568 -0.0578   0.9541
Satisfaction.etudmanquant                   -0.65343  5.65957 -0.1155   0.9084
Equilibre.grpage[25,45)                      0.60710  0.63340  0.9585   0.3409
Satisfaction.grpage[25,45)                   0.46941  0.71325  0.6581   0.5124
Equilibre.grpage[45,65)                      0.57989  0.68468  0.8469   0.3997
Satisfaction.grpage[45,65)                   0.52873  0.74963  0.7053   0.4828
Equilibre.grpage[65,99]                     -3.20673 24.42192 -0.1313   0.8959
Satisfaction.grpage[65,99]                  11.60092 25.53337  0.4543   0.6509
Equilibre.trav.impAussi                      0.42551  1.66328  0.2558   0.7988
Satisfaction.trav.impAussi                   0.54156  1.68645  0.3211   0.7490
Equilibre.trav.impMoins                      0.73727  1.63309  0.4515   0.6529
Satisfaction.trav.impMoins                   0.66621  1.66346  0.4005   0.6899
Equilibre.trav.impPeu                       -1.54722  1.72341 -0.8978   0.3721
Satisfaction.trav.impPeu                    -1.47191  1.77820 -0.8278   0.4104

confint(regm)

                                               2.5 %     97.5 %
Equilibre.(Intercept)                     -2.6270085  4.5059460
Satisfaction.(Intercept)                  -2.9144426  4.3969816
Equilibre.sexeHomme                       -0.8263604  0.3402774
Satisfaction.sexeHomme                    -0.8652610  0.3176483
Equilibre.etudSecondaire                  -1.6232612  0.7872165
Satisfaction.etudSecondaire               -1.3743479  0.8427585
Equilibre.etudTechnique/Professionnel     -1.5690225  0.6148616
Satisfaction.etudTechnique/Professionnel  -1.2594308  0.7877238
Equilibre.etudSupérieur                   -1.7161773  0.5646004
Satisfaction.etudSupérieur                -0.7111056  1.4471296
Equilibre.etudmanquant                   -11.4682300 10.8213108
Satisfaction.etudmanquant                -11.9254421 10.6185783
Equilibre.grpage[25,45)                   -0.6544188  1.8686207
Satisfaction.grpage[25,45)                -0.9511533  1.8899648
Equilibre.grpage[45,65)                   -0.7837766  1.9435542
Satisfaction.grpage[45,65)                -0.9642809  2.0217400
Equilibre.grpage[65,99]                  -51.8472111 45.4337414
Satisfaction.grpage[65,99]               -39.2532028 62.4550416
Equilibre.trav.impAussi                   -2.8871952  3.7382160
Satisfaction.trav.impAussi                -2.8172965  3.9004214
Equilibre.trav.impMoins                   -2.5153112  3.9898428
Satisfaction.trav.impMoins                -2.6468716  3.9792818
Equilibre.trav.impPeu                     -4.9796878  1.8852516
Satisfaction.trav.impPeu                  -5.0135074  2.0696817

library(broom)
tidy(regm, exponentiate = TRUE, conf.int = TRUE)

rego <- svyolr(trav.satisf ~ sexe + etud + trav.imp, design = dw)
summary(rego)

Call:
svyolr(trav.satisf ~ sexe + etud + trav.imp, design = dw)

Coefficients:
                                    Value Std. Error      t value
sexeHomme                   -0.1328026100  0.1545640 -0.859207986
etudSecondaire              -0.0427151234  0.2751085 -0.155266460
etudTechnique/Professionnel  0.0004261287  0.2492533  0.001709621
etudSupérieur                0.6424698737  0.2593446  2.477282672
etudmanquant                -0.4694599623  0.5033490 -0.932672851
trav.impAussi                0.1207125976  0.5044264  0.239306656
trav.impMoins                0.0526740003  0.4895782  0.107590577
trav.impPeu                 -1.5303889517  0.7215699 -2.120915790

Intercepts:
                         Value   Std. Error t value
Insatisfaction|Equilibre -2.0392  0.5427    -3.7577
Equilibre|Satisfaction    0.2727  0.5306     0.5140
(952 observations deleted due to missingness)

confint(rego)

                                 2.5 %     97.5 %
sexeHomme                   -0.4357425  0.1701372
etudSecondaire              -0.5819179  0.4964876
etudTechnique/Professionnel -0.4881013  0.4889536
etudSupérieur                0.1341638  1.1507759
etudmanquant                -1.4560059  0.5170860
trav.impAussi               -0.8679450  1.1093702
trav.impMoins               -0.9068816  1.0122296
trav.impPeu                 -2.9446399 -0.1161380
Insatisfaction|Equilibre    -3.1027595 -0.9755811
Equilibre|Satisfaction       0.1607965  0.3846345

library(JLutils)
library(broom)
tidy(rego, exponentiate = TRUE, conf.int = TRUE)

dwr <- as.svrepdesign(dw, type = "bootstrap", replicates = 100)
library(svrepmisc)
rego_alt <- svyclm(trav.satisf ~ sexe + etud + trav.imp, COLLI Sneaker BIANCO 2018 DEI Primavera Scarpe 00015 Donna SOUND 110 Estate design = dwr)

Warning: (-1) Model failed to converge with max|grad| = 1.04205e-06 (tol = 1e-06) 
In addition: step factor reduced below minimum

Warning: (-1) Model failed to converge with max|grad| = 3.22346e-05 (tol = 1e-06) 
In addition: step factor reduced below minimum

rego_alt

                            Coefficient          SE t value  Pr(>|t|)    
Insatisfaction|Equilibre    -2.03917466  0.58882413 -3.4631 0.0008193 ***
Equilibre|Satisfaction       0.27271800  0.56130300  0.4859 0.6282431    
sexeHomme                   -0.13280477  0.17795218 -0.7463 0.4574344    
etudSecondaire              -0.04271403  0.28310986 -0.1509 0.8804124    
etudTechnique/Professionnel  0.00042809  0.27715652  0.0015 0.9987710    
etudSupérieur                0.64247607  0.27536180  2.3332 0.0218659 *  
etudmanquant                -0.46945975  0.54203306 -0.8661 0.3887334    
trav.impAussi                0.12071087  0.52291049  0.2308 0.8179598    
trav.impMoins                0.05267146  0.52268594  0.1008 0.9199566    
trav.impPeu                 -1.53039056  0.72723032 -2.1044 0.0381308 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

confint(rego_alt)

                                  2.5 %      97.5 %
Insatisfaction|Equilibre    -3.20897657 -0.86937274
Equilibre|Satisfaction      -0.84240838  1.38784439
sexeHomme                   -0.48633783  0.22072829
etudSecondaire              -0.60516118  0.51973312
etudTechnique/Professionnel -0.55019173  0.55104790
etudSupérieur                0.09542179  1.18953035
etudmanquant                -1.54630303  0.60738352
trav.impAussi               -0.91814209  1.15956382
trav.impMoins               -0.98573540  1.09107832
trav.impPeu                 -2.97516053 -0.08562059

tidy(rego_alt, exponentiate = TRUE, conf.int = TRUE)