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.

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

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.

library(srvyr)
library(survey)
dp <- d |> 
  as_survey_design(weights = poids)

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.

    • mod <- svyglm(
  sport ~ sexe + groupe_ages + etudes + relig + heures.tv,
  family = quasibinomial,
  design = dp
)

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

    mod2 <- step(mod)
    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).

    library(gtsummary)
theme_gtsummary_language("fr", decimal.mark = ",", big.mark = " ")
    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)]
    Table 30.1: Facteurs associés à la pratique d’un sport (régression logistique pondérée)

    Pour un graphique des coefficients, nous pouvons utiliser ggstats::ggcoef_model().

    mod2 |> 
  ggstats::ggcoef_model(exponentiate = TRUE)

    Figure 30.1: Facteurs associés à la pratique d’un sport (régression logistique pondérée)

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

    mod2 |> 
  broom.helpers::plot_marginal_predictions(type = "response") |> 
  patchwork::wrap_plots() &
  scale_y_continuous(
    limits = c(0, .8),
    labels = scales::label_percent()
  )

    Figure 30.2: Prédictions marginales moyennes du modèle pondéré