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().
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)
Warning: le package 'gtsummary' a été compilé avec la version R 4.4.3
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.
Warning in svymean.survey.design2(data, design[byfactor %in% byfactor[i], :
Sample size greater than population size: are weights correctly scaled?
Warning in svymean.survey.design2(data, design[byfactor %in% byfactor[i], :
Sample size greater than population size: are weights correctly scaled?
Warning in svymean.survey.design2(data, design[byfactor %in% byfactor[i], :
Sample size greater than population size: are weights correctly scaled?
Warning in svymean.survey.design2(data, design[byfactor %in% byfactor[i], :
Sample size greater than population size: are weights correctly scaled?
Warning in svymean.survey.design2(data, design[byfactor %in% byfactor[i], :
Sample size greater than population size: are weights correctly scaled?
Caractéristique
urbain N = 10261
rural N = 10021
Overall N = 20271
Â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%)
1095 (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%)
1351 (67%)
Manquant
5
1
6
1 Médiane (Q1 – Q3); 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.
Warning in svymean.survey.design2(data, design[byfactor %in% byfactor[i], :
Sample size greater than population size: are weights correctly scaled?
Warning in svymean.survey.design2(data, design[byfactor %in% byfactor[i], :
Sample size greater than population size: are weights correctly scaled?
Warning in svymean.survey.design2(data, design[byfactor %in% byfactor[i], :
Sample size greater than population size: are weights correctly scaled?
Warning in svymean.survey.design2(data, design[byfactor %in% byfactor[i], :
Sample size greater than population size: are weights correctly scaled?
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
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().
# On peut forcer l'utilisation des effectifs non pondérés à certains endroits via une personnalisation du thème, pour éviter d'avoir à le faire tableau par tableau.# NB : tbl_svysummary-str:missing_stat n'est actuellement disponible que dans la version de développement de gtsummary. cf. https://github.com/ddsjoberg/gtsummary/issues/2176list("tbl_svysummary-str:header-noby"="{N_unweighted} obs.","tbl_svysummary-str:header-withby"="**{level}** ({n_unweighted} obs.)",#"tbl_svysummary-str:missing_stat" = "{N_miss_unweighted} obs.","tbl_svysummary-arg:statistic"=list(all_continuous() ~"{median} ({p25} - {p75})",all_categorical() ~"{p}% ({n_unweighted} obs.)" ), "tbl_svysummary-arg:digits"=list(all_continuous() ~1,all_categorical() ~c(p =1, n_unweighted =0) )) |>set_gtsummary_theme()dp |>tbl_svysummary(by = milieu,include =c(educ, travail),missing_stat ="{N_miss_unweighted} obs.") |>add_overall(last =TRUE) |>bold_labels()
Warning in svymean.survey.design2(data, design[byfactor %in% byfactor[i], :
Sample size greater than population size: are weights correctly scaled?
Warning in svymean.survey.design2(data, design[byfactor %in% byfactor[i], :
Sample size greater than population size: are weights correctly scaled?
Warning in svymean.survey.design2(data, design[byfactor %in% byfactor[i], :
Sample size greater than population size: are weights correctly scaled?
Warning in svymean.survey.design2(data, design[byfactor %in% byfactor[i], :
Sample size greater than population size: are weights correctly scaled?
Caractéristique
urbain (912 obs.)1
rural (1088 obs.)1
Overall N = 20271
Niveau d'éducation
aucun
40,3% (375 obs.)
68,0% (763 obs.)
54,0% (1138 obs.)
primaire
24,4% (213 obs.)
25,6% (247 obs.)
25,0% (460 obs.)
secondaire
29,5% (275 obs.)
6,1% (73 obs.)
18,0% (348 obs.)
supérieur
5,7% (49 obs.)
0,3% (5 obs.)
3,0% (54 obs.)
A un emploi ?
non
39,2% (370 obs.)
26,9% (296 obs.)
33,1% (666 obs.)
oui
60,8% (537 obs.)
73,1% (790 obs.)
66,9% (1327 obs.)
Manquant
5 obs.
2 obs.
7 obs.
1
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 :
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.
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() :
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.
Âge révolu (en années) à la date de passation du questionnaire
28,5 (9,3)
28, 29
Région de résidence
Nord
30,1% (707 obs.)
28%, 33%
Est
8,7% (324 obs.)
7,7%, 9,8%
Sud
16,2% (407 obs.)
14%, 18%
Ouest
44,9% (562 obs.)
42%, 48%
Abréviation: IC = intervalle de confiance
1Moyenne (ET); % (n (unweighted) obs.)
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()
Warning in svymean.survey.design2(data, design[byfactor %in% byfactor[i], :
Sample size greater than population size: are weights correctly scaled?
Warning in svymean.survey.design2(data, design[byfactor %in% byfactor[i], :
Sample size greater than population size: are weights correctly scaled?
Warning in svymean.survey.design2(data, design[byfactor %in% byfactor[i], :
Sample size greater than population size: are weights correctly scaled?
Caractéristique
urbain (912 obs.)1
rural (1088 obs.)1
p-valeur2
Âge révolu (en années) à la date de passation du questionnaire
26,0 (20,0 - 33,0)
28,0 (22,0 - 36,0)
<0,001
Région de résidence
<0,001
Nord
25,9% (333 obs.)
34,5% (374 obs.)
Est
4,7% (109 obs.)
12,7% (215 obs.)
Sud
7,7% (111 obs.)
25,0% (296 obs.)
Ouest
61,7% (359 obs.)
27,8% (203 obs.)
1 Médiane (Q1 - Q3); % (n (unweighted) obs.)
2 Design-based KruskalWallis test; Pearson’s X^2: Rao & Scott adjustment
Table 30.3: Tests de comparaison avec prise en compte du plan d’échantillonnage
30.4 Calcul manuel avec {srvyr}
On peut avoir besoin de calculer des moyennes et/ou des proportions par sous-groupe, avec leurs intervalles de confiance, de manière manuelle, par exemple en amont d’un graphique à représenter avec ggplot2. Dans ce cas de figure, les fonctions natives de survey ne sont pas toujours très facile d’emploi et l’on pourra avantageusement recourir à srvyr pour bénéficier d’une syntaxe à la dplyr.
When `proportion` is unspecified, `survey_prop()` now defaults to `proportion = TRUE`.
ℹ This should improve confidence interval coverage.
This message is displayed once per session.
Là encore, on peut passer l’option vartpe = "ci" pour obtenir les intervalles de confiance3.
3 Par défaut, les intervalles de confiance sont calculés avec survey::svymean() et peuvent générer des valeurs inférieures à 0 ou supérieures à 1. Pour un calcul plus précis reposant sur survey::svyciprop(), on précisera proportion = TRUE. Plusieurs méthodes existent pour ce calcul, voir l’aide de survey::svyciprop().
Si l’on passe plusieurs variables dans le group_by(), les proportions sont calculées pour la dernière variable pour chaque combinaison des autres variables. Par exemple, pour la distribution du niveau d’éducation par milieu de résidence (i.e. la somme des proportions est de 100% pour le milieu urbain et de 100% pour celles du milieu rural, soit 200% au total) :
Si l’on souhaite les pourcentages que représentent chaque combinaison au sein de l’ensemble de l’échantillon (i.e. que la somme de toutes les proportions soit de 100%), on aura recours à srvyr::interact().
On peut utiliser, alternativement, la fonction guideR::proportion(), qui utilise de manière sous-jacente les fonctions de srvyr, mais dont la syntaxe est plus simple.
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.
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.
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.
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.