30  Analyses uni- et bivariées pondérées

30.1 La fonction tbl_svysummary()

Dans les chapitres sur la statistique univariée (cf. Chapitre 18) et la statistique bivariée (cf. Chapitre 19), nous avons abordé la fonction gtsummary::tbl_summary() qui permet de générer des tris à plats et des tableaux croisés prêts à être publiés.

Son équivalent pour les objets survey existe : il s’agit de la fonction gtsummary::tbl_svysummary() qui fonctionne de manière similaire.

Pour illustrer son fonctionnement, nous allons utiliser le jeu de données fecondite fournit dans le package questionr. Ce jeu de données fournit un tableau de données femmes comportant une variable poids de pondération que nous allons utiliser. Les données catégorielles étant stockées sous forme de vecteurs numériques avec étiquettes de valeurs (cf. Chapitre 12), nous allons les convertir en facteurs avec labelled::unlabelled(). De même, certaines valeurs manquantes sont indiquées sous formes de user NAs (cf. Section 13.2) : nous allons les convertir en valeurs manquantes classiques (regular NAs) avec labelled::user_na_to_na().

data("fecondite", package = "questionr")
library(srvyr)

Attachement du package : 'srvyr'
L'objet suivant est masqué depuis 'package:stats':

    filter
dp <- femmes |> 
  labelled::user_na_to_na() |> 
  labelled::unlabelled() |> 
  as_survey_design(weights = poids)
dp
Independent Sampling design (with replacement)
Called via srvyr
Sampling variables:
 - ids: `1`
 - weights: poids
Data variables: id_femme (dbl), id_menage (dbl), poids (dbl), date_entretien
  (date), date_naissance (date), age (dbl), milieu (fct), region (fct), educ
  (fct), travail (fct), matri (fct), religion (fct), journal (fct), radio
  (fct), tv (fct), nb_enf_ideal (dbl), test (fct)

Chargeons gtsummary et définissons le français comme langue de rendu des tableaux.

library(gtsummary)
theme_gtsummary_language(
  language = "fr", 
  decimal.mark = ",", 
  big.mark = " "
)
Setting theme `language: fr`

Pour réaliser un tableau croisé, il nous suffit d’appeler gtsummary::tbl_svysummary() de la même manière que l’on aurait procédé avec gtsummary::tbl_summary(). En arrière plan, gtsummary::tbl_svysummary() appellera les différentes fonctions statistiques de survey : la pondération ainsi que les spécificités du plan d’échantillonnage seront donc correctement prises en compte.

dp |> 
  tbl_svysummary(
    by = milieu,
    include = c(age, educ, travail)
  ) |> 
  add_overall(last = TRUE) |> 
  bold_labels()
Caractéristique urbain, N = 1 0261 rural, N = 1 0021 Total, N = 2 0271
Âge révolu (en années) à la date de passation du questionnaire 26 (20 – 33) 28 (22 – 36) 27 (21 – 35)
Niveau d'éducation
    aucun 414 (40%) 681 (68%) 1 095 (54%)
    primaire 251 (24%) 257 (26%) 507 (25%)
    secondaire 303 (30%) 61 (6,1%) 364 (18%)
    supérieur 58 (5,7%) 3 (0,3%) 61 (3,0%)
A un emploi ?
    non 401 (39%) 269 (27%) 670 (33%)
    oui 621 (61%) 731 (73%) 1 351 (67%)
    Manquant 5 1 6
1 Médiane (EI); n (%)
Table 30.1: Tableau croisé sur des données pondérées
Important

Par défaut, les effectifs (ainsi que les pourcentages et autres statistiques) affichés sont pondérés. Il est important de bien comprendre ce que représentent ces effectifs pondérés pour les interpréter correctement. Pour cela, il faut savoir comment les poids de l’enquête ont été calculés.

Dans certains cas, lorsque la population totale est connue, la somme des poids est égale à cette population totale dans laquelle l’échantillon a été tiré au sort. Les effectifs pondérés représentent donc une estimation des effectifs dans la population totale et ne représentent en rien le nombre d’observations dans l’enquête.

Dans d’autres enquêtes, les poids sont générés de telle manière que la somme des poids correspondent au nombre total de personnes enquêtées. Dans ce genre de situation, on a souvent tendance, à tort, à interpréter les effectifs pondérés comme un nombre d’observations. Or, il peut y avoir un écart important entre le nombre d’observations dans l’enquête et les effectifs pondérés.

On pourra éventuellement présenter séparément le nombre d’observations (i.e. les effectifs non pondérés) et les proportions pondérées. gtsummary::tbl_svysummary() fournit justement à la fois ces données pondérées et non pondérées. Il est vrai que cela nécessite quand même quelques manipulations. Pour les cellules, on précisera le type d’effectifs à afficher avec l’argument statistic. Pour personnaliser l’affiche du nombre de valeurs manquantes, cela doit se faire à un niveau plus global via gtsummary::set_gtsummary_theme(). Enfin, on passera par gtsummary::modify_header() pour personnaliser les en-têtes de colonne.

set_gtsummary_theme(
  list("tbl_summary-str:missing_stat" = "{N_miss_unweighted} obs.")
)
dp |> 
  tbl_svysummary(
    by = milieu,
    include = c(educ, travail),
    statistic = all_categorical() ~ "{p}% ({n_unweighted} obs.)",
    digits = all_categorical() ~ c(1, 0)
  ) |> 
  modify_header(
    all_stat_cols() ~ "**{level}** ({n_unweighted} obs.)"
  ) |> 
  bold_labels()
Caractéristique urbain (912 obs.)1 rural (1088 obs.)1
Niveau d'éducation
    aucun 40,3% (375 obs.) 68,0% (763 obs.)
    primaire 24,4% (213 obs.) 25,6% (247 obs.)
    secondaire 29,5% (275 obs.) 6,1% (73 obs.)
    supérieur 5,7% (49 obs.) 0,3% (5 obs.)
A un emploi ?
    non 39,2% (370 obs.) 26,9% (296 obs.)
    oui 60,8% (537 obs.) 73,1% (790 obs.)
    Manquant 5 obs. 2 obs.
1 % (n (unweighted) obs.)

Il faut noter qu’une modification du thème impactera tous les tableaux suivants, jusqu’à ce que le thème soit à nouveau modifié ou bien que l’on fasse appel à gtsummary::reset_gtsummary_theme().

30.2 Calcul manuel avec {survey}

Lorsque l’on travail avec un plan d’échantillonnage, on ne peut utiliser les fonctions statistiques classiques de R. On aura recours à leurs équivalents fournis par survey :

Ces fonctions prennent leurs arguments sous forme de formules pour spécifier les variables d’intérêt.

survey::svymean(~ age, dp)
      mean     SE
age 28.468 0.2697
survey::svymean(~ age, dp) |> confint()
       2.5 %   97.5 %
age 27.93931 28.99653
survey::svyquantile(~age, dp, quantile = c(0.25, 0.5, 0.75), ci = TRUE)
$age
     quantile ci.2.5 ci.97.5        se
0.25       21     21      22 0.2549523
0.5        27     27      28 0.2549523
0.75       35     35      37 0.5099045

attr(,"hasci")
[1] TRUE
attr(,"class")
[1] "newsvyquantile"

Les tris à plat se déclarent en passant comme argument le nom de la variable précédé d’un tilde (~), tandis que les tableaux croisés utilisent les noms des deux variables séparés par un signe plus (+) et précédés par un tilde (~)1.

  • 1 Cette syntaxe est similaire à celle de xtabs().

  • survey::svytable(~region, dp)
    region
        Nord      Est      Sud    Ouest 
    611.0924 175.7404 329.2220 911.2197 
    survey::svytable(~milieu + educ, dp)
            educ
    milieu        aucun   primaire secondaire  supérieur
      urbain 413.608780 250.665214 303.058978  58.412688
      rural  681.131096 256.694363  61.023980   2.679392

    La fonction questionr::freq() peut être utilisée si on lui passe en argument non pas la variable elle-même, mais son tri à plat obtenu avec survey::svytable() :

    survey::svytable(~region, dp) |> 
      questionr::freq(total = TRUE)
               n     %  val%
    Nord   611.1  30.1  30.1
    Est    175.7   8.7   8.7
    Sud    329.2  16.2  16.2
    Ouest  911.2  44.9  44.9
    Total 2027.3 100.0 100.0

    Les fonctions questionr::rprop() et questionr::cprop() peuvent être utilisées pour calculer les pourcentages en ligne ou en colonne.

    survey::svytable(~milieu + educ, dp) |> 
      questionr::cprop()
            educ
    milieu   aucun primaire secondaire supérieur Ensemble
      urbain  37.8  49.4     83.2       95.6      50.6   
      rural   62.2  50.6     16.8        4.4      49.4   
      Total  100.0 100.0    100.0      100.0     100.0   

    Le principe de la fonction survey::svyby() est similaire à celui de tapply() (cf. Section 19.2.3). Elle permet de calculer des statistiques selon plusieurs sous-groupes définis par un facteur.

    survey::svyby(~age, ~region, dp, survey::svymean)
          region      age        se
    Nord    Nord 29.03299 0.4753268
    Est      Est 27.54455 0.5261669
    Sud      Sud 28.96830 0.6148223
    Ouest  Ouest 28.08626 0.4458201

    30.3 Intervalles de confiance et tests statistiques

    La fonction gtsummary::add_ci() peut être appliquée à des tableaux produits avec gtsummary::tbl_svysummary()2. Les méthodes utilisées sont adaptées à la prise en compte d’un plan d’échantillonnage. On se référera à la document de la fonction pour plus de détails sur les méthodes statistiques utilisées. Rappel : pour les variables continues, on sera vigilant à ce que la statistique affichée (médiane par défaut) corresponde au type d’intervalle de confiance calculé (moyenne par défaut).

  • 2 Cela requiert une version récente (≥1.7.0) de gtsummary.

  • dp |> 
      tbl_svysummary(
        include = c(age, region),
        statistic = all_continuous() ~ "{mean} ({sd})"
      ) |> 
      add_ci() |> 
      bold_labels()
    Caractéristique N = 2 0271 95% CI2
    Âge révolu (en années) à la date de passation du questionnaire 28 (9) 28, 29
    Région de résidence
        Nord 611 (30%) 28%, 33%
        Est 176 (8,7%) 7,7%, 9,8%
        Sud 329 (16%) 14%, 18%
        Ouest 911 (45%) 42%, 48%
    1 Moyenne (ET); n (%)
    2 IC = intervalle de confiance
    Table 30.2: Intervalles de confiance avec prise en compte du plan d’échantillonnage

    De même, on peut aisément effectuer des tests de comparaison avec gtsummary::add_p(). Là aussi, les tests utilisés sont des adaptations des tests classiques avec différentes corrections pour tenir compte à la fois de la pondération et du plan d’échantillonnage.

    dp |> 
      tbl_svysummary(
        include = c(age, region),
        by = milieu
      ) |> 
      add_p() |> 
      bold_labels()
    Caractéristique urbain, N = 1 0261 rural, N = 1 0021 p-valeur2
    Âge révolu (en années) à la date de passation du questionnaire 26 (20 – 33) 28 (22 – 36) <0,001
    Région de résidence <0,001
        Nord 265 (26%) 346 (35%)
        Est 48 (4,7%) 128 (13%)
        Sud 79 (7,7%) 250 (25%)
        Ouest 633 (62%) 278 (28%)
    1 Médiane (EI); n (%)
    2 test de Wilcoxon sur la somme des rangs adapté aux plans d'échantillonnage complexes; test du Chi² avec la correction du second ordre de Rao & Scott
    Table 30.3: Tests de comparaison avec prise en compte du plan d’échantillonnage

    30.4 Impact du plan d’échantillonnage

    Lorsque l’on calcul des proportions, moyennes ou médianes pondérées, seuls les poids entrent en ligne de compte. Le plan d’échantillonnage (strates et/ou grappes) n’a de son côté pas d’effet. Par contre, le plan d’échantillonnage a un impact important sur le calcul des variances et, par extension, sur le calcul des intervalles de confiance et des tests de comparaison.

    Pour illustrer cela, nous allons considérer un même jeu de données, avec la même variable de poids, mais en faisant varier la présence de strates et de grappes.

    Commençons par regarder le jeu de données apistrat fourni par survey.

    data("api", package = "survey")
    nrow(apistrat)
    [1] 200
    summary(apistrat$pw)
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      15.10   19.05   32.28   30.97   44.21   44.21 
    sum(apistrat$pw)
    [1] 6194

    Nous avons ici un tableau de données de 200 lignes, avec des poids variant entre 15 et 44. Nous pouvons définir une pondération simple et croiser deux variables.

    d_ponderation_simple <- apistrat |> 
      as_survey_design(weights = pw)
    tbl <- survey::svytable(~ awards + yr.rnd, design = d_ponderation_simple)
    tbl
          yr.rnd
    awards      No     Yes
       No  2068.34  168.09
       Yes 3274.06  683.51

    Réalisons un test du Chi² entre ces deux variables. Si nous appliquions la fonction classique chisq.test() sur ce tableau, cette fonction considérait que nous avons 6194 observations (somme des poids) et dès lors nous obtiendrions une p-valeur très faible.

    tbl |> chisq.test()
    
        Pearson's Chi-squared test with Yates' continuity correction
    
    data:  tbl
    X-squared = 113.84, df = 1, p-value < 2.2e-16

    Le calcul précédent ne tient pas compte que nous n’avons que 200 observations dans notre échantillon. Refaisons le calcul survey::svychisq() qui est adaptée aux plans d’échantillonnage.

    survey::svychisq(~ awards + yr.rnd, design = d_ponderation_simple)
    
        Pearson's X^2: Rao & Scott adjustment
    
    data:  NextMethod()
    F = 2.9162, ndf = 1, ddf = 199, p-value = 0.08926

    Le résultat est ici tout autre et notre test n’est plus significatif au seuil de 5% ! Ici, les corrections de Rao & Scott permettent justement de tenir compte que nous avons un échantillon de seulement 200 observations.

    Regardons maintenant si, à poids égal, il y a une différence entre une enquête stratifiée et une enquête en grappes.

    # Pondération simple
    survey::svytable(~ awards + yr.rnd, design = d_ponderation_simple)
          yr.rnd
    awards      No     Yes
       No  2068.34  168.09
       Yes 3274.06  683.51
    survey::svychisq(~ awards + yr.rnd, design = d_ponderation_simple)
    
        Pearson's X^2: Rao & Scott adjustment
    
    data:  NextMethod()
    F = 2.9162, ndf = 1, ddf = 199, p-value = 0.08926
    # Enquête stratifiée
    d_strates <- apistrat |> 
      as_survey_design(weights = pw, strata = stype)
    survey::svytable(~ awards + yr.rnd, design = d_strates)
          yr.rnd
    awards      No     Yes
       No  2068.34  168.09
       Yes 3274.06  683.51
    survey::svychisq(~ awards + yr.rnd, design = d_strates)
    
        Pearson's X^2: Rao & Scott adjustment
    
    data:  NextMethod()
    F = 2.9007, ndf = 1, ddf = 197, p-value = 0.09012
    # Enquête en grappes
    d_grappes <- apistrat |> 
      as_survey_design(weights = pw, ids = dnum)
    survey::svytable(~ awards + yr.rnd, design = d_grappes)
          yr.rnd
    awards      No     Yes
       No  2068.34  168.09
       Yes 3274.06  683.51
    survey::svychisq(~ awards + yr.rnd, design = d_grappes)
    
        Pearson's X^2: Rao & Scott adjustment
    
    data:  NextMethod()
    F = 3.1393, ndf = 1, ddf = 134, p-value = 0.0787

    On le constate : dans les trois cas les tableaux croisés sont identiques, mais pour autant les trois p-valeurs diffèrent.

    Dès lors qu’un calcul de variance est impliqué, la simple prise en compte des poids est insuffisante : il faut appliquer des corrections en fonction du plan d’échantillonnage !

    Pas d’inquiétude, survey s’en occupe pour vous, dès lors que le plan d’échantillonnage a correctement été défini.