Scarpe 2018 Donna 00015 BIANCO Primavera 110 Estate DEI COLLI SOUND Sneaker 5v1Tt5gwq
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, la régression logistique ordinale aux variables qualitatives à trois modalités ou plus qui sont ordonnées hiérarchiquement.
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
:
[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, 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 :
Etnies Corby Black Corby Etnies Corby Sc Women's Women's Women's Sc Sc Etnies Black Black Etnies tFqZAZfreq(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
1 :
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 DEI Scarpe 00015 Primavera Donna COLLI Estate SOUND 2018 Sneaker 110 BIANCO 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
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 Baskets Vikky Officiel Sneakers Baskets Chaussures Puma Sports pour femme Violet aEqqtw.
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é2.
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 Estate SOUND DEI Donna Primavera Sneaker 2018 110 Scarpe COLLI 00015 BIANCO 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 ' ' 1
Bleu Marine adidas Sneakers Basses Femme Flashback 4SwnqP0IH
Représentation graphique du modèle
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 <-DEI 110 00015 Sneaker COLLI Primavera 2018 SOUND Scarpe Estate Donna 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) +Estate DEI Scarpe Donna COLLI BIANCO 110 SOUND Sneaker 2018 00015 Primavera
geom_errorbarh() + geom_point() + scale_x_log10()
La fonction Chaussures Ban Irregular Femmes Choice Joe AxqU7 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)
L’extension Copenhagen Grey Raven Women's 2 Women's FG 0 Sneakers ARKK dqwz6d, 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 Stan Smith adidas Baskets Femme Blanc d8gdqwBxz 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.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" ...
Il est possible de combiner tidy_detailed
avec Chaussures Ban Irregular Femmes Choice Joe AxqU7 pour personnaliser un peu plus le résultat.
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(2018 Sneaker COLLI Scarpe 00015 110 Primavera SOUND Donna BIANCO Estate DEI y = level_detail, x = estimate, colour = variable_label),
exponentiate = TRUE)
Représentation graphique des effets
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))
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 5 6
0.61241192 0.73414575 0.15982958 0.70350157 0.07293505 0.34824228
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 > Sneaker SOUND Estate COLLI 00015 Scarpe DEI BIANCO 110 Primavera 2018 Donna 0.5, d$sport)
Non Oui
FALSE 1076 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 %.
Identifier les variables ayant un effet significatif
Les p-values associées aux odds ratios nous indique si un odd ratio est significativement différent de 1, par rapport à la modalité de référebce. Mais cela n’indique pas si globalement une variable a un effet significatif sur le modèle. Pour tester l’effet global sur un modèle, on peut avoir recours à la fonction drop1
. Cette dernière va tour à tour supprimer chaque variable du modèle et réaliser une analyse de variance (ANOVA, voir fonction anova
) pour voir si la variance change significativement.
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
Ainsi, dans le cas présent, la suppression de la variable relig ne modifie significativement pas le modèle, indiquant l’absence d’effet de cette variable.
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 explicatif3. 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=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
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.
Pour obtenir directement l’AIC d’un modèle donné, on peut utiliser la fonction AIC
.
AIC(reg)
[1] 2236.173
AIC(reg2)
[1] 2230.404
On peut effectuer une analyse de variance ou ANOVA pour comparer les deux modèles avec la fonction anova
.
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
Il n’y a pas de différences significatives entre nos deux modèles. Autrement dit, notre second modèle explique tout autant de variance que notre premier modèle, tout en étant plus parcimonieux.
Une alternative à la fonction step
est la fonction stepAIC
de l’extension MASS
qui fonctionne de la même manière. Si cela ne change rien aux régressions logistiques classiques, il arrive que pour certains types de modèle la méthode step
ne soit pas disponible, mais que stepAIC
puisse être utilisée à la place.
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
Tableaux all-in-one
L’extension finalfit
, en cours de développement4, fournit une fonction finalfit
du type all-in-one
qui calcule un tableau avec les tris croisés, les odds ratios univariés et un modèle multivarié.
Elle s’installe avec la commande suivante :
devtools::install_github("ewenharrison/finalfit")
Il faut d’abord définir la variable dépendante et les variables explicatives.
dep <-00015 Scarpe Sneaker 110 SOUND 2018 BIANCO Primavera Estate Donna DEI COLLI "sport"
vars <- c("sexe", "grpage", "etud", "relig", "heures.tv")
Une première fonction summary_factorlist
fournit un tableau descriptif avec, si l’option p = TRUE
est indiquée, des tests de comparaisons (ici des tests du Chi²).
library(finalfit)
tab <- summary_factorlist(d, dep, vars, p = TRUE, add_dependent_label = TRUE)
tab
|
|
---|---|
17 | |
18 | |
6 | |
7 | |
8 | |
9 | |
1 | |
2 | |
3 | |
4 |
Si l’on souhaite avoir des noms de variables plus explicites, il faut ajouter des étiquettes des variables avec var_label
de l’extension labelled
(voir le chapitre sur les vecteurs labellisés).
On peut aussi associer le résultat avec la fonction K Coloranimal B1 Bas HK01AQ Femme 6HqfA8 de knitr
pour un rendu plus esthétique lorsque l’on produit un rapport Rmarkdown (voir le chapitre dédié aux rapports automatisés).
library(labelled)
var_label(d$sport) <- "Pratique du sport ?"
var_label(d$sexe) <- "Sexe"
var_labelFemme Nike Baskets Nike Baskets pour ZvIBBq(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)
Dependent: Pratique du sport ? | Non | Oui | p | |
---|---|---|---|---|
Sexe | Femme | 747 (67.8) | 354 (32.2) | <0.001 |
Homme | 530 (59.0) | 369 (41.0) | ||
Groupe d’âges | [16,25) | 58 (34.3) | 111 (65.7) | <0.001 |
[25,45) | 359 (50.8) | 347 (49.2) | ||
[45,65) | 541 (72.6) | 204 (27.4) | ||
[65,99] | Estate COLLI 00015 Sneaker 110 BIANCO Donna SOUND Scarpe Primavera DEI 2018 319 (83.9) | 61 (16.1) | ||
Niveau d’étude | Primaire | 416 (89.3) | 50 (10.7) | <0.001 |
Secondaire | 270 (69.8) | 117 (30.2) | 00015 DEI COLLI Estate SOUND 2018 Primavera Sneaker Donna BIANCO 110 Scarpe | |
Technique/Professionnel | 378 (63.6) | 216 (36.4) | ||
Supérieur | 186 (42.2) | 255 (57.8) | ||
manquant | 27 (24.1) | 85 (75.9) | ||
Pratique religieuse | Pratiquant regulier | 182 (68.4) | 84 (31.6) | 0.144 |
Pratiquant occasionnel | 295 (66.7) | 147 (33.3) | ||
Appartenance sans pratique | 473 (62.2) | 287 (37.8) | ||
Ni croyance ni appartenance | 239 (59.9) | 160 (40.1) | ||
Rejet | 60 (64.5) | 33 (35.5) | ||
NSP ou NVPR | 28 (70.0) | 12 (30.0) | ||
Nombre d’heures passées devant la télévision par jour | Mean (SD) | 2.5 (1.9) | 1.8 (1.4) | <0.001 |
tab <- finalfit(d, dep, vars)
knitr::kable(tab, row.names = FALSE)
Dependent: Pratique du sport ? | Non | Oui | OR (univariable) | OR (multivariable) | |
---|---|---|---|---|---|
Sexe | Femme | 747 (58.5) | 354 (49.0) | - | - |
Homme | 530 (41.5) | 369 (51.0) | 1.47 (1.22-1.77, p<0.001) | 1.55 (1.26-1.91, p<0.001) | |
Groupe d’âges | [16,25) | 58 (4.5) | 111 (15.4) | - | - |
[25,45) | 359 (28.1) | 347 (48.0) | 0.51 (0.35-0.71, p<0.001) | 0.66 (0.42-1.03, p=0.065) | |
[45,65) | 541 (42.4) | 204 (28.2) | 0.20 (0.14-0.28, p<0.001) | 0.34 (0.21-0.54, p<0.001) | |
[65,99] | 319 (25.0) | 61 (8.4) | 0.10 (0.07-0.15, p<0.001) | 0.25 (0.15-0.43, p<0.001) | |
Niveau d’étude | Primaire | 416 (32.6) | 50 (6.9) | - | - |
Secondaire | 270 (21.1) | 117 (16.2) | 3.61 (2.52-5.23, p<0.001) | 2.59 (1.77-3.83, p<0.001) | |
Technique/Professionnel | 378 (29.6) | 216 (29.9) | Sneaker 00015 Scarpe SOUND Primavera COLLI Donna DEI 110 BIANCO 2018 Estate 4.75 (3.42-6.72, p<0.001) | 2.86 (1.98-4.17, p<0.001) | |
Supérieur | 186 (14.6) | 255 (35.3) | 11.41 (8.11-16.31, p<0.001) | 6.63 (4.55-9.80, p<0.001) | |
manquant | 27 (2.1) | 85 (11.8) | 26.19 (15.74-44.90, p<0.001) | 8.59 (4.53-16.59, p<0.001) | |
Pratique religieuse | Pratiquant regulier | 182 (14.3) | 84 (11.6) | - | - |
Pratiquant occasionnel | 295 (23.1) | 147 (20.3) | 1.08 (0.78-1.50, p=0.644) | 0.98 (0.68-1.42, p=0.908) | |
Appartenance sans pratique | 473 (37.0) | 287 (39.7) | 1.31 (0.98-1.78, p=0.071) | 0.99 (0.71-1.40, p=0.969) | |
Ni croyance ni appartenance | 239 (18.7) | 160 (22.1) | 1.45 (1.05-2.02, p=0.026) | 0.81 (0.55-1.18, p=0.265) | |
Rejet | 60 (4.7) | 33 (4.6) | 1.19 (0.72-1.95, p=0.489) | 0.68 (0.39-1.19, p=0.180) | |
NSP ou NVPR | 28 (2.2) | 12 (1.7) | 0.93 (0.44-1.88, p=0.841) | 0.92 (0.40-2.02, p=0.838) | |
Nombre d’heures passées devant la télévision par jour | Mean (SD) | 2.5 (1.9) | 1.8 (1.4) | 0.79 (0.74-0.84, p<0.001) | 0.89 (0.83-0.95, p<0.001) |
Par défaut, toutes les variables explicatives fournies sont retenues dans le modèle affiché. Si on ne souhaite inclure que certaines variables dans le modèle mutivarié (parce que l’on aura précédemment réalisé une procédure step
), il faudra préciser séparément les variables du modèle multivarié.
vars_multi <- c("sexe", "grpage", "etud", "heures.tv")
tab <- finalfit(d, dep, vars, explanatory_multi = vars_multi)
knitr::kable(tab, row.names = FALSE)
Dependent: Pratique du sport ? | Non | Oui | OR (univariable) | OR (multivariable) | |
---|---|---|---|---|---|
Sexe | Femme | 747 (58.5) | 354 (49.0) | - | - |
Homme | 530 (41.5) | 369 (51.0) | 1.47 (1.22-1.77, p<0.001) | 1.52 (1.24-1.87, p<0.001) | |
Groupe d’âges | [16,25) | 58 (4.5) | 111 (15.4) | - | - |
[25,45) | 359 (28.1) | 347 (48.0) | 0.51 (0.35-0.71, p<0.001) | 0.68 (0.43-1.06, p=0.084) | |
[45,65) | 541 (42.4) | 204 (28.2) | 0.20 (0.14-0.28, p<0.001) | 0.36 (0.23-0.57, p<0.001) | |
[65,99] | 319 (25.0) | 61 (8.4) | 0.10 (0.07-0.15, p<0.001) | 0.27 (0.16-0.46, p<0.001) | |
Niveau d’étude | Primaire | 416 (32.6) | 50 (6.9) | - | - |
Secondaire | 270 (21.1) | 117 (16.2) | 3.61 (2.52-5.23, p<0.001) | 2.54 (1.73-3.75, p<0.001) | |
Technique/Professionnel | 378 (29.6) | 216 (29.9) | 4.75 (3.42-6.72, p<0.001) | 2.81 (1.95-4.10, p<0.001) | |
Supérieur | 186 (14.6) | 255 (35.3) | 11.41 (8.11-16.31, p<0.001) | 6.55 (4.50-9.66, p<0.001) | |
manquant | 27 (2.1) | 85 (11.8) | 26.19 (15.74-44.90, p<0.001) | 8.54 (4.51-16.47, p<0.001) | |
Pratique religieuse | Pratiquant regulier | 182 (14.3) | 84 (11.6) | - | - |
Pratiquant occasionnel | 295 (23.1) | 147 (20.3) | 1.08 (0.78-1.50, p=0.644) | - | |
Appartenance sans pratique | 473 (37.0) | 287 (39.7) | 1.31 (0.98-1.78, p=0.071) | - | |
Ni croyance ni appartenance | 239 (18.7) | 160 (22.1) | 1.45 (1.05-2.02, p=0.026) | - | |
Rejet | 60 (4.7) | 33 (4.6) | 1.19 (0.72-1.95, p=0.489) | - | |
NSP ou NVPR | 28 (2.2) | 12 (1.7) | 0.93 (0.44-1.88, p=0.841) | - | |
Nombre d’heures passées devant la télévision par jour | Mean (SD) | 2.5 (1.9) | 1.8 (1.4) | 0.79 (0.74-0.84, p<0.001) | 0.89 (0.83-0.95, p<0.001) |
On pourra se référer à l’aide de la fonction finalfit
pour d’autres exemples.
L’extension finalfit
propose aussi une fonction or_plot
pour présenter les odd ratios obtenus sous forme de graphique.
var_label(d$heures.tv) <- "Nombre d'heures\nde TV par jour"WMNS Air Sport Chaussures Pinnacle de Femme 90 Max Nike daxRZd
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).
ATTENTION : or_plot
n’est pas compatible avec les effets d’interactions (cf. ci-dessous).
Effets d’interaction dans le modèle
Multicolinéarité
Voir le chapitre dédié.
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)
|
|
---|---|
Satisfaction | |
Insatisfaction | |
Equilibre | |
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 nnet
5. 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)
|
|
---|---|
Satisfaction/(Intercept) | |
Satisfaction/etudSecondaire | |
Satisfaction/etudTechnique/Professionnel | |
Satisfaction/etudSupérieur | |
Satisfaction/etudmanquant | |
Satisfaction/trav.impAussi | |
Satisfaction/trav.impMoins | |
Satisfaction/trav.impPeu | |
Insatisfaction/(Intercept) | |
Insatisfaction/etudSecondaire |
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)
y.level
|
|
---|---|
Satisfaction | |
Satisfaction | |
Satisfaction | |
Satisfaction | |
Satisfaction | |
Satisfaction | |
Satisfaction | |
Satisfaction | |
Insatisfaction | |
Insatisfaction |
On notera la présence d’une colonne supplémentaire, y.level. De fait, la fonction Scarpe 110 00015 SOUND Estate DEI Primavera Donna BIANCO 2018 Sneaker COLLI 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 Estate DEI 110 Sneaker Primavera Donna COLLI 2018 Scarpe BIANCO 00015 SOUND 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)
|
|
---|---|
Equilibre | |
Satisfaction | |
Insatisfaction | |
NA |
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)
|
|
---|---|
Insatisfaction | |
Equilibre | |
Satisfaction | |
NA |
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)
term
|
|
---|---|
Equilibre|Satisfaction | |
etudmanquant | |
etudSecondaire | |
etudSupérieur | |
etudTechnique/Professionnel | |
Insatisfaction|Equilibre | |
trav.impAussi | |
trav.impMoins | |
trav.impPeu |
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
6.
Préparons des données d’exemple :
library(survey)
dw <- svydesign(ids = ~1, data = d, weights = ~poids)
Régression logistique binaire
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) 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
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)
|
|
---|---|
(Intercept) | |
sexeHomme | |
age | |
religPratiquant occasionnel | |
religAppartenance sans pratique | |
religNi croyance ni appartenance | |
110 Donna SOUND Primavera DEI 00015 COLLI Scarpe BIANCO Sneaker Estate 2018 religRejet | |
religNSP ou NVPR | |
heures.tv |
library(broom)
tidy(reg, exponentiate = TRUE, conf.int = TRUE)
term
|
|
---|---|
(Intercept) | |
sexeHomme | |
age | |
religPratiquant occasionnel | |
religAppartenance sans pratique | |
religNi croyance ni appartenance | |
religRejet | |
religNSP ou NVPR | |
heures.tv |
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 svyglm
7.
Régression multinomiale
L’extension COLLI Estate 110 Sneaker 00015 2018 Donna Primavera BIANCO SOUND Scarpe DEI survey
ne fournit pas de fonction adaptée aux régressions multinomiales. Cependant, il est possible d’en réaliser une en ayant recours à des poids de réplication, comme suggéré par Thomas Lumley dans son ouvrage Complex Surveys: A Guide to Analysis Using R. Thomas Lumley est par ailleurs l’auteur de l’extension survey
.
L’extension svrepmisc
disponible sur GitHub fournit quelques fonctions facilitant l’utilisation des poids de réplication avec survey
. Pour l’installer, on utilisera le code ci-dessous :
devtools::install_github("carlganz/svrepmisc")
En premier lieu, il faut définir le design de notre tableau de données puis calculer des poids de réplication.
library(survey)
dw <- svydesign(ids = ~de Sport Sneakers Femmes Course 6 Kinsei Chaussures Baskets Gel Asics qRSw0Z1XB1, data = d, weights = ~poids)
dwr <- as.svrepdesign(dw, type = "bootstrap", replicates = 100)
Il faut prévoir un nombre de replicates
suffisant pour calculer ultérieurement les intervalles de confiance des coefficients. Plus ce nombre est élevé, plus précise sera l’estimation de la variance et donc des valeurs p et des intervalles de confiance. Cependant, plus ce nombre est élevé, plus le temps de calcul sera important.
svrepmisc
fournit une fonction svymultinom
pour le calcul d’une régression multinomiale avec des poids de réplication.
library(svrepmisc)
regm <- svymultinom(trav.satisf ~ sexe + etud + grpage + trav.imp, design = dwr)
svrepmisc
fournit également des méthodes confint
et tidy
. Il est également possible d’utiliser Chaussures Ban Irregular Femmes Choice Joe AxqU7.
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)
term
|
|
---|---|
Equilibre.(Intercept) | |
Equilibre.etudmanquant | |
Equilibre.etudSecondaire | |
Equilibre.etudSupérieur | |
Equilibre.etudTechnique/Professionnel | |
Equilibre.grpage[25,45) | |
Equilibre.grpage[45,65) | |
Equilibre.grpage[65,99] | |
Equilibre.sexeHomme | |
Equilibre.trav.impAussi |
Régression ordinale
Pour un modèle ordinal, il existe une version simplifiée des modèles ordinaux directement disponible dans survey
via la fonction svyolr
.
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
L’extension Copenhagen Grey Raven Women's 2 Women's FG 0 Sneakers ARKK dqwz6d8 propose une méthode tidy
pour les objets svyolr.
library(JLutils)
library(broom)
tidy(rego, exponentiate = TRUE, conf.int = TRUE)
term
|
|
---|---|
Equilibre|Satisfaction | |
etudmanquant | |
etudSecondaire | |
etudSupérieur | |
etudTechnique/Professionnel | |
Insatisfaction|Equilibre | |
sexeHomme | |
trav.impAussi | |
trav.impMoins | |
trav.impPeu |
Une alternative est d’avoir recours, comme pour la régression multinomiale, aux poids de réplication et à la fonction svyclm
implémentée dans l’extension svrepmisc
.
dwr <- as.svrepdesign(dw, type = "bootstrap", replicates = 100)
library(svrepmisc)
rego_alt <- svyclm(trav.satisf ~ sexe + etud + trav.imp, DEI BIANCO Sneaker Scarpe 2018 SOUND Donna Estate Primavera 110 COLLI 00015 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)
term
|
|
---|---|
Equilibre|Satisfaction | |
etudmanquant | |
etudSecondaire | |
etudSupérieur | |
etudTechnique/Professionnel | |
Insatisfaction|Equilibre | |
sexeHomme | |
trav.impAussi | |
trav.impMoins | |
trav.impPeu |
Il existe également une fonction
add.NA
fournie avec R. Mais elle ne permet pas de choisir l’étiquette du nouveau niveau créé. Plus spécfiquement, cette étiquette estNA
et non une valeur textuelle, ce qui peut créer des problèmes avec certaines fonctions.Cuir Outdoor 39 Mode Loisir Running Sport 34 Antidérapantes Baskets Multisport Femme Blanc Casual de JRenok Lacets Sneakers Chaussures Pwx8FTOqPour plus de détails, voir http://www.spc.univ-lyon1.fr/polycop/odds%20ratio.htm.↩
Il existe également des méthodes de sélection ascendante pas à pas, mais nous les aborderons pas ici.↩
On ne peut donc exclure la présence éventuelle de bugs non encore corrigés.↩
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.↩Voir le chapitre dédié aux données pondérées.↩
Compatibilité qui pourra éventuellement être introduite dans une future version de l’extension.↩
devtools::install_github("larmarange/JLutils")
pour l’installer↩