library(tidyverse)
library(labelled)
data(hdv2003, package = "questionr")
d <-
hdv2003 |>
mutate(
sexe = sexe |> fct_relevel("Femme"),
groupe_ages = age |>
cut(
c(18, 25, 45, 65, 99),
right = FALSE,
include.lowest = TRUE,
labels = c("18-24 ans", "25-44 ans",
"45-64 ans", "65 ans et plus")
),
etudes = nivetud |>
fct_recode(
"Primaire" = "N'a jamais fait d'etudes",
"Primaire" = "A arrete ses etudes, avant la derniere annee d'etudes primaires",
"Primaire" = "Derniere annee d'etudes primaires",
"Secondaire" = "1er cycle",
"Secondaire" = "2eme cycle",
"Technique / Professionnel" = "Enseignement technique ou professionnel court",
"Technique / Professionnel" = "Enseignement technique ou professionnel long",
"Supérieur" = "Enseignement superieur y compris technique superieur"
) |>
fct_na_value_to_level("Non documenté")
) |>
set_variable_labels(
sport = "Pratique un sport ?",
sexe = "Sexe",
groupe_ages = "Groupe d'âges",
etudes = "Niveau d'études",
relig = "Rapport à la religion",
heures.tv = "Heures de télévision / jour",
poids = "Pondération de l'enquête"
)
30 Régression logistique binaire pondérée
Nous avons abordé la régression logistique binaire non pondérée dans un chapitre dédié, cf. Chapitre 21. Elle se réalise classiquement avec la fonction glm()
en spécifiant family = binomial
.
Lorsque l’on utilise des données d’enquêtes, l’approche est similaire sauf que l’on aura recours à la fonction survey::svyglm()
qui sait gérer des objets survey : non seulement la pondération sera prise en compte, mais le calcul des intervalles de confiance et des p-valeurs sera ajusté en fonction du plan d’échantillonnage.
30.1 Données des exemples
Nous allons reprendre les même données issues de l’enquête Histoire de vie 2003, mais en tenant compte cette fois-ci des poids de pondération fourni dans la variable poids.
Il ne nous reste qu’à définir notre objet survey en spécifiant la pondération fournie avec l’enquête. La documentation ne mentionne ni strates ni grappes.
30.2 Calcul de la régression logistique binaire
La syntaxe de survey::svyglm()
est similaire à celle de glm()
sauf qu’elle a un argument design
au lieu de data
.
La plupart du temps, les poids de pondération ne sont pas des nombres entiers, mais des nombres décimaux. Dès lors, on ne peut plus utiliser la famille de modèles binomiale (qui repose sur des nombres entiers de succès et d’échecs)1. On aura plutôt recours à la famille quasi-binomiale, que l’on spécifie avec family = quasibinomial
et qui constitue une extension de la famille binomiale pouvant gérer des poids non entiers.
1 Si l’on indique family = binomial
, vous obtiendrez avec une version récente de R un message d’avertissement du type Avis : nombre de succès non entier dans un glm binomial !
. Avec une version plus ancienne de R, vous devriez même avoir un message d’erreur.
Simple, non ?
30.3 Sélection de modèle
Comme précédemment, il est possible de procéder à une sélection de modèle pas à pas, par minimisation de l’AIC, avec step()
.
Start: AIC=2309.89
sport ~ sexe + groupe_ages + etudes + relig + heures.tv
Df Deviance AIC
- relig 5 2266.3 2302.2
<none> 2263.9 2309.9
- heures.tv 1 2276.2 2320.2
- sexe 1 2276.4 2320.4
- groupe_ages 3 2313.9 2353.8
- etudes 4 2383.5 2421.2
Step: AIC=2296.28
sport ~ sexe + groupe_ages + etudes + heures.tv
Df Deviance AIC
<none> 2266.3 2296.3
- heures.tv 1 2278.4 2306.4
- sexe 1 2279.0 2307.0
- groupe_ages 3 2318.3 2342.1
- etudes 4 2387.2 2408.8
30.4 Affichage des résultats
Nous pouvons tout à fait utiliser gtsumarry::tbl_regression()
avec ce type de modèles. De même, on peut utiliser gtsummary::add_global_p()
pour calculer les p-valeurs globales des variables ou encore gtsummary::add_vif()
pour vérifier la multicolinéarité (cf. Chapitre 25).
mod2 |>
tbl_regression(exponentiate = TRUE) |>
add_global_p(keep = TRUE) |>
add_vif() |>
bold_labels()
Caractéristique | OR1 | 95% IC1 | p-valeur | GVIF1 | Adjusted GVIF2,1 |
---|---|---|---|---|---|
Sexe | 0,005 | 1,0 | 1,0 | ||
Femme | — | — | |||
Homme | 1,44 | 1,12 – 1,87 | 0,005 | ||
Groupe d'âges | <0,001 | 2,1 | 1,1 | ||
18-24 ans | — | — | |||
25-44 ans | 0,85 | 0,48 – 1,51 | 0,6 | ||
45-64 ans | 0,40 | 0,22 – 0,73 | 0,003 | ||
65 ans et plus | 0,37 | 0,19 – 0,72 | 0,004 | ||
Niveau d'études | <0,001 | 2,2 | 1,1 | ||
Primaire | — | — | |||
Secondaire | 2,66 | 1,62 – 4,38 | <0,001 | ||
Technique / Professionnel | 3,09 | 1,90 – 5,00 | <0,001 | ||
Supérieur | 6,54 | 3,99 – 10,7 | <0,001 | ||
Non documenté | 10,3 | 4,60 – 23,0 | <0,001 | ||
Heures de télévision / jour | 0,89 | 0,82 – 0,97 | 0,006 | 1,1 | 1,0 |
1 OR = rapport de cotes, IC = intervalle de confiance, GVIF = Generalized Variance Inflation Factor | |||||
2 GVIF^[1/(2*df)] |
Pour un graphique des coefficients, nous pouvons utiliser ggstats::ggcoef_model()
.
30.5 Prédictions marginales
Pour visualiser les prédictions marginales moyennes du modèle (cf. Section 22.3), nous pouvons utiliser broom.helpers::plot_marginal_predictions()
.